Skip to content

Commit

Permalink
Cleaned up a tad. Added some diagnostics for debuggind in SCM.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dustinswales
dustinswales committed Dec 3, 2019
1 parent d2799f4 commit 2752142
Show file tree
Hide file tree
Showing 6 changed files with 104 additions and 111 deletions.
17 changes: 2 additions & 15 deletions physics/GFS_rrtmgp_lw_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,28 +71,15 @@ subroutine GFS_rrtmgp_lw_pre_run (Model, Grid, Sfcprop, Statein, ncol, p_lay,
errflg ! Error flag

! Local
integer :: iSFC, iTOA
logical :: top_at_1
real(kind_phys), dimension(ncol, Model%levs, Model%rrtmgp_nBandsSW, NF_AESW) :: &
aerosolssw2

! Initialize CCPP error handling variables
errmsg = ''
errflg = 0

if (.not. Model%lslwr) return

! #######################################################################################
! What is vertical ordering?
! #######################################################################################
top_at_1 = (Statein%prsi(1,1) .lt. Statein%prsi(1, Model%levs))
if (top_at_1) then
iSFC = Model%levs
iTOA = 1
else
iSFC = 1
iTOA = Model%levs
endif

! #######################################################################################
! Call module_radiation_surface::setemis(),to setup surface emissivity for LW radiation.
! #######################################################################################
Expand All @@ -106,7 +93,7 @@ subroutine GFS_rrtmgp_lw_pre_run (Model, Grid, Sfcprop, Statein, ncol, p_lay,
! #######################################################################################
! Call module_radiation_aerosols::setaer(),to setup aerosols property profile
! #######################################################################################
call setaer(p_lev, p_lay, Statein%prslk(1:NCOL,iSFC:iTOA), tv_lay, relhum, &
call setaer(p_lev, p_lay, Statein%prslk(1:NCOL,:), tv_lay, relhum, &
Sfcprop%slmsk, tracer, Grid%xlon, Grid%xlat, ncol, Model%levs, Model%levs+1, &
.true., Model%lslwr, aerosolssw2, aerosolslw, aerodp)

Expand Down
53 changes: 29 additions & 24 deletions physics/GFS_rrtmgp_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,38 +213,43 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
! #######################################################################################
! Compute some fields needed by RRTMGP
! #######################################################################################


! Water-vapor mixing-ratio
q_lay(1:ncol,:) = max( 1.e-6, Statein%qgrs(:,:,1))

! Pressure at layer-interface
p_lev(1:NCOL,iSFC:iTOA+1) = Statein%prsi(1:NCOL,1:Model%levs+1)
!
p_lev(1:NCOL,:) = Statein%prsi(1:NCOL,:)

! Pressure at layer-center
p_lay(1:NCOL,iSFC:iTOA) = Statein%prsl(1:NCOL,1:Model%levs)
!
p_lay(1:NCOL,:) = Statein%prsl(1:NCOL,:)

! Temperature at layer-center
t_lay(1:NCOL,iSFC:iTOA) = Statein%tgrs(1:NCOL,1:Model%levs)
t_lay(1:NCOL,:) = Statein%tgrs(1:NCOL,:)

!
! Temperature at layer-interfaces
t_lev(1:NCOL,iSFC) = Sfcprop%tsfc(1:NCOL)
do iCol=1,NCOL
do iLay=iSFC+1,iTOA
t_lev(iCol,iLay) = (t_lay(iCol,iLay)+t_lay(iCol,iLay-1))/2._kind_phys
enddo
t_lev(iCol,iTOA+1) = t_lay(iCol,iTOA)
enddo

if (top_at_1) then
t_lev(1:NCOL,1) = t_lay(1:NCOL,iTOA)
t_lev(1:NCOL,2:iSFC) = (t_lay(1:NCOL,2:iSFC)+t_lay(1:NCOL,1:iSFC-1))/2._kind_phys
t_lev(1:NCOL,iSFC+1) = Sfcprop%tsfc(1:NCOL)
else
t_lev(1:NCOL,1) = Sfcprop%tsfc(1:NCOL)
t_lev(1:NCOL,2:iTOA) = (t_lay(1:NCOL,2:iTOA)+t_lay(1:NCOL,1:iTOA-1))/2._kind_phys
t_lev(1:NCOL,iTOA+1) = t_lay(1:NCOL,iTOA)
endif

! Compute layer pressure thicknes
deltaP = p_lev(:,iSFC:iTOA)-p_lev(:,iSFC+1:iTOA+1)
deltaP = abs(p_lev(:,2:model%levs+1)-p_lev(:,1:model%levs))

! Compute a bunch of thermodynamic fields needed by the macrophysics schemes. Relative humidity,
! saturation mixing-ratio, vapor mixing-ratio, virtual temperature, layer thickness,...
do iCol=1,NCOL
do iLay=iSFC,iTOA
es = min( Statein%prsl(iCol,iLay), fpvs( Statein%tgrs(iCol,iLay) ) ) ! fpvs and prsl in pa
qs = max( QMIN, eps * es / (Statein%prsl(iCol,iLay) + epsm1*es) )
relhum(iCol,iLay) = max( 0._kind_phys, min( 1._kind_phys, max(QMIN, Statein%qgrs(iCol,iLay,1))/qs ) )
do iLay=1,Model%levs
es = min( p_lay(iCol,iLay), fpvs( t_lay(iCol,iLay) ) ) ! fpvs and prsl in pa
qs = max( QMIN, eps * es / (p_lay(iCol,iLay) + epsm1*es) )
relhum(iCol,iLay) = max( 0._kind_phys, min( 1._kind_phys, max(QMIN, q_lay(iCol,iLay))/qs ) )
qs_lay(iCol,iLay) = qs
q_lay(iCol,iLay) = max( 1.e-6, Statein%qgrs(iCol,iLay,1) )
tv_lay(iCol,iLay) = Statein%tgrs(iCol,iLay) * (1._kind_phys + fvirt*q_lay(iCol,iLay))
tv_lay(iCol,iLay) = t_lay(iCol,iLay) * (1._kind_phys + fvirt*q_lay(iCol,iLay))
deltaZ(iCol,iLay) = (rog*0.001) * (log(p_lev(iCol,iLay)) - log(p_lev(iCol,iLay+1))) * tv_lay(iCol,iLay)
enddo
enddo
Expand All @@ -254,18 +259,18 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
! #######################################################################################
! First recast remaining all tracers (except sphum) forcing them all to be positive
do j = 2, model%NTRAC
tracer(1:NCOL,1:Model%levs,j) = max(0.0, Statein%qgrs(1:NCOL,1:Model%levs,j))
tracer(1:NCOL,:,j) = max(0.0, Statein%qgrs(1:NCOL,:,j))
enddo

if (Model%ntoz > 0) then
do iLay=iSFC,iTOA
do iLay=1,Model%levs
do iCol=1,NCOL
o3_lay(iCol,iLay) = max( QMIN, tracer(iCol,iLay,Model%ntoz) )
enddo
enddo
! OR Use climatological ozone data
else
call getozn (Statein%prslk(1:NCOL,iSFC:iTOA), Grid%xlat, NCOL, Model%levs, o3_lay)
call getozn (Statein%prslk(1:NCOL,:), Grid%xlat, NCOL, Model%levs, o3_lay)
endif

! #######################################################################################
Expand Down
44 changes: 29 additions & 15 deletions physics/GFS_rrtmgp_sw_post.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,38 +144,41 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou
! #######################################################################################
! Initialize
hswc = 0
Diag%topfsw = topfsw_type ( 0., 0., 0. )
! sfcflx_sw = sfcfsw_type ( 0., 0., 0., 0. )
Diag%topfsw = topfsw_type ( 0., 0., 0. )
Radtend%sfcfsw = sfcfsw_type ( 0., 0., 0., 0. )
if (l_clrskysw_hr) then
hsw0(:,:) = 0.
endif
if (l_fluxessw2D) then
flxprf_sw = profsw_type ( 0., 0., 0., 0. )
endif
!if (l_sfcfluxessw1D) then
! scmpsw = cmpfsw_type (0.,0.,0.,0.,0.,0.)
!endif
if (l_sfcfluxessw1D) then
scmpsw = cmpfsw_type (0.,0.,0.,0.,0.,0.)
endif

if (Model%lsswr .and. nDay .gt. 0) then
! Clear-sky heating-rate (optional)
if (l_clrskysw_HR) then
call check_error_msg('GFS_rrtmgp_post',compute_heating_rate( &
fluxswUP_clrsky(idxday,:), &
fluxswDOWN_clrsky(idxday,:), &
p_lev(idxday,1:Model%levs+1), &
p_lev(idxday,:), &
thetaTendClrSky))
hsw0(idxday,:)=thetaTendClrSky
endif
! All-sky heating-rate (mandatory)
call check_error_msg('GFS_rrtmgp_post',compute_heating_rate( &
fluxswUP_allsky(idxday,:), &
fluxswDOWN_allsky(idxday,:), &
p_lev(idxday,1:Model%levs+1), &
p_lev(idxday,:), &
thetaTendAllSky))
hswc(idxday,:) = thetaTendAllSky

! Copy fluxes from RRTGMP types into model radiation types.
! Mandatory outputs
write(*,"(a11,2i8)") "iTOA/iSFC: ",iTOA,iSFC
write(*,*) "fluxswDOWN_allsky: ",fluxswDOWN_allsky(idxday,:)
write(*,*) "fluxswDOWN_clrsky: ",fluxswDOWN_clrsky(idxday,:)
Diag%topfsw(idxday)%upfxc = fluxswUP_allsky(idxday,iTOA)
Diag%topfsw(idxday)%upfx0 = fluxswUP_clrsky(idxday,iTOA)
Diag%topfsw(idxday)%dnfxc = fluxswDOWN_allsky(idxday,iTOA)
Expand All @@ -202,7 +205,7 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou
do k = 1, Model%levs
Radtend%htrsw(1:im,k) = hswc(1:im,k)
enddo
! Clear-sk heating rate
! Clear-sky heating rate
if (Model%swhtr) then
do k = 1, Model%levs
Radtend%swhc(1:im,k) = hsw0(1:im,k)
Expand Down Expand Up @@ -272,23 +275,34 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou
if (Radtend%coszen(i) > 0.) then
! SW all-sky fluxes
tem0d = Model%fhswr * Radtend%coszdg(i) / Radtend%coszen(i)
Diag%fluxr(i,2 ) = Diag%fluxr(i,2) + fluxswUP_allsky( idxday(i),iTOA) * tem0d ! total sky top sw up
Diag%fluxr(i,3 ) = Diag%fluxr(i,3) + fluxswUP_allsky( idxday(i),iSFC) * tem0d ! total sky sfc sw up
Diag%fluxr(i,4 ) = Diag%fluxr(i,4) + fluxswDOWN_allsky(idxday(i),iSFC) * tem0d ! total sky sfc sw dn
!write(*,"(a23,3f10.6)") 'In GFS_rrtmgp_sw_post: ',Diag%topfsw(i)%dnfxc, tem0d,Diag%fluxr(i,23)
!write(*,"(a23,f20.15)") 'In GFS_rrtmgp_sw_post: ',Model%fhswr
!Diagfluxr(i,2 ) = Diag%fluxr(i,2) + fluxswUP_allsky( i,iTOA) * tem0d ! total sky top sw up
!Diag%fluxr(i,3 ) = Diag%fluxr(i,3) + fluxswUP_allsky( i,iSFC) * tem0d ! total sky sfc sw up
!Diag%fluxr(i,4 ) = Diag%fluxr(i,4) + fluxswDOWN_allsky(i,iSFC) * tem0d ! total sky sfc sw dn
Diag%fluxr(i,2 ) = Diag%fluxr(i,2) + Diag%topfsw(i)%upfxc * tem0d ! total sky top sw up
Diag%fluxr(i,3 ) = Diag%fluxr(i,3) + Radtend%sfcfsw(i)%upfxc * tem0d
Diag%fluxr(i,4 ) = Diag%fluxr(i,4) + Radtend%sfcfsw(i)%dnfxc * tem0d ! total sky sfc sw dn
! SW uv-b fluxes
Diag%fluxr(i,21) = Diag%fluxr(i,21) + scmpsw(i)%uvbfc * tem0d ! total sky uv-b sw dn
Diag%fluxr(i,22) = Diag%fluxr(i,22) + scmpsw(i)%uvbf0 * tem0d ! clear sky uv-b sw dn
! SW TOA incoming fluxes
Diag%fluxr(i,23) = Diag%fluxr(i,23) + fluxswDOWN_allsky(idxday(i),iTOA) * tem0d ! top sw dn
!temiag%fluxr(i,23) = Diag%fluxr(i,23) + fluxswDOWN_allsky(i,iTOA) * tem0d ! top sw dn
Diag%fluxr(i,23) = Diag%fluxr(i,23) + Diag%topfsw(i)%dnfxc * tem0d ! top sw dn
write(*,"(a23,3f10.6)") 'In GFS_rrtmgp_sw_post: ',Diag%topfsw(i)%dnfxc, tem0d,Diag%fluxr(i,23)
! SW SFC flux components
Diag%fluxr(i,24) = Diag%fluxr(i,24) + scmpsw(i)%visbm * tem0d ! uv/vis beam sw dn
Diag%fluxr(i,25) = Diag%fluxr(i,25) + scmpsw(i)%visdf * tem0d ! uv/vis diff sw dn
Diag%fluxr(i,26) = Diag%fluxr(i,26) + scmpsw(i)%nirbm * tem0d ! nir beam sw dn
Diag%fluxr(i,27) = Diag%fluxr(i,27) + scmpsw(i)%nirdf * tem0d ! nir diff sw dn
! SW clear-sky fluxes
Diag%fluxr(i,29) = Diag%fluxr(i,29) + fluxswUP_clrsky( idxday(i),iTOA) * tem0d ! clear sky top sw up
Diag%fluxr(i,31) = Diag%fluxr(i,31) + fluxswUP_clrsky( idxday(i),iSFC) * tem0d ! clear sky sfc sw up
Diag%fluxr(i,32) = Diag%fluxr(i,32) + fluxswDOWN_clrsky(idxday(i),iSFC) * tem0d ! clear sky sfc sw dn
!Diag%fluxr(i,29) = Diag%fluxr(i,29) + fluxswUP_clrsky( i,iTOA) * tem0d ! clear sky top sw up
!Diag%fluxr(i,31) = Diag%fluxr(i,31) + fluxswUP_clrsky( i,iSFC) * tem0d ! clear sky sfc sw up
!Diag%fluxr(i,32) = Diag%fluxr(i,32) + fluxswDOWN_clrsky(i,iSFC) * tem0d ! clear sky sfc sw dn
Diag%fluxr(i,29) = Diag%fluxr(i,29) + Diag%topfsw(i)%upfx0 * tem0d
Diag%fluxr(i,31) = Diag%fluxr(i,31) + Radtend%sfcfsw(i)%upfx0 * tem0d
Diag%fluxr(i,32) = Diag%fluxr(i,32) + Radtend%sfcfsw(i)%dnfx0 * tem0d

endif
enddo

Expand Down
31 changes: 9 additions & 22 deletions physics/GFS_rrtmgp_sw_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,31 +87,18 @@ subroutine GFS_rrtmgp_sw_pre_run (Model, Interstitial, Grid, Sfcprop, Statein,
errflg ! Error flag

! Local variables
integer :: i, j, iCol, iBand, iSFC, iTOA, iLay
integer :: i, j, iCol, iBand, iLay
real(kind_phys), dimension(ncol, NF_ALBD) :: sfcalb
real(kind_phys), dimension(ncol, Model%levs, sw_gas_props%get_nband(), NF_AESW) :: &
aerosolssw2
real(kind_phys), dimension(ncol, Model%levs, Model%rrtmgp_nBandsLW, NF_AELW) :: &
aerosolslw
logical :: top_at_1

! Initialize CCPP error handling variables
errmsg = ''
errflg = 0

if (.not. Model%lsswr) return

! #######################################################################################
! What is vertical ordering?
! #######################################################################################
top_at_1 = (Statein%prsi(1,1) .lt. Statein%prsi(1, Model%levs))
if (top_at_1) then
iSFC = Model%levs
iTOA = 1
else
iSFC = 1
iTOA = Model%levs
endif

! #######################################################################################
! Compute cosine of zenith angle (only when SW is called)
Expand Down Expand Up @@ -166,20 +153,20 @@ subroutine GFS_rrtmgp_sw_pre_run (Model, Interstitial, Grid, Sfcprop, Statein,
! #######################################################################################
! Call module_radiation_aerosols::setaer(),to setup aerosols property profile
! #######################################################################################
call setaer(p_lev, p_lay, Statein%prslk(1:NCOL,iSFC:iTOA), tv_lay, relhum, &
call setaer(p_lev, p_lay, Statein%prslk(1:NCOL,:), tv_lay, relhum, &
Sfcprop%slmsk, tracer, Grid%xlon, Grid%xlat, NCOL, Model%levs, Model%levs+1, &
Model%lsswr, .true., aerosolssw2, aerosolslw, aerodp)

! Store aerosol optical properties
! SW.
! For RRTMGP SW the bands are now ordered from [IR(band) -> nIR -> UV], in RRTMG the
! band ordering was [nIR -> UV -> IR(band)]
aerosolssw(1:NCOL,1:Model%levs,1,1) = aerosolssw2(1:NCOL,1:Model%levs,sw_gas_props%get_nband(),1)
aerosolssw(1:NCOL,1:Model%levs,1,2) = aerosolssw2(1:NCOL,1:Model%levs,sw_gas_props%get_nband(),2)
aerosolssw(1:NCOL,1:Model%levs,1,3) = aerosolssw2(1:NCOL,1:Model%levs,sw_gas_props%get_nband(),3)
aerosolssw(1:NCOL,1:Model%levs,2:sw_gas_props%get_nband(),1) = aerosolssw2(1:NCOL,1:Model%levs,1:sw_gas_props%get_nband()-1,1)
aerosolssw(1:NCOL,1:Model%levs,2:sw_gas_props%get_nband(),2) = aerosolssw2(1:NCOL,1:Model%levs,1:sw_gas_props%get_nband()-1,2)
aerosolssw(1:NCOL,1:Model%levs,2:sw_gas_props%get_nband(),3) = aerosolssw2(1:NCOL,1:Model%levs,1:sw_gas_props%get_nband()-1,3)
aerosolssw(1:NCOL,:,1,1) = aerosolssw2(1:NCOL,:,sw_gas_props%get_nband(),1)
aerosolssw(1:NCOL,:,1,2) = aerosolssw2(1:NCOL,:,sw_gas_props%get_nband(),2)
aerosolssw(1:NCOL,:,1,3) = aerosolssw2(1:NCOL,:,sw_gas_props%get_nband(),3)
aerosolssw(1:NCOL,:,2:sw_gas_props%get_nband(),1) = aerosolssw2(1:NCOL,:,1:sw_gas_props%get_nband()-1,1)
aerosolssw(1:NCOL,:,2:sw_gas_props%get_nband(),2) = aerosolssw2(1:NCOL,:,1:sw_gas_props%get_nband()-1,2)
aerosolssw(1:NCOL,:,2:sw_gas_props%get_nband(),3) = aerosolssw2(1:NCOL,:,1:sw_gas_props%get_nband()-1,3)

end subroutine GFS_rrtmgp_sw_pre_run

Expand Down
32 changes: 16 additions & 16 deletions physics/rrtmgp_lw_rte.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,10 +79,10 @@ subroutine rrtmgp_lw_rte_run(Model, Statein, Interstitial, ncol, lw_gas_props, p
! Local variables
type(ty_fluxes_byband) :: &
flux_allsky, flux_clrsky
real(kind_phys), dimension(ncol,model%levs+1),target :: &
fluxLW_up_allsky, fluxLW_up_clrsky, fluxLW_dn_allsky, fluxLW_dn_clrsky
real(kind_phys), dimension(ncol,model%levs+1,lw_gas_props%get_nband()),target :: &
fluxLWBB_up_allsky, fluxLWBB_dn_allsky
fluxLW_up_allsky, fluxLW_up_clrsky, fluxLW_dn_allsky, fluxLW_dn_clrsky
! real(kind_phys), dimension(ncol,model%levs+1,lw_gas_props%get_nband()),target :: &
! fluxLWBB_up_allsky, fluxLWBB_dn_allsky
logical :: &
l_ClrSky_HR, l_AllSky_HR_byband, top_at_1

Expand All @@ -92,22 +92,22 @@ subroutine rrtmgp_lw_rte_run(Model, Statein, Interstitial, ncol, lw_gas_props, p
if (.not. lslwr) return

! Vertical ordering?
top_at_1 = (Statein%prsi(1,1) .lt. Statein%prsi(1, Model%levs))
top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs))

! Are any optional outputs requested? Need to know now to compute correct fluxes.
l_ClrSky_HR = present(hlw0)
l_AllSky_HR_byband = present(hlwb)

! Initialize RRTMGP DDT containing 2D(3D) fluxes
flux_allsky%flux_up => fluxLW_up_allsky
flux_allsky%flux_dn => fluxLW_dn_allsky
flux_clrsky%flux_up => fluxLW_up_clrsky
flux_clrsky%flux_dn => fluxLW_dn_clrsky
flux_allsky%bnd_flux_up => fluxLW_up_allsky
flux_allsky%bnd_flux_dn => fluxLW_dn_allsky
flux_clrsky%bnd_flux_up => fluxLW_up_clrsky
flux_clrsky%bnd_flux_dn => fluxLW_dn_clrsky
! Only calculate fluxes by-band, only when heating-rate profiles by band are requested.
if (l_AllSky_HR_byband) then
flux_allsky%bnd_flux_up => fluxLWBB_up_allsky
flux_allsky%bnd_flux_dn => fluxLWBB_dn_allsky
endif
!if (l_AllSky_HR_byband) then
! flux_allsky%bnd_flux_up => fluxLWBB_up_allsky
! flux_allsky%bnd_flux_dn => fluxLWBB_dn_allsky
!endif

! Compute clear-sky fluxes (if requested)
! Clear-sky fluxes are gas+aerosol
Expand All @@ -120,8 +120,8 @@ subroutine rrtmgp_lw_rte_run(Model, Statein, Interstitial, ncol, lw_gas_props, p
Interstitial%sfc_emiss_byband, & ! IN - surface emissivity in each LW band
flux_clrsky))
! Store fluxes
fluxlwUP_clrsky = flux_clrsky%flux_up
fluxlwDOWN_clrsky = flux_clrsky%flux_dn
fluxlwUP_clrsky = sum(flux_clrsky%bnd_flux_up,dim=3)
fluxlwDOWN_clrsky = sum(flux_clrsky%bnd_flux_dn,dim=3)
endif

! All-sky fluxes
Expand All @@ -134,8 +134,8 @@ subroutine rrtmgp_lw_rte_run(Model, Statein, Interstitial, ncol, lw_gas_props, p
Interstitial%sfc_emiss_byband, & ! IN - surface emissivity in each LW band
flux_allsky))
! Store fluxes
fluxlwUP_allsky = flux_allsky%flux_up
fluxlwDOWN_allsky = flux_allsky%flux_dn
fluxlwUP_allsky = sum(flux_allsky%bnd_flux_up,dim=3)
fluxlwDOWN_allsky = sum(flux_allsky%bnd_flux_dn,dim=3)

end subroutine rrtmgp_lw_rte_run

Expand Down
Loading

0 comments on commit 2752142

Please sign in to comment.