Skip to content

Commit

Permalink
Fixed a few bugs, some housekeeping.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dustinswales
dustinswales committed Jun 20, 2019
1 parent 9e5405c commit c445658
Show file tree
Hide file tree
Showing 7 changed files with 56 additions and 152 deletions.
43 changes: 0 additions & 43 deletions physics/GFS_rrtmg_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -734,49 +734,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input
cldcov = 0.0
endif


! DJS2019: START Hack
! Compute layer cloud fraction.
clwmin = 0.0
cldcov(:,:) = 0.0
if (.not. Model%lmfshal) then
do k = 1, LMK
do i = 1, IM
clwt = 1.0e-6 * (plyr(i,k)*0.001)
if (ccnd(i,k,1) > 0.) then
onemrh= max( 1.e-10, 1.0-rhly(i,k) )
clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
tem1 = 2000.0 / tem1
value = max( min( tem1*(ccnd(i,k,1)-clwm), 50.0 ), 0.0 )
tem2 = sqrt( sqrt(rhly(i,k)) )
cldcov(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
endif
enddo
enddo
else
do k = 1, LMK
do i = 1, IM
clwt = 1.0e-6 * (plyr(i,k)*0.001)
if (ccnd(i,k,1) > 0.) then
onemrh= max( 1.e-10, 1.0-rhly(i,k) )
clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan
if (Model%lmfdeep2) then
tem1 = xrc3 / tem1
else
tem1 = 100.0 / tem1
endif
value = max( min( tem1*(ccnd(i,k,1)-clwm), 50.0 ), 0.0 )
tem2 = sqrt( sqrt(rhly(i,k)) )
cldcov(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
endif
enddo
enddo
endif
! DJS2019: END


if (Model%imp_physics == 99 .or. Model%imp_physics == 10) then ! zhao/moorthi's prognostic cloud scheme
! or unified cloud and/or with MG microphysics

Expand Down
4 changes: 2 additions & 2 deletions physics/GFS_rrtmgp_lw_post.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,11 +120,11 @@ subroutine GFS_rrtmgp_lw_post_run (Model, Grid, Radtend, &
! #######################################################################################
top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs))
if (top_at_1) then
iSFC = Model%levs
iSFC = Model%levs+1
iTOA = 1
else
iSFC = 1
iTOA = Model%levs
iTOA = Model%levs+1
endif

! #######################################################################################
Expand Down
61 changes: 9 additions & 52 deletions physics/GFS_rrtmgp_pre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,7 @@ end subroutine GFS_rrtmgp_pre_init
!! | faersw | aerosol_optical_properties_for_shortwave_bands_01-16 | aerosol optical properties for shortwave bands 01-16 | various | 4 | real | kind_phys | out | F |
!! | mtopa | model_layer_number_at_cloud_top | vertical indices for low, middle and high cloud tops | index | 2 | integer | | out | F |
!! | mbota | model_layer_number_at_cloud_base | vertical indices for low, middle and high cloud bases | index | 2 | integer | | out | F |
!! | de_lgth | cloud_decorrelation_length | cloud decorrelation length | km | 1 | real | kind_phys | out | F |
!! | cldsa | cloud_area_fraction_for_radiation | fraction of clouds for low, middle, high, total and BL | frac | 2 | real | kind_phys | out | F |
!! | aerodp | atmosphere_optical_thickness_due_to_ambient_aerosol_particles | vertical integrated optical depth for various aerosol species | none | 2 | real | kind_phys | out | F |
!! | alb1d | surface_albedo_perturbation | surface albedo perturbation | frac | 1 | real | kind_phys | out | F |
Expand All @@ -114,7 +115,7 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
ncol, lw_gas_props, sw_gas_props, & ! IN
raddt, p_lay, t_lay, p_lev, t_lev, tsfg, tsfa, alb1d, cld_frac, cld_lwp, & ! OUT
cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, faerlw, & ! OUT
faersw, cldsa, mtopa, mbota, aerodp, nday, idxday, gas_concentrations, errmsg, errflg)
faersw, cldsa, mtopa, mbota, de_lgth, aerodp, nday, idxday, gas_concentrations, errmsg, errflg)

! Inputs
type(GFS_control_type), intent(in) :: &
Expand Down Expand Up @@ -161,7 +162,7 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
errmsg ! Error message
integer, intent(out) :: &
errflg ! Error flag
real(kind_phys), dimension(ncol,Model%levr),intent(out) :: &
real(kind_phys), dimension(ncol,Model%levs),intent(out) :: &
cld_frac, & ! Total cloud fraction
cld_lwp, & ! Cloud liquid water path
cld_reliq, & ! Cloud liquid effective radius
Expand All @@ -180,15 +181,16 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
mtopa ! Vertical indices for cloud bases
real(kind_phys), dimension(ncol,5), intent(out) :: &
cldsa ! Fraction of clouds for low, middle, high, total and BL
real(kind_phys), dimension(ncol), intent(out) :: &
de_lgth !
real(kind_phys), dimension(ncol,NSPC1), intent(out) :: &
aerodp ! Vertical integrated optical depth for various aerosol species

! Local variables
integer :: i, j, k, iCol, iBand, iSFC, iTOA, iLay
integer :: i, j, iCol, iBand, iSFC, iTOA, iLay
logical :: top_at_1
real(kind_phys),dimension(NCOL,Model%levs) :: vmr_o3, vmr_h2o
real(kind_phys) :: es, qs
real(kind_phys), dimension(ncol) :: de_lgth
real(kind_phys), dimension(ncol, NF_ALBD) :: sfcalb
real(kind_phys), dimension(ncol, Model%levs) :: relhum, qs_lay, q_lay, deltaZ, tv_lay,&
deltaP, o3_lay
Expand Down Expand Up @@ -239,10 +241,10 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop,
enddo

! Guard against case when model uppermost model layer higher than rrtmgp allows.
where(p_lev(1:nCol,iTOA+1) .lt. lw_gas_props%get_press_min())
where(p_lev(1:nCol,iTOA+1) .lt. sw_gas_props%get_press_min())
! Set to RRTMGP min(pressure/temperature)
p_lev(1:nCol,iTOA+1) = spread(lw_gas_props%get_press_min(),dim=1,ncopies=ncol)
t_lev(1:nCol,iTOA+1) = spread(lw_gas_props%get_temp_min(),dim=1,ncopies=ncol)
p_lev(1:nCol,iTOA+1) = spread(sw_gas_props%get_press_min(),dim=1,ncopies=ncol)
t_lev(1:nCol,iTOA+1) = spread(sw_gas_props%get_temp_min(),dim=1,ncopies=ncol)
! Recompute layer pressure/temperature.
p_lay(1:NCOL,iTOA) = 0.5_kind_phys*(p_lev(1:NCOL,iTOA) + p_lev(1:NCOL,iTOA+1))
t_lay(1:NCOL,iTOA) = 0.5_kind_phys*(t_lev(1:NCOL,iTOA) + t_lev(1:NCOL,iTOA+1))
Expand Down Expand Up @@ -468,7 +470,6 @@ subroutine cloud_microphysics(Model, Tbd, Grid, Sfcprop, ncol, tracer, p_lay, t_
! Local variables
real(kind_phys), dimension(ncol, Model%levs, Model%ncnd) :: cld_condensate
integer :: i,k
real(kind_phys) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2
real(kind_phys), parameter :: xrc3 = 100.
real(kind_phys), dimension(ncol, Model%levs) :: delta_q, cnv_w, cnv_c, effr_l, effr_i, effr_r, effr_s, cldcov

Expand Down Expand Up @@ -571,50 +572,6 @@ subroutine cloud_microphysics(Model, Tbd, Grid, Sfcprop, ncol, tracer, p_lay, t_
cldcov = 0.0
endif

! #######################################################################################
! This is a hack to get the first-column in a file to contain a cloud.
! #######################################################################################
! DJS2019: START
! Compute layer cloud fraction.
clwmin = 0.0
cldcov(:,:) = 0.0
if (.not. Model%lmfshal) then
do k = 1, Model%levs
do i = 1, NCOL
clwt = 1.0e-6 * (p_lay(i,k)*0.1)
if (cld_condensate(i,k,1) > 0.) then
onemrh= max( 1.e-10, 1.0-relhum(i,k) )
clwm = clwmin / max( 0.01, p_lay(i,k)*0.1 )
tem1 = min(max(sqrt(sqrt(onemrh*qs_lay(i,k))),0.0001),1.0)
tem1 = 2000.0 / tem1
value = max( min( tem1*(cld_condensate(i,k,1)-clwm), 50.0 ), 0.0 )
tem2 = sqrt( sqrt(relhum(i,k)) )
cldcov(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
endif
enddo
enddo
else
do k = 1, Model%levs
do i = 1, NCOL
clwt = 1.0e-6 * (p_lay(i,k)*0.1)
if (cld_condensate(i,k,1) .gt. 0) then
onemrh= max( 1.e-10, 1.0-relhum(i,k) )
clwm = clwmin / max( 0.01, p_lay(i,k)*0.1 )
tem1 = min(max((onemrh*qs_lay(i,k))**0.49,0.0001),1.0) !jhan
if (Model%lmfdeep2) then
tem1 = xrc3 / tem1
else
tem1 = 100.0 / tem1
endif
value = max( min( tem1*(cld_condensate(i,k,1)-clwm), 50.0 ), 0.0 )
tem2 = sqrt( sqrt(relhum(i,k)) )
cldcov(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
endif
enddo
enddo
endif
! DJS2019: END

! #######################################################################################
! MICROPHYSICS
! #######################################################################################
Expand Down
38 changes: 19 additions & 19 deletions physics/GFS_rrtmgp_sw_post.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,51 +149,51 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Grid, Diag, Radtend, Coupling, &
! Initialize outputs
hswc(:,:) = 0.
topflx_sw = topfsw_type ( 0., 0., 0. )
sfcflx_sw = sfcfsw_type ( 0., 0., 0., 0. )
! sfcflx_sw = 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, &
fluxswDOWN_clrsky, &
fluxswUP_clrsky(idxday,:), &
fluxswDOWN_clrsky(idxday,:), &
p_lev(idxday,1:Model%levs+1), &
thetaTendClrSky))
hsw0(idxday,:)=thetaTendClrSky
endif
! All-sky heating-rate (mandatory)
call check_error_msg('GFS_rrtmgp_post',compute_heating_rate( &
fluxswUP_allsky, &
fluxswDOWN_allsky, &
fluxswUP_allsky(idxday,:), &
fluxswDOWN_allsky(idxday,:), &
p_lev(idxday,1:Model%levs+1), &
thetaTendAllSky))
hswc(idxday,:) = thetaTendAllSky

! Copy fluxes from RRTGMP types into model radiation types.
! Mandatory outputs
topflx_sw%upfxc = fluxswUP_allsky(:,iTOA)
topflx_sw%upfx0 = fluxswUP_clrsky(:,iTOA)
topflx_sw%dnfxc = fluxswDOWN_allsky(:,iTOA)
sfcflx_sw%upfxc = fluxswUP_allsky(:,iSFC)
sfcflx_sw%upfx0 = fluxswUP_clrsky(:,iSFC)
sfcflx_sw%dnfxc = fluxswDOWN_allsky(:,iSFC)
sfcflx_sw%dnfx0 = fluxswDOWN_clrsky(:,iSFC)
topflx_sw(idxday)%upfxc = fluxswUP_allsky(idxday,iTOA)
topflx_sw(idxday)%upfx0 = fluxswUP_clrsky(idxday,iTOA)
topflx_sw(idxday)%dnfxc = fluxswDOWN_allsky(idxday,iTOA)
sfcflx_sw(idxday)%upfxc = fluxswUP_allsky(idxday,iSFC)
sfcflx_sw(idxday)%upfx0 = fluxswUP_clrsky(idxday,iSFC)
sfcflx_sw(idxday)%dnfxc = fluxswDOWN_allsky(idxday,iSFC)
sfcflx_sw(idxday)%dnfx0 = fluxswDOWN_clrsky(idxday,iSFC)

! Optional output
if(l_fluxessw2D) then
flxprf_sw%upfxc = fluxswUP_allsky
flxprf_sw%dnfxc = fluxswDOWN_allsky
flxprf_sw%upfx0 = fluxswUP_clrsky
flxprf_sw%dnfx0 = fluxswDOWN_clrsky
flxprf_sw(idxday,:)%upfxc = fluxswUP_allsky(idxday,:)
flxprf_sw(idxday,:)%dnfxc = fluxswDOWN_allsky(idxday,:)
flxprf_sw(idxday,:)%upfx0 = fluxswUP_clrsky(idxday,:)
flxprf_sw(idxday,:)%dnfx0 = fluxswDOWN_clrsky(idxday,:)
endif
endif

Expand Down
22 changes: 2 additions & 20 deletions physics/radlw_main.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -354,9 +354,6 @@ module rrtmg_lw
! ================
subroutine rrtmg_lw_init ()
open(59,file='rrtmg_aux_dump.txt',status='unknown')
open(60,file='rrtmg_aux_tautot.txt',status='unknown')
open(61,file='rrtmg_aux_taucld.txt',status='unknown')
end subroutine rrtmg_lw_init
!> \defgroup module_radlw_main GFS radlw Main
Expand Down Expand Up @@ -1289,28 +1286,14 @@ subroutine rrtmg_lw_run &
endif
endif ! if_ivflip
write(59,*) "#"
write(60,*) "#"
do j=1,nLay
write(59,"(9F10.3)") plyr(1,j),tlyr(1,j),cld_lwp(1,j), &
& cld_iwp(1,j), cld_cf(1,j), sum(totuclfl(j-1:j))/2., &
& sum(totdclfl(j-1:j))/2., sum(totuflux(j-1:j))/2., &
& sum(totdflux(j-1:j))/2.
write(60,*) tautot(:,j)
write(61,*) taucld(:,j)
enddo
enddo lab_do_iplon
enddo lab_do_iplon
!...................................
end subroutine rrtmg_lw_run
!-----------------------------------
!> @}
subroutine rrtmg_lw_finalize ()
close(59)
close(60)
close(61)
end subroutine rrtmg_lw_finalize
Expand Down Expand Up @@ -1742,7 +1725,6 @@ subroutine cldprop &
do ib = 1, nbands
tauliq(ib) = max(f_zero, cldliq*(absliq1(index,ib) &
& + fint*(absliq1(index+1,ib)-absliq1(index,ib)) ))

enddo
endif ! end if_ilwcliq_block
endif ! end if_cldliq_block
Expand Down Expand Up @@ -1802,7 +1784,7 @@ subroutine cldprop &
endif ! end if_cldice_block
do ib = 1, nbands
taucld(ib,k) = tauice(ib) + tauliq(ib) + tauran + tausnw
taucld(ib,k) = tauice(ib) + tauliq(ib) + tauran + tausnw
enddo
endif lab_if_cld
Expand Down
18 changes: 9 additions & 9 deletions physics/rrtmg_sw_cloud_optics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2069,11 +2069,11 @@ subroutine rrtmg_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cld
asy_cld ! In-cloud asymmetry parameter (1)

! Local variables
integer :: iCol, iLay, iBand, index, ia, istr
integer :: iCol, iLay, iBand, index, ia
real(kind_phys) :: tau_rain, tau_snow, factor, fint, cld_ref_iceTemp,asyw,ssaw,za1,za2

real(kind_phys), dimension(nBandsSW) :: ssa_rain, ssa_snow, asy_rain, asy_snow, &
tau_liq, ssa_liq, asy_liq, tau_ice, ssa_ice, asy_ice, forwliq, asycoliq, &
tau_liq, ssa_liq, asy_liq, tau_ice, ssa_ice, asy_ice, asycoliq, &
forwice, extcoice, asycoice, ssacoice, fdelta, extcoliq, ssacoliq

! Initialize
Expand All @@ -2098,7 +2098,7 @@ subroutine rrtmg_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cld
asy_ice(:) = 0._kind_phys
asy_rain(:) = 0._kind_phys
asy_snow(:) = 0._kind_phys
if (cld_frac(iCol,iLay) .gt. 0._kind_phys) then
if (cld_frac(iCol,iLay) .gt. 1.e-12_kind_phys) then
! ###########################################################################
! Rain clouds
! ###########################################################################
Expand Down Expand Up @@ -2231,20 +2231,20 @@ subroutine rrtmg_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cld
! ###########################################################################
! Compute total cloud radiative properties (tau, omega, and g)
! ###########################################################################
if (cld_frac(iCol,iLay) .gt. 0._kind_phys) then
if (cld_frac(iCol,iLay) .gt. 1.e-12_kind_phys) then
do iBand = 1,nBandsSW
! Sum up radiative properties by type.
tau_cld(iCol,iLay,iBand) = tau_liq(iBand) + tau_ice(iBand) + tau_rain + tau_snow
ssa_cld(iCol,iLay,iBand) = ssa_liq(iBand) + ssa_ice(iBand) + ssa_rain(iBand) + ssa_snow(iBand)
asy_cld(iCol,iLay,iBand) = asy_liq(iBand) + asy_ice(iBand) + asy_rain(iBand) + asy_snow(iBand)
tau_cld(iCol,iLay,iBand) = max(1.e-12_kind_phys, tau_liq(iBand) + tau_ice(iBand) + tau_rain + tau_snow)
ssa_cld(iCol,iLay,iBand) = max(1.e-12_kind_phys, ssa_liq(iBand) + ssa_ice(iBand) + ssa_rain(iBand) + ssa_snow(iBand))
asy_cld(iCol,iLay,iBand) = max(1.e-12_kind_phys, asy_liq(iBand) + asy_ice(iBand) + asy_rain(iBand) + asy_snow(iBand))
! Delta-scale
asyw = asy_cld(iCol,iLay,iBand)/max(0._kind_phys, ssa_cld(iCol,iLay,iBand))
asyw = asy_cld(iCol,iLay,iBand)/max(1.e-12_kind_phys, ssa_cld(iCol,iLay,iBand))
ssaw = min(1._kind_phys-0.000001, ssa_cld(iCol,iLay,iBand)/tau_cld(iCol,iLay,iBand))
za1 = asyw * asyw
za2 = ssaw * za1
tau_cld(iCol,iLay,iBand) = (1._kind_phys - za2) * tau_cld(iCol,iLay,iBand)
ssa_cld(iCol,iLay,iBand) = (ssaw - za2) / (1._kind_phys - za2)
asy_cld(iCol,iLay,iBand) = (asyw - za2/ssaw)/(1-za2/ssaw)
asy_cld(iCol,iLay,iBand) = asyw/(1+asyw)
enddo ! Loop over SW bands
endif ! END sum cloudy properties
!
Expand Down
Loading

0 comments on commit c445658

Please sign in to comment.