Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adjust GFS diagnostic Aerosol Optical Depth (AOD) output to the exact 550nm #16

Merged
merged 3 commits into from
Nov 14, 2022
Merged

Adjust GFS diagnostic Aerosol Optical Depth (AOD) output to the exact 550nm #16

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1cb95b6
Adjust GFS diagnostic AOD output to the exact 550nm
ChunxiZhang-NOAA Oct 26, 2022
bc2828a
Modify a few lines of code in radiation_aerosols.f to make them prope…
ChunxiZhang-NOAA Nov 4, 2022
40d2db4
Merge remote-tracking branch 'upstream/ufs/dev' into bug/aod550
ChunxiZhang-NOAA Nov 10, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 62 additions & 17 deletions physics/radiation_aerosols.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -405,7 +405,7 @@ module module_radiation_aerosols !
! ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw band
! asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw band
real (kind=kind_phys),allocatable,save,dimension(:,:) :: &
& extrhi_grt, scarhi_grt, ssarhi_grt, asyrhi_grt
& extrhi_grt, extrhi_grt_550, scarhi_grt, ssarhi_grt, asyrhi_grt
!
!> \name relative humidity dependent aerosol optical properties:
!! species : ss001, ss002, ss003, ss004, ss005, so4,
Expand All @@ -417,7 +417,7 @@ module module_radiation_aerosols !

!> extinction coefficient
real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: &
& extrhd_grt, scarhd_grt, ssarhd_grt, asyrhd_grt
& extrhd_grt, extrhd_grt_550, scarhd_grt, ssarhd_grt, asyrhd_grt

!> gocart species
integer, parameter :: num_gc = 5
Expand Down Expand Up @@ -453,6 +453,7 @@ module module_radiation_aerosols !
real (kind=kind_phys), parameter :: wvn550 = 1.0e4/0.55
!> the sw spectral band covering wvn550 (comp in aer_init)
integer, save :: nv_aod = 1
integer :: i550,id550

! --- public interfaces

Expand Down Expand Up @@ -3488,10 +3489,12 @@ subroutine gocart_aerinit &

if ( .not. allocated( extrhi_grt ) ) then
allocate ( extrhi_grt ( kcm1,nswlwbd) )
allocate ( extrhi_grt_550 ( kcm1,1) )
allocate ( scarhi_grt ( kcm1,nswlwbd) )
allocate ( ssarhi_grt ( kcm1,nswlwbd) )
allocate ( asyrhi_grt ( kcm1,nswlwbd) )
allocate ( extrhd_grt (krhlev,kcm2,nswlwbd) )
allocate ( extrhd_grt_550 (krhlev,kcm2,1) )
allocate ( scarhd_grt (krhlev,kcm2,nswlwbd) )
allocate ( ssarhd_grt (krhlev,kcm2,nswlwbd) )
allocate ( asyrhd_grt (krhlev,kcm2,nswlwbd) )
Expand Down Expand Up @@ -3834,6 +3837,9 @@ subroutine rd_gocart_luts
do i = 1, kaerbndi ! convert from m to micron
j = kaerbndi -i + 1 ! flip i-index
wavelength_du(j) = 1.e6 * lambda_du(i)
if (int(wavelength_du(j)*100) == 55) then
id550=j
endif
enddo
do k = 1, iradius
do i = 1, kaerbndi
Expand Down Expand Up @@ -3898,6 +3904,9 @@ subroutine rd_gocart_luts
do i = 1, kaerbndd ! convert from m to micron
j = kaerbndd -i + 1 ! flip i-index
wavelength(j) = 1.e6 * lambda(i)
if (int(wavelength(j)*100) == 55) then
i550=j
endif
enddo
do k = 1, iradius
ik = ibeg + k - 1
Expand Down Expand Up @@ -4016,6 +4025,9 @@ subroutine optavg_gocart
refb = sumreft * rsolbd

extrhi_grt(nc,nb) = sumk * rsolbd
if (nb==nv_aod) then
extrhi_grt_550(nc,1) = rhidext0_grt(id550,nc)
endif
scarhi_grt(nc,nb) = sums * rsolbd
asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10)
ssarhi_grt(nc,nb) = 4.0*refb &
Expand Down Expand Up @@ -4048,6 +4060,9 @@ subroutine optavg_gocart
refb = sumreft * rsolbd

extrhd_grt(nh,nc,nb) = sumk * rsolbd
if (nb==nv_aod) then
extrhd_grt_550(nh,nc,1) = rhdpext0_grt(i550,nh,nc)
endif
scarhd_grt(nh,nc,nb) = sums * rsolbd
asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
ssarhd_grt(nh,nc,nb) = 4.0*refb &
Expand Down Expand Up @@ -4253,7 +4268,8 @@ subroutine aer_property_gocart &

! --- locals:
real (kind=kind_phys), dimension(nlay,nswlwbd):: tauae,ssaae,asyae
real (kind=kind_phys), dimension(nspc) :: spcodp
real (kind=kind_phys), dimension(nlay,1):: tauae_550
real (kind=kind_phys), dimension(nlay,nspc) :: spcodp

real (kind=kind_phys),dimension(nlay,kcm) :: aerms
real (kind=kind_phys),dimension(nlay) :: dz1, rh1
Expand All @@ -4269,20 +4285,25 @@ subroutine aer_property_gocart &
do m = 1, NSWLWBD
do k = 1, NLAY
tauae(k,m) = f_zero
if (m==nv_aod) then
tauae_550(k,1) = f_zero
endif
ssaae(k,m) = f_one
asyae(k,m) = f_zero
enddo
enddo

! --- set floor value for aerms (kg/m3)
do k = 1, NLAY
do m = 1, kcm
aerms(k,m) = 1.e-15
enddo
do m = 1, kcm
aerms(k,m) = 1.e-15
enddo
enddo

do m = 1, nspc
spcodp(m) = f_zero
do k = 1, NLAY
ChunxiZhang-NOAA marked this conversation as resolved.
Show resolved Hide resolved
do m = 1, nspc
spcodp(k,m) = f_zero
enddo
enddo

do k = 1, NLAY
Expand Down Expand Up @@ -4320,11 +4341,10 @@ subroutine aer_property_gocart &

! --- update diagnostic aod arrays
do k = 1, NLAY
aerodp(i,1) = aerodp(i,1) + tauae(k,nv_aod)
enddo

do m = 1, NSPC
aerodp(i,m+1) = spcodp(m)
aerodp(i,1) = aerodp(i,1) + tauae_550(k,1)
do m = 1, NSPC
aerodp(i,m+1) = aerodp(i,m+1)+spcodp(k,m)
enddo
enddo

endif ! end if_larsw_block
Expand Down Expand Up @@ -4390,9 +4410,10 @@ subroutine aeropt

! --- locals:
real (kind=kind_phys) :: drh0, drh1, rdrh
real (kind=kind_phys) :: cm, ext01, sca01, asy01, ssa01
real (kind=kind_phys) :: ext1, asy1, ssa1, sca1
real (kind=kind_phys) :: cm, ext01, ext01_550, sca01,asy01,ssa01
real (kind=kind_phys) :: ext1, ext1_550, asy1, ssa1,sca1,tau_550
real (kind=kind_phys) :: sum_tau,sum_asy,sum_ssa,tau,asy,ssa
real (kind=kind_phys) :: sum_tau_550
integer :: ih1, ih2, nbin, ib, ntrc, ktrc

! --- linear interp coeffs for rh-dep species
Expand All @@ -4416,6 +4437,10 @@ subroutine aeropt
do ib = 1, nswlwbd

sum_tau = f_zero
if (ib == nv_aod ) then
sum_tau_550 = f_zero
ext1_550 = f_zero
endif
sum_ssa = f_zero
sum_asy = f_zero

Expand All @@ -4429,6 +4454,9 @@ subroutine aeropt
do m = 1, kcm1
cm = max(aerms(k,m),0.0) * dz1(k)
ext1 = ext1 + cm*extrhi_grt(m,ib)
if (ib == nv_aod) then
ext1_550 = ext1_550 + cm*extrhi_grt_550(m,1)
endif
sca1 = sca1 + cm*scarhi_grt(m,ib)
ssa1 = ssa1 + cm*extrhi_grt(m,ib) * ssarhi_grt(m,ib)
asy1 = asy1 + cm*scarhi_grt(m,ib) * asyrhi_grt(m,ib)
Expand All @@ -4439,7 +4467,9 @@ subroutine aeropt

! --- update aod from individual species
if ( ib==nv_aod ) then
spcodp(1) = spcodp(1) + tau
tau_550 = ext1_550
spcodp(k,1) = tau_550
sum_tau_550 = sum_tau_550 + tau_550
endif
! --- update sum_tau, sum_ssa, sum_asy
sum_tau = sum_tau + tau
Expand All @@ -4449,6 +4479,9 @@ subroutine aeropt
! --- determine tau, ssa, asy for non-dust aerosols
do ntrc = 2, nspc
ext1 = f_zero
if ( ib==nv_aod ) then
ext1_550 = f_zero
endif
asy1 = f_zero
sca1 = f_zero
ssa1 = f_zero
Expand All @@ -4459,13 +4492,20 @@ subroutine aeropt
cm = max(aerms(k,m1),0.0) * dz1(k)
ext01 = extrhd_grt(ih1,m,ib) + &
& rdrh * (extrhd_grt(ih2,m,ib)-extrhd_grt(ih1,m,ib))
if ( ib==nv_aod ) then
ext01_550 = extrhd_grt_550(ih1,m,1) + &
& rdrh * (extrhd_grt_550(ih2,m,1)-extrhd_grt_550(ih1,m,1))
endif
sca01 = scarhd_grt(ih1,m,ib) + &
& rdrh * (scarhd_grt(ih2,m,ib)-scarhd_grt(ih1,m,ib))
ssa01 = ssarhd_grt(ih1,m,ib) + &
& rdrh * (ssarhd_grt(ih2,m,ib)-ssarhd_grt(ih1,m,ib))
asy01 = asyrhd_grt(ih1,m,ib) + &
& rdrh * (asyrhd_grt(ih2,m,ib)-asyrhd_grt(ih1,m,ib))
ext1 = ext1 + cm*ext01
if ( ib==nv_aod ) then
ext1_550 = ext1_550 + cm*ext01_550
endif
sca1 = sca1 + cm*sca01
ssa1 = ssa1 + cm*ext01 * ssa01
asy1 = asy1 + cm*sca01 * asy01
Expand All @@ -4475,7 +4515,9 @@ subroutine aeropt
if (sca1 > f_zero) asy=min(f_one, asy1/sca1)
! --- update aod from individual species
if ( ib==nv_aod ) then
spcodp(ktrc) = spcodp(ktrc) + tau
tau_550 = ext1_550
spcodp(k,ktrc) = tau_550
sum_tau_550 = sum_tau_550 + tau_550
endif
! --- update sum_tau, sum_ssa, sum_asy
sum_tau = sum_tau + tau
Expand All @@ -4485,6 +4527,9 @@ subroutine aeropt

! --- determine total tau, ssa, asy for aerosol mixture
tauae(k,ib) = sum_tau
if ( ib==nv_aod ) then
tauae_550(k,1) = sum_tau_550
endif
if (sum_tau > f_zero) ssaae(k,ib) = sum_ssa / sum_tau
if (sum_ssa > f_zero) asyae(k,ib) = sum_asy / sum_ssa

Expand Down