Adding VIIRS-M Radiance Data Assimilation in GEOS

Author: ahrem1234

GEOSaana_GridComp/GEOSgsi_Coupler/cplr_nst.F90

@@ -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

GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt

@@ -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
@@ -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

GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_satinfo.txt

@@ -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

GEOSaana_GridComp/GSI_GridComp/etc/gsi.rc.tmpl

@@ -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
 ::

GEOSaana_GridComp/GSI_GridComp/etc/gsidiags.rc

@@ -143,6 +143,8 @@ ssu_tirosn
 tmi_trmm
 abi_g16
 abi_g18
+viirs-m_npp
+viirs-m_j1
 ::
 
 ozlist::

GEOSaana_GridComp/GSI_GridComp/gsi_nstcouplermod.f90

@@ -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
 !-------------------------------------------------------------------------
@@ -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

GEOSaana_GridComp/GSI_GridComp/gsi_obOperTypeManager.F90

@@ -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