Skip to content

Commit

Permalink
Icarus: missing working precision and initialisation. (#78)
Browse files Browse the repository at this point in the history 
* Reals to wp.

* Changes after Steve's review.
  • Loading branch information
@alejandrobodas
alejandrobodas authored Jul 21, 2023
1 parent 49c0d7e commit e0558d4
Showing 1 changed file with 22 additions and 22 deletions.
44 changes: 22 additions & 22 deletions src/simulator/icarus/icarus.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ MODULE MOD_ICARUS
isccp_taumin = 0.3_wp, & ! Minimum optical depth for joint-hostogram
pstd = 1013250._wp, & ! Mean sea-level pressure (Pa)
isccp_t0 = 296._wp, & ! Mean surface temperature (K)
output_missing_value = -1.E+30 ! Missing values
output_missing_value = -1.E+30_wp ! Missing values

contains
! ##########################################################################
Expand Down Expand Up @@ -231,8 +231,8 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
itrop(1:npoints) = 1

do ilev=1,nlev
where(pfull(1:npoints,ilev) .lt. 40000. .and. &
pfull(1:npoints,ilev) .gt. 5000. .and. &
where(pfull(1:npoints,ilev) .lt. 40000._wp .and. &
pfull(1:npoints,ilev) .gt. 5000._wp .and. &
at(1:npoints,ilev) .lt. attropmin(1:npoints))
ptrop(1:npoints) = pfull(1:npoints,ilev)
attropmin(1:npoints) = at(1:npoints,ilev)
Expand All @@ -257,16 +257,16 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
! ############################################################################
do ilev=1,nlev
press(1:npoints) = pfull(1:npoints,ilev)*10._wp
dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10
dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10._wp
atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp)
rvh20(1:npoints) = qv(1:npoints,ilev)*amd/amw
wk(1:npoints) = rvh20(1:npoints)*avo*atmden/amd
rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev))
rh20s(1:npoints) = rvh20(1:npoints)*rhoave(1:npoints)
rfrgn(1:npoints) = rhoave(1:npoints)-rh20s(1:npoints)
tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0))
tauwv(1:npoints) = wk(1:npoints)*1.e-20*((0.0224697_wp*rh20s(1:npoints)* &
tmpexp(1:npoints))+(3.41817e-7*rfrgn(1:npoints)))*0.98_wp
tauwv(1:npoints) = wk(1:npoints)*1.e-20_wp*((0.0224697_wp*rh20s(1:npoints)* &
tmpexp(1:npoints))+(3.41817e-7_wp*rfrgn(1:npoints)))*0.98_wp
dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints))
enddo

Expand All @@ -282,7 +282,7 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,

! Update trans_layers_above with transmissivity from this layer for next time around loop
trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*&
(1.-dem_wv(1:npoints,ilev))
(1._wp-dem_wv(1:npoints,ilev))
enddo

! Add in surface emission
Expand All @@ -307,7 +307,7 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
! Emissivity
dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
demIN(1:npoints,ibox,ilev) .eq. 0)
demIN(1:npoints,ibox,ilev) .eq. 0._wp)

! Increase TOA flux emitted from layer
fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)*bb*trans_layers_above(1:npoints,ibox)
Expand Down Expand Up @@ -349,20 +349,20 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
taumin(1:npoints) = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
if (isccp_top_height .eq. 1) then
do j=1,npoints
if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then
if (transmax(j) .gt. 0.001_wp .and. transmax(j) .le. 0.9999999_wp) then
fluxtopinit(j) = fluxtop(j,ibox)
tauir(j) = tau(j,ibox)/2.13_wp
endif
enddo
do icycle=1,2
do j=1,npoints
if (tau(j,ibox) .gt. (tauchk)) then
if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then
if (transmax(j) .gt. 0.001_wp .and. transmax(j) .le. 0.9999999_wp) then
emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j) )
fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
if (tb(j,ibox) .gt. 260.) then
if (tb(j,ibox) .gt. 260._wp) then
tauir(j) = tau(j,ibox) / 2.56_wp
end if
end if
Expand All @@ -373,7 +373,7 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,

! Cloud-top temperature
where(tau(1:npoints,ibox) .gt. tauchk)
tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox))))
tb(1:npoints,ibox)= 1307.27_wp/ (log(1._wp + (1._wp/fluxtop(1:npoints,ibox))))
where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints))
tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp
tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
Expand All @@ -396,6 +396,8 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
! top pressure (isccp_top_height = 2) or the radiatively determined cloud top
! pressure (isccp_top_height = 1 or 3)
! ####################################################################################
ptop(:,:)=0._wp
levmatch(:,:)=0
do ibox=1,ncol
!segregate according to optical thickness
if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
Expand Down Expand Up @@ -438,17 +440,16 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
end if
enddo
else
ptop(1:npoints,ibox)=0.
do ilev=1,nlev
where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0))
where((ptop(1:npoints,ibox) .eq. 0._wp) .and.(frac_out(1:npoints,ibox,ilev) > 0._wp))
ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
levmatch(1:npoints,ibox)=ilev
endwhere
end do
end if
where(tau(1:npoints,ibox) .le. tauchk)
ptop(1:npoints,ibox)=0._wp
levmatch(1:npoints,ibox)=0._wp
levmatch(1:npoints,ibox)=0
endwhere
enddo

Expand All @@ -459,7 +460,7 @@ SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,
boxptop(1:npoints,1:ncol) = output_missing_value
do ibox=1,ncol
do j=1,npoints
if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0.) then
if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0._wp) then
if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
boxtau(j,ibox) = tau(j,ibox)
boxptop(j,ibox) = ptop(j,ibox)!/100._wp
Expand Down Expand Up @@ -535,7 +536,7 @@ SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,
do ilev2=1,7
do j=1,npoints !
if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
fq_isccp(j,ilev,ilev2)= 0.
fq_isccp(j,ilev,ilev2)= 0._wp
else
fq_isccp(j,ilev,ilev2)= output_missing_value
end if
Expand All @@ -560,7 +561,7 @@ SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,
! Compute column quantities and joint-histogram
do j=1,npoints
! Subcolumns that are cloudy(true) and not(false)
box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0.)
box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0._wp)

! Compute joint histogram and column quantities for points that are sunlit and cloudy
if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then
Expand All @@ -571,7 +572,7 @@ SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,
fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol

! Column cloud area
totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol
totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin),wp)/ncol

! Subcolumn cloud albedo
!albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),&
Expand All @@ -591,10 +592,10 @@ SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,
enddo

! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
where(totalcldarea(1:npoints) .gt. 0)
where(totalcldarea(1:npoints) .gt. 0._wp)
meanptop(1:npoints) = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
meantaucld(1:npoints) = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp)
meantaucld(1:npoints) = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1._wp))**(1._wp/0.895_wp)
elsewhere
meanptop(1:nPoints) = output_missing_value
meanalbedocld(1:nPoints) = output_missing_value
Expand Down Expand Up @@ -642,4 +643,3 @@ subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
enddo
end subroutine cosp_simulator_optics
end module MOD_ICARUS

0 comments on commit e0558d4

Please sign in to comment.