Skip to content

all-sky ATMS radiance assimilation and using surface-sensitive ATMS radiances #207

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

Open
gmao-yzhu wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

all-sky ATMS radiance assimilation and using surface-sensitive ATMS radiances #207

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5e19ae6
commit changes for all-sky and surface-sensitive ATMS
gmao-yzhu Jul 14, 2025
14c91de
commit more changes for all-sky and surface-sensitive ATMS
gmao-yzhu Aug 14, 2025
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
15 changes: 12 additions & 3 deletions GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1285,6 +1285,7 @@ subroutine Run ( gc, import, export, clock, rc )
real(r_single),dimension(:,:,:), pointer :: dqlmr ! cloud liquid mixing ratio
real(r_single),dimension(:,:,:), pointer :: dqrmr ! rain mixing ratio
real(r_single),dimension(:,:,:), pointer :: dqsmr ! snow mixing ratio
real(r_single),dimension(:,:,:), pointer :: dclfr ! cloud fraction
real(r_single),dimension(:,: ), pointer :: dfrland ! land fraction
real(r_single),dimension(:,: ), pointer :: dfrlandice ! land-ice fraction
real(r_single),dimension(:,: ), pointer :: dfrlake ! lake fraction
Expand Down Expand Up @@ -1764,6 +1765,9 @@ subroutine GSI_GridCompGetPointers_()
case ('qstot')
call ESMFL_StateGetPointerToData(export, dqsmr, trim(cvar), alloc=.true., rc=STATUS)
VERIFY_(STATUS)
case ('cloud')
call ESMFL_StateGetPointerToData(export, dclfr, trim(cvar), alloc=.true., rc=STATUS)
VERIFY_(STATUS)
case ('ozone')
call ESMFL_StateGetPointerToData(export, doz, trim(cvar), alloc=.true., rc=STATUS)
VERIFY_(STATUS)
Expand Down Expand Up @@ -3010,6 +3014,8 @@ subroutine GSI_GridCompCopyInternal2Export_(lit)
CALL GSI_GridCompSwapJI_(dqrmr,GSI_MetGuess_Bundle(it)%r3(ipnt)%q)
case ('qs')
CALL GSI_GridCompSwapJI_(dqsmr,GSI_MetGuess_Bundle(it)%r3(ipnt)%q)
case ('cf')
CALL GSI_GridCompSwapJI_(dclfr,GSI_MetGuess_Bundle(it)%r3(ipnt)%q)
end select
endif
enddo
Expand Down Expand Up @@ -4062,21 +4068,24 @@ subroutine GSI_GridCompSetupSpecs (GC, opthw, rc)

! Declare export 3d-fields (extra met-guess)
! ------------------------
integer, parameter :: nex3dx=4
integer, parameter :: nex3dx=5
character(len=16), parameter :: exsname3dx(nex3dx) = (/ &
'qitot ', &
'qltot ', &
'qrtot ', &
'qstot ' /)
'qstot ', &
'cloud ' /)
character(len=40), parameter :: exlname3dx(nex3dx) = (/ &
'mass fraction of cloud ice water inc ', &
'mass fraction of cloud liquid water inc ', &
'mass fraction of rain inc ', &
'mass fraction of snow inc ' /)
'mass fraction of snow inc ', &
'mass fraction of cloud inc ' /)
character(len=16), parameter :: exunits3dx(nex3dx) = (/ &
'1 ', &
'1 ', &
'1 ', &
'1 ', &
'1 ' /)

! Declare export 2d-fields for trace gases - this needs
Expand Down
114 changes: 72 additions & 42 deletions GEOSaana_GridComp/GSI_GridComp/crtm_interface.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1014,6 +1014,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
! 2015-09-10 zhu - generalize enabling all-sky and aerosol usage in radiance assimilation,
! use n_clouds_fwd_wk,n_aerosols_fwd_wk,cld_sea_only_wk, cld_sea_only_wk,cw_cv,etc
! 2019-03-22 Wei/Martin - added VIIRS AOD obs in addition to MODIS AOD obs
! 2025-07-10 zhu/wei/jinajun - adapt all-sky atms, cloud fraction, cloud overlap into GEOS
!
! input argument list:
! obstype - type of observations for which to get profile
Expand Down Expand Up @@ -1079,6 +1080,11 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
use obsmod, only: iadate
use aeroinfo, only: nsigaerojac
use chemmod, only: lread_ext_aerosol !for separate aerosol input file
use crtm_cloudcover_define, only: cloudcover_maximum_overlap, &
cloudcover_random_overlap, &
cloudcover_maxran_overlap, &
cloudcover_average_overlap, & !default
cloudcover_overcast_overlap

implicit none

Expand Down Expand Up @@ -2001,22 +2007,38 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
! Include cloud guess profiles in mw radiance computation

if (n_clouds_fwd_wk>0) then
atmosphere(1)%cloud_fraction(k) = zero
do ii=1,n_clouds_fwd_wk
atmosphere(1)%cloud(ii)%water_content(:) = zero
end do

kgkg_kgm2=(atmosphere(1)%level_pressure(k)-atmosphere(1)%level_pressure(k-1))*r100/grav
if (cw_cv.or.ql_cv) then
if (icmask) then
c6(k) = kgkg_kgm2
auxdp(k)=abs(prsi_rtm(kk+1)-prsi_rtm(kk))*r10
auxq (k)=q(kk2)

if (icmask) then

! In CRTM, if cloud fraction of the layer < 1.0E-12, set cloud content and
! effective radius of all hydrometer types in that layer to zero
! CRTM minimum thresholds: cloud content=1.0E-6 and cloud fraction=1.E-12
! print*, 'crtm_interface: icfs=',icfs
if (icfs==0 .and. trim(obstype)=='atms') then
atmosphere(1)%cloud_fraction(k) = cf(kk2)
end if

c6(k) = kgkg_kgm2
auxdp(k)=abs(prsi_rtm(kk+1)-prsi_rtm(kk))*r10
auxq (k)=q(kk2)

do ii=1,n_clouds_fwd_wk
cloud_cont(k,ii)=cloud(kk2,ii)*c6(k)
cloud_efr (k,ii)=zero
end do

if ((cw_cv.or.ql_cv).and.(.not. lprecip_wk)) then

if (regional .and. (.not. wrf_mass_regional) .and. (.not. cold_start)) then
do ii=1,n_clouds_fwd_wk
cloud_cont(k,ii)=cloud(kk2,ii)*c6(k)
cloud_efr (k,ii)=cloudefr(kk2,ii)
end do
else
do ii=1,n_clouds_fwd_wk
cloud_cont(k,ii)=cloud(kk2,ii)*c6(k)
end do
end if

! total column q, tpw
Expand All @@ -2027,63 +2049,70 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
if(n_clouds_fwd_wk > 3) snow_guess = snow_guess + cloud_cont(k,4)
do ii=1,n_clouds_fwd_wk
if ( trim(obstype) == 'amsr2' .and. atmosphere(1)%level_pressure(k) <= 50.0_r_kind ) cycle
hwp_guess(ii) = hwp_guess(ii) + cloud_cont(k,ii)
hwp_guess(ii) = hwp_guess(ii) + cloud_cont(k,ii)
enddo
atmosphere(1)%cloud_fraction(k) = zero

do ii=1,n_clouds_fwd_wk
if (ii==1 .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) &
cloud_cont(k,1)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,1))
if (ii==2 .and. atmosphere(1)%temperature(k)<t0c) &
cloud_cont(k,2)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,2))
if(cloud_cont(k,ii) > 1.000_r_kind*1.0E-6_r_kind) then
atmosphere(1)%cloud_fraction(k) = one
end if
end do
!crtm2.3.x if (.not. regional .and. icfs==0 ) atmosphere(1)%cloud_fraction(k) = cf(kk2)
endif
else
if (icmask) then
c6(k) = kgkg_kgm2
do ii=1,n_clouds_fwd_wk
!cloud_cont(k,ii)=cloud(kk2,ii)*kgkg_kgm2
cloud_cont(k,ii)=cloud(kk2,ii)*c6(k)
if (imp_physics==11 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then
else ! .not. (cw_cv.or.ql_cv)
if (imp_physics==11 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then
do ii=1,n_clouds_fwd_wk
cloud_efr (k,ii)=cloudefr(kk2,ii)
else
cloud_efr (k,ii)=zero
endif
enddo
enddo
endif

if (cloud_cont(k,1) >= 1.0e-6_r_kind) clw_guess = clw_guess + cloud_cont(k,1)
! total column q, tpw
if(present(tpwc_guess)) tpwc_guess = tpwc_guess + q(kk2)*c6(k)
clw_guess = clw_guess + cloud_cont(k,1)
ciw_guess = ciw_guess + cloud_cont(k,2)
if(n_clouds_fwd_wk > 2) rain_guess = rain_guess + cloud_cont(k,3)
if(n_clouds_fwd_wk > 3) snow_guess = snow_guess + cloud_cont(k,4)
if(present(tcwv)) &
tcwv = tcwv + (atmosphere(1)%absorber(k,1)*0.001_r_kind)*c6(k)
do ii=1,n_clouds_fwd_wk
if (cloud_cont(k,ii) >= 1.0e-6_r_kind) hwp_guess(ii) = hwp_guess(ii) + cloud_cont(k,ii)
hwp_guess(ii) = hwp_guess(ii) + cloud_cont(k,ii)
enddo

!Add lower bound to all hydrometers
!note: may want to add lower bound value for effective radius

atmosphere(1)%cloud_fraction(k) = zero
do ii=1,n_clouds_fwd_wk
if (trim(cloud_names_fwd(ii))=='ql' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) &
if (trim(cloud_names_fwd(ii))=='ql' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
if (trim(cloud_names_fwd(ii))=='qi' .and. atmosphere(1)%temperature(k)<t0c) &
end if
if (trim(cloud_names_fwd(ii))=='qi' .and. atmosphere(1)%temperature(k)<t0c) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
if (trim(cloud_names_fwd(ii))=='qr' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) &
end if
if (trim(cloud_names_fwd(ii))=='qr' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
if (trim(cloud_names_fwd(ii))=='qs' .and. atmosphere(1)%temperature(k)<t0c) &
end if
if (trim(cloud_names_fwd(ii))=='qs' .and. atmosphere(1)%temperature(k)<t0c) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
if (trim(cloud_names_fwd(ii))=='qg' .and. atmosphere(1)%temperature(k)<t0c) &
end if
if (trim(cloud_names_fwd(ii))=='qg' .and. atmosphere(1)%temperature(k)<t0c) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
if(cloud_cont(k,ii) > 1.000_r_kind*1.0E-6_r_kind) then
atmosphere(1)%cloud_fraction(k) = one
end if
end do
!crtm2.3.x if (.not. regional .and. icfs==0 ) atmosphere(1)%cloud_fraction(k) = cf(kk2)
end if
endif
endif
endif

do ii=1,n_clouds_fwd_wk
if (icfs==0 .and. trim(obstype)=='atms') then
if (cloud_cont(k,ii) > 1.000_r_kind*1.0E-6_r_kind .and. &
atmosphere(1)%cloud_fraction(k)<1.001_r_kind*1.0E-12_r_kind) then
atmosphere(1)%cloud_fraction(k)=1.001_r_kind*1.0E-12_r_kind
endif
else
if (cloud_cont(k,ii) > 1.000_r_kind*1.0E-6_r_kind) then
atmosphere(1)%cloud_fraction(k) = one
end if
end if
end do
endif !icmask
endif !if (n_clouds_fwd_wk>0)

! Add in a drop-off to absorber amount in the stratosphere to be in more
! agreement with ECMWF profiles. The drop-off is removed when climatological CO2 fields
Expand Down Expand Up @@ -2113,6 +2142,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &

! Set clouds for CRTM
if(n_clouds_fwd_wk>0) then
if (icfs==0 .and. trim(obstype)=='atms') options(1)%Overlap_Id = cloudcover_average_overlap()
atmosphere(1)%n_clouds = n_clouds_fwd_wk
call Set_CRTM_Cloud (msig,n_actual_clouds_wk,cloud_names,icmask,n_clouds_fwd_wk,cloud_cont,cloud_efr,jcloud,auxdp, &
atmosphere(1)%temperature,atmosphere(1)%pressure,auxq,atmosphere(1)%cloud,lprecip_wk)
Expand Down
22 changes: 22 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_cloudy_radiance_info.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ radiance_mod_instr_input::
gmi sea ex_obserr3 .true. .false.
mhs global ex_obserr3 .true. .false.
amsr2 sea ex_obserr1 .true. .false.
atms sea ex_obserr1 .true. .true.
::

obs_gmi::
Expand Down Expand Up @@ -106,3 +107,24 @@ obs_ssmi::
6 0.050 0.10 0.00
7 0.050 0.10 0.00
::

obs_atms::
! Parameters for the observation error model
! cclr [kg/m2] & ccld [kg/m2]: range of cloud amounts
! over which the main increase in error take place
! ch cclr ccld
1 0.030 0.35
2 0.030 0.35
3 0.030 0.40
4 0.030 0.45
5 0.030 0.50
6 0.100 1.00
7 0.150 1.00
16 0.020 0.35
17 0.030 0.40
18 0.030 0.50
19 0.030 0.50
20 0.030 0.50
21 0.030 0.50
22 0.030 0.50
::
Loading
Loading