diff --git a/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_nst.F90 b/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_nst.F90 index 57770f0..2857653 100644 --- a/GEOSaana_GridComp/GEOSgsi_Coupler/cplr_nst.F90 +++ b/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 diff --git a/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt b/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt index 2e4f811..fb27617 100644 --- a/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt +++ b/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 diff --git a/GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_satinfo.txt b/GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_satinfo.txt index 74a6a93..6c3984a 100644 --- a/GEOSaana_GridComp/GSI_GridComp/etc/gmao_global_satinfo.txt +++ b/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 diff --git a/GEOSaana_GridComp/GSI_GridComp/etc/gsi.rc.tmpl b/GEOSaana_GridComp/GSI_GridComp/etc/gsi.rc.tmpl index b28b145..0cca58c 100755 --- a/GEOSaana_GridComp/GSI_GridComp/etc/gsi.rc.tmpl +++ b/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 :: diff --git a/GEOSaana_GridComp/GSI_GridComp/etc/gsidiags.rc b/GEOSaana_GridComp/GSI_GridComp/etc/gsidiags.rc index 1e39875..3bc3b35 100755 --- a/GEOSaana_GridComp/GSI_GridComp/etc/gsidiags.rc +++ b/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:: diff --git a/GEOSaana_GridComp/GSI_GridComp/gsi_nstcouplermod.f90 b/GEOSaana_GridComp/GSI_GridComp/gsi_nstcouplermod.f90 index 1ed6fef..f5c582b 100644 --- a/GEOSaana_GridComp/GSI_GridComp/gsi_nstcouplermod.f90 +++ b/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 !------------------------------------------------------------------------- @@ -48,6 +49,7 @@ module GSI_NSTCouplerMod public GSI_NSTCoupler_read public GSI_NSTCoupler_skindepth public GSI_NSTCoupler_deter +public GSI_NSTCoupler_deter_viirs public GSI_NSTCoupler_final public :: nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl @@ -110,6 +112,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 diff --git a/GEOSaana_GridComp/GSI_GridComp/gsi_obOperTypeManager.F90 b/GEOSaana_GridComp/GSI_GridComp/gsi_obOperTypeManager.F90 index 4f817c1..19ad354 100644 --- a/GEOSaana_GridComp/GSI_GridComp/gsi_obOperTypeManager.F90 +++ b/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 diff --git a/GEOSaana_GridComp/GSI_GridComp/qcmod.f90 b/GEOSaana_GridComp/GSI_GridComp/qcmod.f90 index c488a63..09029e2 100755 --- a/GEOSaana_GridComp/GSI_GridComp/qcmod.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/qcmod.f90 @@ -79,6 +79,7 @@ module qcmod ! 2020-01-13 mkim - add allsky MHS to qc_mhs ! 2020-01-17 j.jin - updated qc_amsr2 for all-sky assimilation. ! 2022-03-04 j.jin/akkraoui - added qc for tmi and amsre +! 2026-05-22 a.lee - add qc_viirs ! ! subroutines included: ! sub init_qcvars @@ -103,6 +104,7 @@ module qcmod ! sub qc_amsr2 - qc amsr2 data ! sub qc_saphir - qc saphir data ! sub qc_abi - qc abi data +! sub qc_viirs - qc viirs-m data ! ! remarks: variable definitions below ! def dfact - factor for duplicate obs at same location for conv. data @@ -167,6 +169,7 @@ module qcmod public :: qc_msu public :: qc_irsnd public :: qc_avhrr + public :: qc_viirs public :: qc_amsua public :: qc_mhs public :: qc_atms @@ -572,7 +575,7 @@ subroutine setup_tzr_qc(obstype) tzchk = 0.50_r_kind elseif ( obstype == 'amsua' .or. obstype == 'ssmis' .or. obstype == 'ssmi' ) then tzchk = 0.12_r_kind - elseif ( obstype == 'avhrr' .or. obstype == 'avhrr_navy' ) then + elseif ( obstype == 'avhrr' .or. obstype == 'avhrr_navy' .or. obstype == 'viirs-m' ) then tzchk = 0.85_r_kind elseif ( obstype == 'hirs2' .or. obstype == 'hirs3' .or. obstype == 'hirs4' .or. & obstype == 'sndr' .or. obstype == 'sndrd1' .or. obstype == 'sndrd2'.or. & @@ -3039,6 +3042,304 @@ subroutine qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & return end subroutine qc_avhrr +subroutine qc_viirs(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & + zsges,cenlat,frac_sea,pangs,trop5,tzbgr,tsavg5,tbc,tb_obs,tnoise, & + wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, & + id_qc,aivals,errf,varinv,varinv_use,cld,cldp) +! id_qc,aivals,errf,varinv,varinv_use,cld,cldp,radmod) ! all-sky + +!$$$ subprogram documentation block +! . . . +! subprogram: qc_viirs QC for viirs +! +! prgmmr: li org: np23 date: 2011-04-08 +! +! abstract: set quality control criteria for viirs +! +! program history log: +! 2026-05-21 A.lee - added qc_viirs in GEOS system. +! +! input argument list: +! nchanl - number of channels per obs +! is - integer counter for number of observation types to process +! ndat - total number of observations types to process +! nsig - number of vertical levels +! ich - channel number +! sea - logical, sea flag +! land - logical, land flag +! ice - logical, ice flag +! snow - logical, snow flag +! luse - logical use flag +! zsges - elevation of guess +! cenlat - latitude of observation +! frac_sea - fraction of grid box covered with water +! pangs - solar zenith angle +! trop5 - tropopause pressure +! tzbgr - Tz over water +! tsavg5 - surface skin temperature +! tbc - simulated - observed BT with bias correction +! tb_obs - observed Brightness temperatures +! tnoise - channel noise array +! wavenumber - array of channel wavenumbers +! ptau5 - transmittances as a function of level and channel +! prsltmp - array of layer pressures in vertical (surface to toa) +! tvp - array of temperatures in vertical (surface to toa) +! temp - temperature sensitivity array +! wmix - moisture sensitivity array +! emissivity_k - surface emissivity sensitivity +! ts - skin temperature sensitivity +! id_qc - qc index - see qcmod definition +! aivals - array holding sums for various statistics as a function of obs type +! errf - criteria of gross error +! varinv - observation weight (modified obs var error inverse) +! varinv_use - observation weight used(modified obs var error inverse) +! +! output argument list: +! id_qc - qc index - see qcmod definition +! aivals - array holding sums for various statistics as a function of obs type +! errf - criteria of gross error +! varinv - observation weight (modified obs var error inverse) +! varinv_use - observation weight used(modified obs var error inverse) +! cld - cloud fraction +! cldp - cloud pressure +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ end documentation block + + use kinds, only: r_kind, i_kind + implicit none + +! Declare passed variables + + logical, intent(in ) :: sea,land,ice,snow,luse + integer(i_kind), intent(in ) :: ndat,nsig,nchanl,is + integer(i_kind),dimension(nchanl), intent(in ) :: ich + integer(i_kind),dimension(nchanl), intent(inout) :: id_qc + real(r_kind), intent(in ) :: zsges,cenlat,frac_sea,pangs,trop5 + real(r_kind), intent(in ) :: tzbgr,tsavg5 + real(r_kind), intent( out) :: cld,cldp + real(r_kind),dimension(40,ndat), intent(inout) :: aivals + real(r_kind),dimension(nchanl), intent(in ) :: tbc,emissivity_k,ts,wavenumber,tb_obs + real(r_kind),dimension(nchanl), intent(in ) :: tnoise + real(r_kind),dimension(nsig,nchanl),intent(in ) :: ptau5,temp,wmix + real(r_kind),dimension(nsig), intent(in ) :: prsltmp,tvp + real(r_kind),dimension(nchanl), intent(inout) :: errf,varinv,varinv_use + +! Declare local parameters + + real(r_kind),parameter:: oneover400=1.0_r_kind/400.0_r_kind + + + real(r_kind) :: demisf,dtempf,efact,dtbf,term,cenlatx,sfchgtfact + real(r_kind) :: sum1,sum2,sum3,cloudp,tmp,dts + real(r_kind),dimension(nchanl,nsig) :: dtb + integer(i_kind) :: i,k,kk,lcloud + integer(i_kind), dimension(nchanl) :: irday + real(r_kind) :: dtz,ts_ave,xindx,tzchks + + + irday = 1 + +! Reduce weight given to obs for shortwave ir if +! solar zenith angle tiny_r_kind + if (pangs <= 120.0_r_kind .and. frac_sea > zero) then + !if (pangs <= 89.0_r_kind .and. frac_sea > zero) then +! QC2 in statsrad + if(luse)aivals(9,is) = aivals(9,is) + one + do i=1,nchanl + if(wavenumber(i) > r2000)then + if(wavenumber(i) > r2400)then + varinv(i)=zero + varinv_use(i)=zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_2400_qc + irday(i) = 0 + else + tmp=one-(wavenumber(i)-r2000)*ptau5(1,i)& + *max(zero,cos(pangs*deg2rad))*oneover400 + varinv(i)=tmp*varinv(i) + varinv_use(i)=tmp*varinv_use(i) + if(id_qc(i) == igood_qc)id_qc(i)=ifail_2000_qc + end if + end if + end do + endif + + if(sea)then + demisf = r0_01 + dtempf = half + else if(land)then + demisf = r0_02 + dtempf = two + else if(ice)then + demisf = r0_03 + dtempf = four + else if(snow)then + demisf = r0_02 + dtempf = two + else + demisf = r0_03 + dtempf = four + end if + +! Reduce weight for obs over higher topography + sfchgtfact=one + if (zsges > r2000) then +! QC1 in statsrad + if(luse)aivals(8,is) = aivals(8,is) + one + sfchgtfact = (r2000/zsges)**4 + endif + +! Generate q.c. bounds and modified variances for height change and ptau5 + sum3=zero + do i=1,nchanl + if (tb_obs(i) > r1000 .or. tb_obs(i) <= zero) then + varinv(i)=zero + varinv_use(i)=zero + end if + varinv(i) = varinv(i)*(one-(one-sfchgtfact)*ptau5(1,i)) + varinv_use(i) = varinv_use(i)*(one-(one-sfchgtfact)*ptau5(1,i)) + +! Modify error based on transmittance at top of model + varinv(i)=varinv(i)*ptau5(nsig,i) + varinv_use(i)=varinv_use(i)*ptau5(nsig,i) + errf(i)=errf(i)*ptau5(nsig,i) + +! QC based on presence/absence of cloud + sum3=sum3+tbc(i)*tbc(i)*varinv_use(i) + end do + sum3=0.75_r_kind*sum3 + lcloud=0 + cld=zero + cldp=r10*prsltmp(1) +! +! Zeroing dtb since it used outside the loop in which is defined +! + dtb=zero + + do k=1,nsig + if(prsltmp(k) > trop5)then + sum1=zero + sum2=zero + do i=1,nchanl + if(varinv_use(i) > tiny_r_kind)then + dtb(i,k)=(tvp(k)-tsavg5)*ts(i) + do kk=1,k-1 + dtb(i,k)=dtb(i,k)+(tvp(k)-tvp(kk))*temp(kk,i) + end do + sum1=sum1+tbc(i)*dtb(i,k)*varinv_use(i) + sum2=sum2+dtb(i,k)*dtb(i,k)*varinv_use(i) + end if + end do + if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) + cloudp=min(max(sum1/sum2,zero),one) + sum1=zero + do i=1,nchanl + if(varinv_use(i) > tiny_r_kind)then + tmp=tbc(i)-cloudp*dtb(i,k) + sum1=sum1+tmp*tmp*varinv_use(i) + end if + end do + if(sum1 < sum3)then + sum3=sum1 + lcloud=k + cld=cloudp + cldp=r10*prsltmp(k) + end if + end if + + end do + + do i=1,nchanl + cld_qc: do k=1,lcloud + if(abs(cld*dtb(i,k)) > tnoise(i))then +! QC4 in statsrad + if(luse)aivals(11,is) = aivals(11,is) + one + varinv(i) = zero + varinv_use(i) = zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc + exit cld_qc + end if + end do cld_qc + end do + +! If no clouds check surface temperature/emissivity + + sum1=zero + sum2=zero + do i=1,nchanl + sum1=sum1+tbc(i)*ts(i)*varinv_use(i) + sum2=sum2+ts(i)*ts(i)*varinv_use(i) + end do + if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) + dts=abs(sum1/sum2) + if(abs(dts) > one)then + if(.not. sea)then + dts=min(dtempf,dts) + else + dts=min(three,dts) + end if + do i=1,nchanl + if(abs(dts*ts(i)) > tnoise(i))then +! QC3 in statsrad + if(luse .and. varinv(i) > zero) & + aivals(10,is) = aivals(10,is) + one + varinv(i) = zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc + end if + end do + end if + +! +! Apply Tz retrieval +! + if(tzr_qc > 0)then + dtz = rmiss_single + if ( sea ) then + call tz_retrieval(nchanl,nsig,ich,irday,temp,wmix,tnoise,varinv,ts,tbc,tzbgr,1,0,dtz,ts_ave) + endif +! +! Apply QC with Tz retrieval +! + if (dtz /= rmiss_single ) then + do i = 1, nchanl + if ( varinv(i) > tiny_r_kind .and. iuse_rad(ich(i)) >= 1 .and. ts(i) > tschk) then + xindx = ((ts(i)-ts_ave)/(one-ts_ave))**3 + tzchks = tzchk*(half)**xindx + + if ( abs(dtz) > tzchks ) then + varinv(i) = zero + if ( id_qc(i) == igood_qc ) id_qc(i) = ifail_tzr_qc + if(luse)aivals(13,is) = aivals(13,is) + one + endif + endif + enddo + endif + end if + + cenlatx=abs(cenlat)/30.0 + if (cenlatx < one) then + if(luse)aivals(6,is) = aivals(6,is) + one + efact = half*half*(cenlatx+one) + do i=1,nchanl + if(varinv(i) > tiny_r_kind)errf(i)=efact*errf(i) + end do + endif + + ! Generate q.c. bounds and modified variances. + do i=1,nchanl + if(varinv(i) > tiny_r_kind)then + dtbf = demisf*abs(emissivity_k(i))+dtempf*abs(ts(i)) + term = dtbf*dtbf + if(term > tiny_r_kind)varinv(i)=varinv(i)/(one+varinv(i)*term) + end if + end do + + return +end subroutine qc_viirs + subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & zsges,cenlat,tb_obsbc1,si_mean,cosza,clw,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs,cldeff_fg,factch6, & diff --git a/GEOSaana_GridComp/GSI_GridComp/radiance_mod.f90 b/GEOSaana_GridComp/GSI_GridComp/radiance_mod.f90 index 8b4c408..feec643 100644 --- a/GEOSaana_GridComp/GSI_GridComp/radiance_mod.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/radiance_mod.f90 @@ -12,6 +12,7 @@ module radiance_mod ! 2015-07-20 Yanqiu Zhu ! 2016-10-27 Yanqiu - add ATMS ! 2020-01-13 mkim - add allsky MHS +! 2026-05-22 a.lee - add viirs-m radiance assimilation ! ! subroutines included: ! sub radiance_mode_init - guess init @@ -404,6 +405,7 @@ subroutine radiance_obstype_init if (index(dtype(i),'sndr') /= 0) rtype(i)='sndr' if (index(dtype(i),'hirs') /= 0) rtype(i)='hirs' if (index(dtype(i),'avhrr') /= 0) rtype(i)='avhrr' + if (index(dtype(i),'viirs-m') /= 0) rtype(i)='viirs' if (index(dtype(i),'modis') /= 0) rtype(i)='modis' if (index(dtype(i),'seviri') /= 0) rtype(i)='seviri' @@ -414,7 +416,8 @@ subroutine radiance_obstype_init rtype(i) == 'avhrr' .or. rtype(i) == 'amsre' .or. rtype(i) == 'ssmis' .or. & rtype(i) == 'ssmi' .or. rtype(i) == 'atms' .or. rtype(i) == 'cris' .or. & rtype(i) == 'amsr2' .or. rtype(i) == 'gmi' .or. rtype(i) == 'saphir' .or. & - rtype(i) == 'cris-fsr' .or. rtype(i) == 'tmi' .or. rtype(i) == 'abi') then + rtype(i) == 'cris-fsr' .or. rtype(i) == 'tmi' .or. rtype(i) == 'abi' .or. & + rtype(i) == 'viirs' ) then drtype(i)='rads' end if end do diff --git a/GEOSaana_GridComp/GSI_GridComp/radinfo.f90 b/GEOSaana_GridComp/GSI_GridComp/radinfo.f90 index e3a27df..aa91325 100644 --- a/GEOSaana_GridComp/GSI_GridComp/radinfo.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/radinfo.f90 @@ -55,6 +55,7 @@ module radinfo ! 2016-11-29 shlyaeva - make nvarjac public ! 2018-07-24 W. Gu - the routines to handle correlated R-covariance moved out ! 2019-06-19 Hu - add option reset_bad_radbc for reset radiance bias correction coefficient if it is bad. +! 2026-05-22 a.lee - add viirs-m radiance assimilation ! ! subroutines included: ! sub init_rad - set satellite related variables to defaults @@ -1607,6 +1608,7 @@ subroutine init_predx ! 2015-10-22 jung - changed from using satinfo information in the radstat file to ! using information from the satinfo file. ! 2016-07-14 jung - mods to make SEVIRI channel numbers consistent with other instruments. +! 2026-05-22 a.lee - use of ONLY water surface for BC pred coeff. update. ! ! attributes: ! language: f90 @@ -1640,7 +1642,7 @@ subroutine init_predx logical mean_only logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,tmi,gmi,amsr2,saphir logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img - logical avhrr,avhrr_navy,goessndr,goes_img,ahi,seviri,abi + logical avhrr,avhrr_navy,viirs,goessndr,goes_img,ahi,seviri,abi character(len=20):: obstype,platid character(len=20):: satsens,satsens_id @@ -1823,6 +1825,7 @@ subroutine init_predx abi = obstype == 'abi' avhrr = obstype == 'avhrr' avhrr_navy = obstype == 'avhrr_navy' + viirs = obstype == 'viirs-m' ssmi = obstype == 'ssmi' amsre_low = obstype == 'amsre_low' amsre_mid = obstype == 'amsre_mid' @@ -1906,7 +1909,8 @@ subroutine init_predx ! angle dependent bias. if( ( abs(data_chan(j)%omgnbc) > 200. .or. & data_chan(j)%tbobs < 50. .or. & - data_chan(j)%tbobs > 500. ) ) cycle loopc + data_chan(j)%tbobs > 500. .or. & + data_fix%water_frac < 0.99 ) ) cycle loopc ! if errinv= (1 /(obs error)) is small (small = less than 1.e-6) ! the observation did not pass quality control. In this diff --git a/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 b/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 index 1b0c53a..21e415b 100644 --- a/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/read_obs.F90 @@ -684,6 +684,7 @@ subroutine read_obs(ndata,mype) ! 2019-01-15 Li - add to handle mbuoyb ! 2019-03-27 h. liu - add abi ! 2021-04-16 j.jin - read tmi and amsre bufr (made at gmao) data. +! 2026-05-22 a.lee - add viirs-m radiance assimilation ! ! ! input argument list: @@ -920,7 +921,8 @@ subroutine read_obs(ndata,mype) obstype == 'ssu' .or. obstype == 'atms' .or. & obstype == 'cris' .or. obstype == 'cris-fsr' .or. & obstype == 'amsr2' .or. obstype == 'tmi' .or. & - obstype == 'gmi' .or. obstype == 'saphir' ) then + obstype == 'gmi' .or. obstype == 'saphir' .or. & + obstype == 'viirs-m' ) then ditype(i) = 'rad' else if (is_extOzone(dfile(i),obstype,dplat(i))) then ditype(i) = 'ozone' @@ -1027,6 +1029,8 @@ subroutine read_obs(ndata,mype) parallel_read(i)= .true. else if(avhrr)then parallel_read(i)= .true. + else if(obstype == 'viirs-m')then + parallel_read(i)= .true. else if(amsre)then ! parallel_read(i)= .true. ! turn parallel read off else if(obstype == 'goes_img' )then @@ -1827,6 +1831,15 @@ subroutine read_obs(ndata,mype) mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) string='READ_AVHRR' + + ! Process SST VIIRS RADIANCE data + else if(obstype == 'viirs-m') then + call read_sst_viirs(mype,val_dat,ithin,rmesh,dplat(i),gstime,& + infile,lunout,obstype,nread,npuse,nouse,twind,sis, & + mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & + nobs_sub1(1,i),read_rec(i),dval_use) + string='READ_VIIRS' + end if rad_obstype_select ! Process ozone data diff --git a/GEOSaana_GridComp/GSI_GridComp/read_viirs.f90 b/GEOSaana_GridComp/GSI_GridComp/read_viirs.f90 new file mode 100644 index 0000000..8039e8d --- /dev/null +++ b/GEOSaana_GridComp/GSI_GridComp/read_viirs.f90 @@ -0,0 +1,521 @@ +subroutine read_sst_viirs(mype,val_viirs,ithin,rmesh,jsatid,& + gstime,infile,lunout,obstype,nread,ndata,nodata,twind,sis, & + mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs, & + nrec_start,dval_use) +!$$$ subprogram documentation block +! . . . . +! subprogram: read_sst_viirs, read M-Band (12, 15 and 16) VIIRS radiance data from the SST VIIRS data file +! +! prgmmr: li, xu org: np23 date: 2019-08-29 +! +! abstract: This routine reads BUFR format VIIRS 1b radiance (brightness +! temperature) from the SST VIIRS files +! +! When running the gsi in regional mode, the code only +! retains those observations that fall within the regional +! domain +! +! program history log: +! 2026-05-22 a.lee - added read_viirs in GEOS system. +! +! input argument list: +! mype - mpi task id +! val_viirs- weighting factor applied to super obs +! ithin - flag to thin data +! rmesh - thinning mesh size (km) +! jsatid - satellite to read +! gstime - analysis time in minutes from reference date +! infile - unit from which to read BUFR data +! lunout - unit to which to write data for further processing +! obstype - observation type to process +! twind - input group time window (hours) +! sis - satellite/instrument/sensor indicator +! +! mype_root - "root" task for sub-communicator +! mype_sub - mpi task id within sub-communicator +! npe_sub - number of data read tasks +! mpi_comm_sub - sub-communicator for data read +! nrec_start - first subset with useful information +! dval_use - true. if any dval weighting is used for satellite data +! +! output argument list: +! nread - number of BUFR VIIRS observations read +! ndata - number of BUFR VIIRS profiles retained for further processing +! nodata - number of BUFR VIIRS observations retained for further processing +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,r_single,r_double,i_kind + use satthin, only: super_val,itxmax,makegrids,map2tgrid,destroygrids, & + checkob, finalcheck,score_crit + use satthin, only: radthin_time_info,tdiff2crit + use obsmod, only: time_window_max, ianldate + use gridmod, only: diagnostic_reg,regional,nlat,nlon,tll2xy,txy2ll,rlats,rlons + use constants, only: deg2rad, zero, one, two, half, rad2deg, r60inv + use radinfo, only: iuse_rad,jpch_rad,nusis + use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen + use deter_sfc_mod, only: deter_sfc + use obsmod, only: bmiss + use gsi_nstcouplermod, only: nst_gsi,nstinfo + use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter_viirs + use sfcobsqc, only: get_sunangle + use mpimod, only: npe + implicit none + + +! Declare passed variables + character(len=*), intent(in ) :: infile,obstype,jsatid + character(len=20), intent(in ) :: sis + integer(i_kind), intent(in ) :: mype,lunout,ithin,nrec_start + integer(i_kind), intent(inout) :: nread + integer(i_kind), dimension(npe), intent(inout) :: nobs + integer(i_kind), intent(inout) :: ndata,nodata + real(r_kind), intent(in ) :: rmesh,gstime,twind + real(r_kind), intent(inout) :: val_viirs + integer(i_kind), intent(in ) :: mype_root + integer(i_kind), intent(in ) :: mype_sub + integer(i_kind), intent(in ) :: npe_sub + integer(i_kind), intent(in ) :: mpi_comm_sub + logical, intent(in ) :: dval_use + +! Declare local parameters +! character(6),parameter:: file_sst='SST_AN' +! integer(i_kind),parameter:: mlat_sst = 3000 +! integer(i_kind),parameter:: mlon_sst = 5000 + real(r_kind),parameter:: r6=6.0_r_kind + real(r_kind),parameter:: scan_start=-56.28_r_kind, scan_inc=1.26472_r_kind + real(r_double),parameter:: r360=360.0_r_double + real(r_kind),parameter:: tbmin=50.0_r_kind + real(r_kind),parameter:: tbmax=550.0_r_kind + + real(r_kind),parameter :: r_earth=6370.0_r_kind, sat_hgt=829.0_r_kind + real(r_kind),parameter :: nviirs=6304.0_r_kind,nfov=90.0_r_kind + character(len=80),parameter :: & + headr='YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAID SAZA' + +! Declare local variables + logical outside,iuse,assim + character(len=8) :: subset + + integer(i_kind) klon1,klatp1,klonp1,klat1 + integer(i_kind) nchanl,iret,ich_m15 + integer(i_kind) idate,maxinfo + integer(i_kind) ilat,ilon + integer(i_kind),dimension(5):: idate5 + integer(i_kind) nmind,isflg,idomsfc + integer(i_kind) itx,k,i,bufsat,n + integer(i_kind) nreal,nele,itt +! integer(i_kind) nlat_sst,nlon_sst + integer(i_kind) ksatid + + real(r_kind) dlon,dlat,rsc + real(r_kind) dlon_earth,dlat_earth,sfcr + real(r_kind) dlon_earth_deg,dlat_earth_deg + real(r_kind) w00,w01,w10,w11,dx1,dy1 + real(r_kind) pred,crit1,tdiff,sstime,dx,dy,dist1 +! real(r_kind) dlat_sst,dlon_sst,sst_hires + real(r_kind) t4dv, t4dv_sun + + real(r_kind),dimension(0:4):: rlndsea + real(r_kind),dimension(0:3):: sfcpct + real(r_kind),dimension(0:3):: ts +! real(r_kind),dimension(mlat_sst):: rlats_sst +! real(r_kind),dimension(mlon_sst):: rlons_sst +! real(r_kind),allocatable,dimension(:,:):: sst_an + real(r_kind),allocatable,dimension(:,:):: data_all + real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10 + real(r_kind) :: zob,tref,dtw,dtc,tz_tr + real(r_kind) :: scan_ang,sat_zen,dfov,r01 + real(r_single) :: sol_zen + integer(i_kind) :: scan_pos + + real(r_double), dimension(10) :: hdr + real(r_double), dimension(2,3) :: bufrf + integer(i_kind) lnbufr,ireadsb,ireadmg,iskip,irec,next + integer(i_kind),allocatable,dimension(:)::nrec + + real(r_kind) cdist,disterr,disterrmax,dlon00,dlat00 + integer(i_kind) ntest + + real(r_kind) :: ptime,timeinflat,crit0 + integer(i_kind) :: ithin_time,n_tbin,it_mesh + + character(len=10) :: str + integer(i_kind) :: idate_sun + +!************************************************************************** + +! Start routine here. Set constants. Initialize variables + maxinfo = 31 + lnbufr = 10 + disterrmax=zero + ntest=0 + ndata = 0 + nodata = 0 + nread = 0 + nchanl = 5 + ich_m15 = 2 + r01 = 0.01_r_kind + + dfov = (nviirs - one)/nfov + + ilon=3 + ilat=4 + + if(nst_gsi>0) then + call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype + endif + + rlndsea(0) = zero + rlndsea(1) = 30._r_kind + rlndsea(2) = 20._r_kind + rlndsea(3) = 30._r_kind + rlndsea(4) = 30._r_kind + + + if (jsatid == 'npp') then + bufsat = 224 + elseif (jsatid == 'n20' .or. jsatid == 'j1') then + bufsat = 225 + elseif (jsatid == 'n21' .or. jsatid == 'j2') then + bufsat = 226 + else + write(*,*) 'READ_SST_VIIRS: Unrecognized value for jsatid '//jsatid//':RETURNING' + return + end if + + +! If all channels of a given sensor are set to monitor or not +! assimilate mode (iuse_rad<1), reset relative weight to zero. +! We do not want such observations affecting the relative +! weighting between observations within a given thinning group. + + assim=.false. + search: do i=1,jpch_rad + if ((nusis(i)==sis) .and. (iuse_rad(i)>0)) then + assim=.true. + exit search + endif + end do search + if (.not.assim) val_viirs=zero + + + call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time) + if( ptime > 0.0_r_kind) then + n_tbin=nint(2*time_window_max/ptime) + else + n_tbin=1 + endif +! Make thinning grids + call makegrids(rmesh,ithin) + +! Read hi-res sst analysis +! allocate(sst_an(mlat_sst, mlon_sst)) +! call rdgrbsst(file_sst,mlat_sst,mlon_sst,& +! sst_an,rlats_sst,rlons_sst,nlat_sst,nlon_sst) + +! Allocate arrays to hold all data for given satellite + if(dval_use) maxinfo = maxinfo + 2 + nreal = maxinfo + nstinfo + nele = nreal + nchanl + allocate(data_all(nele,itxmax),nrec(itxmax)) + + open(lnbufr,file=trim(infile),form='unformatted') ! open bufr data file + +! Associate the tables file with the message file, and identify the +! latter to BUFRLIB software + call openbf (lnbufr,'IN',lnbufr) + + next=0 + nrec=999999 + irec=0 +! Read BUFR VIIRS 1b data + read_msg: do while (ireadmg(lnbufr,subset,idate) >= 0) + irec=irec+1 + if(irec < nrec_start) cycle read_msg + next=next+1 + if(next == npe_sub)next=0 + if(next /= mype_sub)cycle + read_loop: do while (ireadsb(lnbufr) == 0) + call ufbint(lnbufr,hdr,10,1,iret,headr) + call ufbrep(lnbufr,bufrf,2,3,iret,'CHNM TMBR') + if(iret <= 0) cycle read_loop + ksatid = nint(hdr(9)) ! Extract satellite id from bufr file + if(ksatid /= bufsat) cycle read_loop ! If this sat is not the one we want, read next record + + iskip = 0 + do k=1,nchanl-2 + if(bufrf(2,k) < zero .or. bufrf(2,k) > tbmax) then + iskip=iskip+1 + end if + end do + if(iskip >= nchanl)cycle read_loop + +! Extract date information. If time outside window, skip this obs + idate5(1) = nint(hdr(1)) !year + idate5(2) = nint(hdr(2)) !month + idate5(3) = nint(hdr(3)) !day + idate5(4) = nint(hdr(4)) !hour + idate5(5) = nint(hdr(5)) !minute + rsc = hdr(6) !second in real + call w3fs21(idate5,nmind) + write(str,'(I4.4,I2.2,I2.2,I2.2)') idate5(1), idate5(2), idate5(3), idate5(4) + read(str,*) idate_sun + t4dv=(real((nmind-iwinbgn),r_kind) + rsc*r60inv)*r60inv + t4dv_sun=(real(idate5(5),r_kind) + rsc*r60inv)*r60inv + sstime=real(nmind,r_kind) + rsc*r60inv + tdiff=(sstime-gstime)*r60inv + + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) cycle read_loop + else + if (abs(tdiff) > twind) cycle read_loop + endif + +! Convert obs location to radians + if (abs(hdr(7))>90.0_r_double .or. abs(hdr(8))>r360) cycle read_loop + if (hdr(8)==r360) hdr(8)=hdr(8)-r360 + if (hdr(8)< zero) hdr(8)=hdr(8)+r360 + + dlon_earth_deg = hdr(8) + dlat_earth_deg = hdr(7) + + dlon_earth = hdr(8)*deg2rad !convert degrees to radians + dlat_earth = hdr(7)*deg2rad + sat_zen = hdr(10) + +! Regional case + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(diagnostic_reg) then + call txy2ll(dlon,dlat,dlon00,dlat00) + ntest=ntest+1 + cdist=sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* & + (sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00)) + cdist=max(-one,min(cdist,one)) + disterr=acos(cdist)*rad2deg + disterrmax=max(disterrmax,disterr) + end if + if(outside) cycle read_loop + +! Global case: Get relative lat/lon in analysis grids + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + + nread = nread + 1 + + call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc,sfcpct, & + ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) + if(abs(sfcpct(0) - 1.0) > zero) cycle read_loop + + if(nst_gsi>0) then + tref = ts(0) + dtw = zero + dtc = zero + tz_tr = one + if(sfcpct(0)>zero) then + !call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) + call gsi_nstcoupler_deter_viirs(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) + if(tref < 1e-9) cycle read_loop + endif + endif + + crit0 = 0.01_r_kind + timeinflat=6.0_r_kind + call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) + call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh) + if(.not. iuse)cycle read_loop + +! Interpolate hi-res sst analysis to observation location +! dlat_sst = dlat_earth +! dlon_sst = dlon_earth +! call grdcrd1(dlat_sst,rlats_sst,nlat_sst,1) +! call grdcrd1(dlon_sst,rlons_sst,nlon_sst,1) +! +! klon1=int(dlon_sst); klat1=int(dlat_sst) +! dx =dlon_sst-klon1; dy =dlat_sst-klat1 +! dx1 =one-dx; dy1 =one-dy +! w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy +! +! klat1=min(max(1,klat1),nlat_sst); klon1=min(max(0,klon1),nlon_sst) +! if(klon1==0) klon1=nlon_sst +! klatp1=min(nlat_sst,klat1+1); klonp1=klon1+1 +! if(klonp1==nlon_sst+1) klonp1=1 +! +! sst_hires=w00*sst_an(klat1,klon1 ) + w10*sst_an(klatp1,klon1 ) + & +! w01*sst_an(klat1,klonp1) + w11*sst_an(klatp1,klonp1) +! + !####################################################### Only process the obs. at the location where hsst(sst_an, + !in this case?) is in the current processed range +! if ( sst_hires < zero ) then +! print*,' sst_hires,klat1,klon1 : ',sst_hires,klat1,klon1 +! endif + +! "Score" observation. We use this information to id "best" obs. + +! Locate the observation on the analysis grid. Get sst and land/sea/ice +! mask. + +! isflg - surface flag +! 0 sea +! 1 land +! 2 sea ice +! 3 snow +! 4 mixed +! call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc,sfcpct, & +! ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) +! if(isflg /= zero) print*, 'really? in viirs?' +! if(isflg /= zero) cycle read_loop + !if(isflg /= zero .or. abs(sfcpct(0)) - 1.0) .lt. 0.0001) cycle read_loop + + call checkob(dist1,crit1,itx,iuse) + if(.not. iuse)cycle read_loop + + +! Set common predictor parameters + +! Compute "score" for observation. All scores>=0.0. Lowest score is "best" + + pred = (600.0_r_kind - bufrf(2,ich_m15)) * r01 + + crit1=crit1+rlndsea(isflg) + crit1 = crit1+pred + call finalcheck(dist1,crit1,itx,iuse) + + if(.not. iuse)cycle read_loop +! +! calculate scan angle +! + scan_ang = asin(r_earth/(r_earth+sat_hgt)*sin(sat_zen*deg2rad))/deg2rad +! +! get scan position (1 to 90) +! + scan_pos = 1+(scan_ang-scan_start)/scan_inc + if ( scan_pos > nfov ) scan_pos = nfov +! +! calculate solar zenith angle + +! + call get_sunangle(idate_sun,t4dv_sun,dlon_earth,dlat_earth,sol_zen) + !call get_sunangle(ianldate,t4dv,dlon_earth,dlat_earth,sol_zen) + +! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr +! +! if(nst_gsi>0) then +! tref = ts(0) +! dtw = zero +! dtc = zero +! tz_tr = one +! if(sfcpct(0)>zero) then +! call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) +! !call gsi_nstcoupler_deter_viirs(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) +! !if(tref == zero) cycle read_loop +! endif +! endif + +! Transfer information to work array + data_all(1, itx) = hdr(9) ! satellite id + data_all(2, itx) = t4dv ! time (hours) + data_all(3, itx) = dlon ! grid relative longitude + data_all(4, itx) = dlat ! grid relative latitude + data_all(5, itx) = sat_zen*deg2rad ! satellite zenith angle (radians) + data_all(6, itx) = bmiss ! satellite azimuth angle + data_all(7, itx) = scan_ang*deg2rad ! scan angle + data_all(8, itx) = real(scan_pos) ! scan position + data_all(9, itx) = 90.0_r_kind-sol_zen ! solar zenith angle (degree) + data_all(10,itx) = bmiss ! solar azimuth angle (radians) + data_all(11,itx) = sfcpct(0) ! sea percentage of + data_all(12,itx) = sfcpct(1) ! land percentage + data_all(13,itx) = sfcpct(2) ! sea ice percentage + data_all(14,itx) = sfcpct(3) ! snow percentage + data_all(15,itx) = ts(0) ! ocean skin temperature (from surface file: sst_full) + data_all(16,itx) = ts(1) ! land skin temperature (from surface file: sst_full) + data_all(17,itx) = ts(2) ! ice skin temperature (from surface file: sst_full) + data_all(18,itx) = ts(3) ! snow skin temperature (from surface file: sst_full) + data_all(19,itx) = tsavg ! average skin temperature + data_all(20,itx) = vty ! vegetation type + data_all(21,itx) = vfr ! vegetation fraction + data_all(22,itx) = sty ! soil type + data_all(23,itx) = stp ! soil temperature + data_all(24,itx) = sm ! soil moisture + data_all(25,itx) = sn ! snow depth + data_all(26,itx) = zz ! surface height + data_all(27,itx) = idomsfc + 0.001_r_kind ! dominate surface type + data_all(28,itx) = sfcr ! surface roughness + data_all(29,itx) = ff10 ! ten meter wind factor + data_all(30,itx) = dlon_earth_deg ! earth relative longitude (degrees) + data_all(31,itx) = dlat_earth_deg ! earth relative latitude (degrees) + if(dval_use)then + data_all(32,itx) = val_viirs ! weighting factor applied to super obs + data_all(33,itx) = itt ! + end if + + if(nst_gsi>0) then + data_all(maxinfo+1,itx) = tref ! foundation temperature + data_all(maxinfo+2,itx) = dtw ! dt_warm at zob + data_all(maxinfo+3,itx) = dtc ! dt_cool at zob + data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr) + endif + + data_all(1+nreal,itx) = bufrf(2,1) + data_all(2+nreal,itx) = bmiss + data_all(3+nreal,itx) = bmiss + data_all(4+nreal,itx) = bufrf(2,2) + data_all(5+nreal,itx) = bufrf(2,3) + + nrec(itx)=irec + +! End of satellite read block + + enddo read_loop + enddo read_msg + call closbf(lnbufr) + + call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& + nele,itxmax,nread,ndata,data_all,score_crit,nrec) + +! Allow single task to check for bad obs, update superobs sum, +! and write out data to scratch file for further processing. + if (mype_sub==mype_root.and.ndata>0) then + do n=1,ndata + do k=1,nchanl + if(data_all(k+nreal,n) > tbmin .and. & + data_all(k+nreal,n) < tbmax) nodata=nodata+1 + end do + end do + if(dval_use .and. assim)then + do n=1,ndata + itt=nint(data_all(33,n)) + super_val(itt)=super_val(itt)+val_viirs + end do + end if + +! Write retained data to local file + call count_obs(ndata,nele,ilat,ilon,data_all,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon + write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata) + endif +! write(6,*) 'READ_VIIRS: mype, total number of obs info, nread,ndata : ',mype, nread,ndata + +! Deallocate local arrays + deallocate(data_all,nrec) +! deallocate(sst_an) + +! Deallocate arrays + call destroygrids + + if(diagnostic_reg.and.ntest>0 .and. mype_sub==mype_root) & + write(6,*)'READ_VIIRS-m: ',& + 'mype,ntest,disterrmax=',mype,ntest,disterrmax + +! End of routine + return +end subroutine read_sst_viirs diff --git a/GEOSaana_GridComp/GSI_GridComp/setuprad.f90 b/GEOSaana_GridComp/GSI_GridComp/setuprad.f90 index 29c22d3..1865338 100755 --- a/GEOSaana_GridComp/GSI_GridComp/setuprad.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/setuprad.f90 @@ -225,7 +225,8 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ! 2020-08-26 mkim - adjusted MHS QC for all-sky ! 2020-09-27 j.jin - assimilate SSMI, TMI and AMSRE (gmao format) in all-sky conditions. ! 2023-06-30 j.jin - add a qc_flag "ierrret_ges" to flag failed cloud retrieval from simulated Tb. -! 2024-03-27 mkim - changed the way to select all-sky GMI observations to use in BC coefficients updates +! 2024-03-27 mkim - changed the way to select all-sky GMI observations to use in BC coefficients updates +! 2026-05-22 a.lee - add viirs-m radiance assimilation ! ! input argument list: ! lunin - unit from which to read radiance (brightness temperature, tb) obs @@ -302,7 +303,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& use crtm_interface, only: ilzen_ang2,iscan_ang2,iszen_ang2,isazi_ang2, ilazi_ang2 use crtm_module, only: H2O_ID use clw_mod, only: calc_clw, ret_amsua, gmi_37pol_diff,mhs_si - use qcmod, only: qc_ssmi,qc_seviri,qc_abi,qc_ssu,qc_avhrr,qc_goesimg,qc_msu,qc_irsnd,qc_amsua,qc_mhs,qc_atms + use qcmod, only: qc_ssmi,qc_seviri,qc_abi,qc_ssu,qc_avhrr,qc_viirs,qc_goesimg,qc_msu,qc_irsnd,qc_amsua,qc_mhs,qc_atms use qcmod, only: igood_qc,ifail_gross_qc,ifail_interchan_qc,ifail_crtm_qc,ifail_satinfo_qc,qc_noirjaco3,ifail_cloud_qc use qcmod, only: ifail_cao_qc,cao_check use qcmod, only: ifail_iland_det, ifail_isnow_det, ifail_iice_det, ifail_iwater_det, ifail_imix_det, & @@ -379,7 +380,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& logical cao_flag logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb,goes_img,ahi,mhs,abi type(sparr2) :: dhx_dx - logical avhrr,avhrr_navy,lextra,ssu,iasi,cris,seviri,atms + logical avhrr,avhrr_navy,viirs,lextra,ssu,iasi,cris,seviri,atms logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,amsr2,gmi,saphir logical tmi logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img @@ -522,6 +523,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ahi = obstype == 'ahi' avhrr = obstype == 'avhrr' avhrr_navy = obstype == 'avhrr_navy' + viirs = obstype == 'viirs-m' ssmi = obstype == 'ssmi' amsre_low = obstype == 'amsre_low' amsre_mid = obstype == 'amsre_mid' @@ -1539,6 +1541,27 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, & id_qc,aivals,errf,varinv,varinv_use,cld,cldp) +! +! ---------- VIIRS ------------------- +! VIIRS Q C + + else if (viirs) then + + frac_sea=data_s(ifrac_sea,n) + + do i=1,nchanl + m=ich(i) + if (icld_det(m)>0 .and. varinv(i) >= tiny_r_kind) then + varinv_use(i) = varinv(i) + else + varinv_use(i) = zero + end if + end do + + call qc_viirs(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n), & + zsges,cenlat,frac_sea,pangs,trop5,tzbgr,tsavg5,tbc,tb_obs,tnoise, & + wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, & + id_qc,aivals,errf,varinv,varinv_use,cld,cldp) ! ---------- SSM/I , SSMIS, AMSRE ------------------- ! SSM/I, SSMIS, & AMSRE Q C diff --git a/GEOSaana_GridComp/GSI_GridComp/sfcobsqc.f90 b/GEOSaana_GridComp/GSI_GridComp/sfcobsqc.f90 index 8327b53..e238f8d 100644 --- a/GEOSaana_GridComp/GSI_GridComp/sfcobsqc.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/sfcobsqc.f90 @@ -22,6 +22,7 @@ module sfcobsqc ! 2014-10-01 Xue - add GSD surface data uselist ! 2015-07-10 pondeca - add reject list for cldch ! 2018-03-14 pondeca - add station accept list for mesonet visibility +! 2026-05-22 a.lee - correct the use of calcsun in subroutine get_sunangle ! ! subroutines included: ! sub init_rjlists @@ -1126,7 +1127,8 @@ subroutine get_sunangle(idate,dtime,dlon,dlat,sunangle) rlon=real(dlon/deg2rad,kind=r_single) ; if (rlon > s180) rlon=rlon-s360 rlat=real(dlat/deg2rad,kind=r_single) - call calcsun(idayyr3,ndays,time,rlat,rlon,sunangle) + call calcsun(idayyr3,ndays,time,rlon,rlat,sunangle) + !call calcsun(idayyr3,ndays,time,rlat,rlon,sunangle) end subroutine get_sunangle