Skip to content

Adding VIIRS-M Radiance Data Assimilation in GEOS #237

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
ahrem1234 wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Adding VIIRS-M Radiance Data Assimilation in GEOS #237

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
279 changes: 279 additions & 0 deletions GEOSaana_GridComp/GEOSgsi_Coupler/cplr_nst.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,285 @@ subroutine deter_nst_(dlat_earth,dlon_earth,obstime,zob,tref,dtw,dtc,tz_tr)
end subroutine deter_nst_
!*******************************************************************************************

subroutine deter_nst_viirs_(dlat_earth,dlon_earth,obstime,zob,tref,dtw,dtc,tz_tr)
!$$$ subprogram documentation block
! . . . .
! subprogram: deter_nst_viirs determine NSST variable at observation location over water
! prgmmr: Xu Li org: np2 date: 2011-04-08
! abstract: determines NSST variables over water surface type based on surrounding surface types
!
! program history log:
! 2026-05-22 a.lee Modify deter_nst_viirs_ for viirs (needs further investigation)
!
! input argument list:
! dlat_earth - earth latitude in radians
! dlon_earth - earth longitude in radians
! obstime - observation time relative to analysis time
! zob - obs. depth in the water
!
! output argument list:
! tref - oceanic foundation temperature
! dtw - diurnal warming at depth zob
! dtc - sublayer cooling at depth zob
! tz_tr - d(Tz)/d(Tbar); is DIFFERENT to Xu Li. Xu Li: d(Tz)/d(Tr); Tr: foundation or ref. Temp.
! - in GEOS context it is = d(Tz)/d(Ts)
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero, one
use mpimod, only: mype
use gridmod, only: nlat,nlon,regional,tll2xy,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc
use guess_grids, only: nfldsfc,hrdifsfc
use mpimod, only: mype
use gsi_nstcouplermod, only: tref_full,dt_cool_full,dt_warm_full,z_c_full,z_w_full
use satthin, only: isli_full

use GSI_GridCompMod, only: MU_SKIN=>GEOS_MU_SKIN

! use ESMF, only: ESMF_ConfigGetAttribute
! use m_die, only: die
implicit none

real(r_kind), intent(in ) :: dlat_earth,dlon_earth,obstime,zob
real(r_kind), intent(out) :: tref,dtw,dtc,tz_tr

! local variables
character(len=*), parameter:: myname_='deter_nst_viirs_'
real(r_kind):: dt_cool,z_c,z_w,dt_warm
integer(i_kind):: istyp00,istyp01,istyp10,istyp11
integer(i_kind):: itnst,itnstp
integer(i_kind):: ix,iy,ixp,iyp,j,i,k
real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11,dtnst,dtnstp,wgtmin
real(r_kind):: tref_00,tref_01,tref_10,tref_11,tr_tmp
real(r_kind):: dt_cool_00,dt_cool_01,dt_cool_10,dt_cool_11
real(r_kind):: z_c_00,z_c_01,z_c_10,z_c_11,z_c_tmp
real(r_kind):: dt_warm_00,dt_warm_01,dt_warm_10,dt_warm_11
real(r_kind):: z_w_00,z_w_01,z_w_10,z_w_11,z_w_tmp

real(r_kind):: wgtavg,dlat,dlon
logical outside


! Read in the temperature profile exponent: mu_skin used for skin SST analysis
! ----------------------------------------------------------------------------
! CALL ESMF_ConfigGetAttribute(CF, MU_SKIN, label = 'mu_skin:', default=0.2_r_kind, rc=STATUS)
! if (status/=0) then
! call die ( myname_,': failed to get mu_skin' )
! endif
! if(IamRoot) print *,trim(Iam),': Set MU_SKIN= ', MU_SKIN

! Initialize output.
tref = zero
dtw = zero
dtc = zero
tz_tr = one


if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
else
dlat=dlat_earth
dlon=dlon_earth
call grdcrd1(dlat,rlats_sfc,nlat_sfc,1)
call grdcrd1(dlon,rlons_sfc,nlon_sfc,1)
end if

iy=int(dlon); ix=int(dlat)
dy =dlon-iy; dx =dlat-ix
dx1 =one-dx; dy1 =one-dy
w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=one-w00-w10-w01

ix=min(max(1,ix),nlat_sfc); iy=min(max(0,iy),nlon_sfc)
ixp=min(nlat_sfc,ix+1); iyp=iy+1
if(iy==0) iy=nlon_sfc
if(iyp==nlon_sfc+1) iyp=1

! Get time interpolation factors for nst files
if(obstime > hrdifsfc(1) .and. obstime <= hrdifsfc(nfldsfc))then
do j=1,nfldsfc-1
if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then
itnst=j
itnstp=j+1
dtnst=(hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j))
end if
end do
else if(obstime <=hrdifsfc(1))then
itnst=1
itnstp=1
dtnst=one
else
itnst=nfldsfc
itnstp=nfldsfc
dtnst=one
end if
dtnstp=one-dtnst

! Set surface type flag.

istyp00 = isli_full(ix ,iy )
istyp10 = isli_full(ixp,iy )
istyp01 = isli_full(ix ,iyp)
istyp11 = isli_full(ixp,iyp)
!
! Use the time interpolation factors for nst files
!
tref_00 = 0.0_r_kind
tref_01 = 0.0_r_kind
tref_10 = 0.0_r_kind
tref_11 = 0.0_r_kind

if ( tref_full(ix ,iy ,itnst) /= 290.0_r_kind .and. tref_full(ix ,iy ,itnstp) /= 290.0_r_kind .and. &
tref_full(ix ,iyp ,itnst) /= 290.0_r_kind .and. tref_full(ix ,iyp ,itnstp) /= 290.0_r_kind .and. &
tref_full(ixp ,iy ,itnst) /= 290.0_r_kind .and. tref_full(ixp ,iy ,itnstp) /= 290.0_r_kind .and. &
tref_full(ixp ,iyp ,itnst) /= 290.0_r_kind .and. tref_full(ixp ,iyp ,itnstp) /= 290.0_r_kind ) then

tref_00 = tref_full (ix ,iy ,itnst)*dtnst + tref_full (ix ,iy ,itnstp)*dtnstp
tref_01 = tref_full (ix ,iyp,itnst)*dtnst + tref_full (ix ,iyp,itnstp)*dtnstp
tref_10 = tref_full (ixp,iy ,itnst)*dtnst + tref_full (ixp,iy ,itnstp)*dtnstp
tref_11 = tref_full (ixp,iyp,itnst)*dtnst + tref_full (ixp,iyp,itnstp)*dtnstp
endif

dt_cool_00 = dt_cool_full(ix ,iy ,itnst)*dtnst + dt_cool_full(ix ,iy ,itnstp)*dtnstp
dt_cool_01 = dt_cool_full(ix ,iyp,itnst)*dtnst + dt_cool_full(ix ,iyp,itnstp)*dtnstp
dt_cool_10 = dt_cool_full(ixp,iy ,itnst)*dtnst + dt_cool_full(ixp,iy ,itnstp)*dtnstp
dt_cool_11 = dt_cool_full(ixp,iyp,itnst)*dtnst + dt_cool_full(ixp,iyp,itnstp)*dtnstp

z_c_00 = z_c_full (ix ,iy ,itnst)*dtnst + z_c_full (ix ,iy ,itnstp)*dtnstp
z_c_01 = z_c_full (ix ,iyp,itnst)*dtnst + z_c_full (ix ,iyp,itnstp)*dtnstp
z_c_10 = z_c_full (ixp,iy ,itnst)*dtnst + z_c_full (ixp,iy ,itnstp)*dtnstp
z_c_11 = z_c_full (ixp,iyp,itnst)*dtnst + z_c_full (ixp,iyp,itnstp)*dtnstp


dt_warm_00 = dt_warm_full(ix ,iy ,itnst)*dtnst + dt_warm_full(ix ,iy ,itnstp)*dtnstp
dt_warm_01 = dt_warm_full(ix ,iyp,itnst)*dtnst + dt_warm_full(ix ,iyp,itnstp)*dtnstp
dt_warm_10 = dt_warm_full(ixp,iy ,itnst)*dtnst + dt_warm_full(ixp,iy ,itnstp)*dtnstp
dt_warm_11 = dt_warm_full(ixp,iyp,itnst)*dtnst + dt_warm_full(ixp,iyp,itnstp)*dtnstp

z_w_00 = z_w_full (ix ,iy ,itnst)*dtnst + z_w_full (ix ,iy ,itnstp)*dtnstp
z_w_01 = z_w_full (ix ,iyp,itnst)*dtnst + z_w_full (ix ,iyp,itnstp)*dtnstp
z_w_10 = z_w_full (ixp,iy ,itnst)*dtnst + z_w_full (ixp,iy ,itnstp)*dtnstp
z_w_11 = z_w_full (ixp,iyp,itnst)*dtnst + z_w_full (ixp,iyp,itnstp)*dtnstp

! Interpolate nst variables to obs location (water surface only)

wgtavg = zero
tr_tmp = zero
dt_cool = zero
z_c_tmp = zero
dt_warm = zero
z_w_tmp = zero

if (istyp00 == 0)then
wgtavg = wgtavg + w00
tr_tmp = tr_tmp + w00*tref_00
dt_cool = dt_cool + w00*dt_cool_00
z_c_tmp = z_c_tmp + w00*z_c_00
dt_warm = dt_warm + w00*dt_warm_00
z_w_tmp = z_w_tmp + w00*z_w_00
end if
if(istyp01 == 0)then
wgtavg = wgtavg + w01
tr_tmp = tr_tmp + w01*tref_01
dt_cool = dt_cool + w01*dt_cool_01
z_c_tmp = z_c_tmp + w01*z_c_01
dt_warm = dt_warm + w01*dt_warm_01
z_w_tmp = z_w_tmp + w01*z_w_01
end if
if(istyp10 == 0)then
wgtavg = wgtavg + w10
tr_tmp = tr_tmp + w10*tref_10
dt_cool = dt_cool + w10*dt_cool_10
z_c_tmp = z_c_tmp + w10*z_c_10
dt_warm = dt_warm + w10*dt_warm_10
z_w_tmp = z_w_tmp + w10*z_w_10
end if
if(istyp11 == 0)then
wgtavg = wgtavg + w11
tr_tmp = tr_tmp + w11*tref_11
dt_cool = dt_cool + w11*dt_cool_11
z_c_tmp = z_c_tmp + w11*z_c_11
dt_warm = dt_warm + w11*dt_warm_11
z_w_tmp = z_w_tmp + w11*z_w_11
end if

if(wgtavg > zero)then
! tr_tmp = tr_tmp/wgtavg
! tref = tr_tmp

z_w_tmp = z_w_tmp/wgtavg
z_w = z_w_tmp
dt_warm = dt_warm/wgtavg

dt_cool = dt_cool/wgtavg
z_c_tmp = z_c_tmp/wgtavg
z_c = z_c_tmp

if(tref_00 /= 0.0_r_kind .and. &
tref_01 /= 0.0_r_kind .and. &
tref_10 /= 0.0_r_kind .and. &
tref_11 /= 0.0_r_kind )then
tr_tmp = tr_tmp/wgtavg
tref = tr_tmp
endif
! Jacobian calculation: d(T(z))/d(Ts)

tz_tr = one

if ( (zob > z_w) .AND. (zob > z_c) ) then
! should not have obs that is deeper than both z_w & z_c
! once you fix sfcpt(0)- make it agree with frac water & ice in BKG,
! this branch of the if loop should never be exercised.
! For now, set diurnal fields to zero. SA. 02/06/2012
dtc = zero
dtw = zero
tz_tr = zero ! gradient should NOT impact temperature analysis.
! if(mype==0) WRITE(885,771) tref, ix, iy, ixp, iyp
! if(mype==0) WRITE(885,771) zob, z_c, z_w, dt_cool, dt_warm
else

if (zob .le. z_c) then
dtc = (one-zob/z_c) * dt_cool ! linear T(z) profile in cool-layer
dtw = dt_warm ! account for complete warm-layer temp. rise
tz_tr = one
! if(mype==0) WRITE(886,771) zob, z_c, z_w, dt_cool, dt_warm

elseif ( (zob > z_c) .and. (zob <= 0.05)) then ! z_c < zob < 5cm. That is, zob certainly corresponds to satellite measurement (IR, MW)
dtc = zero ! sensor does not feel cool-layer effects.
dtw = dt_warm * (one- ((zob-z_c)/(z_w-z_c))**MU_SKIN) ! ZENG & BELJAARS warm layer T(z) profile
tz_tr = one ! THIS IS AN Approx. for Sattelite (MW) data. Correct way is below. --this MIGHT NEED REVISIT! 10/03/2014.

else ! ((zob > z_c) .AND. (zob .le. z_w))
dtc = zero ! sensor does not feel cool-layer effects.
dtw = dt_warm * (one- ((zob-z_c)/(z_w-z_c))**MU_SKIN) ! ZENG & BELJAARS warm layer T(z) profile
tz_tr = one- ((zob-z_c)/(z_w-z_c))**MU_SKIN ! larger zob => smaller tz_tr for MU_SKIN ~0.2- 0.3. But if MU_SKIN ~1, tz_tr ~1.
! Implies deeper (in water) obs do not change Ts as much. Which makes sense from an ATMOSPHERIC point of view.
! if(mype==0) WRITE(887,771) zob, z_c, z_w, dt_cool, dt_warm
end if
! call cal_tztr_(dt_warm, dt_cool, z_c, z_w, zob, tz_tr)
end if

end if

! z_c >=0 by definition. in GEOS.

! keep Xu Li''s code for future ref.
! dtw = fac_dtl*dt_warm*(one-min(zob,z_w)/z_w)
! if ( z_c > zero ) then
! dtc = fac_tsl*dt_cool*(one-min(zob,z_c)/z_c)
! else
! dtc = zero
! endif
! call cal_tztr_(dt_warm, dt_cool, z_c, z_w, zob, tz_tr)
!--

!771 FORMAT(F10.4, 2x,4I4)
!771 FORMAT(E12.4, 2x, 4E12.4)
end subroutine deter_nst_viirs_


!subroutine cal_tztr_(dt_warm, dt_cool, z_c, z_w, z, tztr)
!
! abstract: calculate d(Tz)/d(Ts) with T-Profile info from NSST Model
Expand Down
2 changes: 2 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -467,6 +467,7 @@ set (SRCS_SOLVER
read_anowbufr.f90
read_atms.f90
read_avhrr.f90
read_viirs.f90
read_avhrr_navy.f90
read_bufrtovs.f90
read_tgas.f90
Expand Down Expand Up @@ -661,6 +662,7 @@ set (OBJS_OPENBUFR
read_atms.f90
read_avhrr.f90
read_avhrr_navy.f90
read_viirs.f90
read_bufrtovs.f90
read_cris.f90
read_fl_hdob.f90
Expand Down
10 changes: 10 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_satinfo.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4490,3 +4490,13 @@
abi_g18 14 -1 2.000 0.000 2.500 10.000 0.000 -2 -1 -1
abi_g18 15 -1 2.000 0.000 2.500 10.000 0.000 -2 -1 -1
abi_g18 16 -1 2.000 0.000 2.500 10.000 0.000 -2 -1 -1
viirs-m_npp 12 1 0.760 0.000 1.000 10.000 0.000 1 -1 -1
viirs-m_npp 13 -1 0.800 0.000 1.000 10.000 0.000 -1 -1 -1
viirs-m_npp 14 -1 0.800 0.000 1.000 10.000 0.000 -1 -1 -1
viirs-m_npp 15 1 0.580 0.000 1.000 10.000 0.000 1 -1 -1
viirs-m_npp 16 1 0.720 0.000 1.000 10.000 0.000 1 -1 -1
viirs-m_j1 12 1 0.760 0.000 1.000 10.000 0.000 1 -1 -1
viirs-m_j1 13 -1 0.800 0.000 1.000 10.000 0.000 -1 -1 -1
viirs-m_j1 14 -1 0.800 0.000 1.000 10.000 0.000 -1 -1 -1
viirs-m_j1 15 1 0.580 0.000 1.000 10.000 0.000 1 -1 -1
viirs-m_j1 16 1 0.720 0.000 1.000 10.000 0.000 1 -1 -1
2 changes: 2 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/etc/gsi.rc.tmpl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,8 @@ OBS_INPUT::
amsr2bufr amsr2 gcom-w1 amsr2_gcom-w1 0.0 1 0 gmao_amsr2_bufr
abibufr abi g16 abi_g16 0.0 1 0 ncep_goescsr_bufr
abibufr abi g18 abi_g18 0.0 1 0 ncep_goescsr_bufr
sstvcwbufr viirs-m npp viirs-m_npp 0.0 1 0 ncep_sstvcw_bufr
sstvcwbufr viirs-m j1 viirs-m_j1 0.0 1 0 ncep_sstvcw_bufr
! The following needs to be set for oneob-test:
! prepqc dummy dummy dummy 1.0 1 0 test
::
Expand Down
2 changes: 2 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/etc/gsidiags.rc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,6 +143,8 @@ ssu_tirosn
tmi_trmm
abi_g16
abi_g18
viirs-m_npp
viirs-m_j1
::

ozlist::
Expand Down
11 changes: 11 additions & 0 deletions GEOSaana_GridComp/GSI_GridComp/gsi_nstcouplermod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
! 2012-03-05 SA- _full fields: tref, dt_cool, dt_warm, z_c, z_w, ... are declared here INSTEAD of satthin
! 2015-05-01 Li- Change the nst fields to be single precision
! 2017-09-14 LI- Change the default value to be 1 for fac_dtl & fac_tsl
! 2026-05-22 a.lee Change gsi_nstcoupler_deter for viirs (needs further investigation)
!
!EOP
!-------------------------------------------------------------------------
Expand Down Expand Up @@ -110,6 +111,16 @@ subroutine deter_nst_(dlat_earth,dlon_earth,obstime,zob,tref,dtw,dtc,tz_tr)
real(r_kind), intent(out) :: tref,dtw,dtc,tz_tr
end subroutine deter_nst_
end interface

interface gsi_nstcoupler_deter_viirs
subroutine deter_nst_viirs_(dlat_earth,dlon_earth,obstime,zob,tref,dtw,dtc,tz_tr)
use kinds, only: r_kind
implicit none

real(r_kind), intent(in ) :: dlat_earth,dlon_earth,obstime,zob
real(r_kind), intent(out) :: tref,dtw,dtc,tz_tr
end subroutine deter_nst_viirs_
end interface
!-------------------

contains
Expand Down
1 change: 1 addition & 0 deletions GEOSaana_GridComp/GSI_GridComp/gsi_obOperTypeManager.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -348,6 +348,7 @@ function dtype2index_(dtype) result(index_)
!
case("avhrr_navy"); index_= iobOper_rad
case("avhrr" ); index_= iobOper_rad
case("viirs-m" ); index_= iobOper_rad

case("tcp" ,"[tcpoper]" ); index_= iobOper_tcp

Expand Down
Loading
Loading