!---------------------------------------------------------------------------- !BOP ! ! !MODULE: GSI_NSTCouplerMod --- ! ! !DESCRIPTION: This stub provides the default interfaces to the ! NST analysis calculations in GSI. ! ! !REVISION HISTORY: ! ! 07Oct2011 Akella/RT - Initial code ! 05Mar2012 Akella - Create_nst and getnst from satthin are now nst_int_ & nst_set_ ! Destroy_nst from satthin is nst_final_ ! 2017--9-27 Li - change nst_set to be nst_read and fix a bug in cal_tztr_ ! !EOP !------------------------------------------------------------------------- subroutine nst_init_() use mpimod, only: mype use gridmod, only: nlat_sfc,nlon_sfc, nlat, nlon use guess_grids, only: nfldnst, ntguesnst use gsi_nstcouplermod, only: tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,& c_0_full,c_d_full,w_0_full,w_d_full use mpeu_util, only: die, perr implicit none if(.not.allocated(tref_full)) allocate(tref_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(dt_cool_full)) allocate(dt_cool_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(z_c_full)) allocate(z_c_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(dt_warm_full)) allocate(dt_warm_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(z_w_full)) allocate(z_w_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(c_0_full)) allocate(c_0_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(c_d_full)) allocate(c_d_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(w_0_full)) allocate(w_0_full (nlat_sfc,nlon_sfc,nfldnst)) if(.not.allocated(w_d_full)) allocate(w_d_full (nlat_sfc,nlon_sfc,nfldnst)) if( ntguesnst < 1 .or. ntguesnst > nfldnst ) then call perr('nst_init','ntguesnst = ',ntguesnst) call die('nst_init') endif if(mype == 0)write(6,*)'GETNST: enter with nlat_sfc,nlon_sfc=',nlat_sfc,nlon_sfc,& ' and nlat,nlon=',nlat,nlon end subroutine nst_init_ !******************************************************************************************* subroutine nst_read_(mype_io) use kinds, only: i_kind use gridmod, only: use_gfs_nemsio, use_gfs_ncio use ncepgfs_io, only: read_gfsnst use netcdfgfs_io, only: read_gfsncnst use ncepnems_io, only: read_nemsnst use gsi_nstcouplermod, only: tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,& c_0_full,c_d_full,w_0_full,w_d_full implicit none integer(i_kind), intent(in ) :: mype_io if ( use_gfs_nemsio ) then call read_nemsnst(mype_io,tref_full,dt_cool_full,z_c_full, & dt_warm_full,z_w_full,c_0_full,c_d_full,w_0_full,w_d_full) else if ( use_gfs_ncio ) then call read_gfsncnst(mype_io,tref_full,dt_cool_full,z_c_full, & dt_warm_full,z_w_full,c_0_full,c_d_full,w_0_full,w_d_full) else call read_gfsnst(mype_io,tref_full,dt_cool_full,z_c_full, & dt_warm_full,z_w_full,c_0_full,c_d_full,w_0_full,w_d_full) endif end subroutine nst_read_ !******************************************************************************************* subroutine nst_final_ () use gsi_nstcouplermod, only: tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,& c_0_full,c_d_full,w_0_full,w_d_full if(allocated(tref_full)) deallocate(tref_full) if(allocated(dt_cool_full)) deallocate(dt_cool_full) if(allocated(z_c_full)) deallocate(z_c_full) if(allocated(dt_warm_full)) deallocate(dt_warm_full) if(allocated(z_w_full)) deallocate(z_w_full) if(allocated(c_0_full)) deallocate(c_0_full) if(allocated(c_d_full)) deallocate(c_d_full) if(allocated(w_0_full)) deallocate(w_0_full) if(allocated(w_d_full)) deallocate(w_d_full) end subroutine nst_final_ !******************************************************************************************* subroutine deter_nst_(dlat_earth,dlon_earth,obstime,zob,tref,dtw,dtc,tz_tr) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_nst 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: ! 2011-04-08 Li ! 2013-01-23 parrish - change from grdcrd to grdcrd1 (to allow successful ! ! input argument list: ! obstime - observation time relative to analysis time ! dlat_earth - earth latitude in radians ! dlon_earth - earth longitude in radians ! 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(tr) ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use constants, only: zero,one,z_w_max use gridmod, only: regional,tll2xy,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc use guess_grids, only: nfldnst,hrdifnst use gsi_nstcouplermod, only: fac_dtl,fac_tsl use gsi_nstcouplermod, only: tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,& c_0_full,c_d_full,w_0_full,w_d_full use satthin, only: isli_full 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 real(r_kind):: dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d integer(i_kind) istyp00,istyp01,istyp10,istyp11 integer(i_kind):: itnst,itnstp integer(i_kind):: ix,iy,ixp,iyp,j real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11,dtnst,dtnstp 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 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):: c_0_00,c_0_01,c_0_10,c_0_11 real(r_kind):: c_d_00,c_d_01,c_d_10,c_d_11 real(r_kind):: w_0_00,w_0_01,w_0_10,w_0_11 real(r_kind):: w_d_00,w_d_01,w_d_10,w_d_11 real(r_kind):: wgtavg,dlat,dlon logical outside 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 > hrdifnst(1) .and. obstime <= hrdifnst(nfldnst))then do j=1,nfldnst-1 if(obstime > hrdifnst(j) .and. obstime <= hrdifnst(j+1))then itnst=j itnstp=j+1 dtnst=(hrdifnst(j+1)-obstime)/(hrdifnst(j+1)-hrdifnst(j)) end if end do else if(obstime <=hrdifnst(1))then itnst=1 itnstp=1 dtnst=one else itnst=nfldnst itnstp=nfldnst 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 = 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 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 c_0_00 = c_0_full (ix ,iy ,itnst)*dtnst + c_0_full (ix ,iy ,itnstp)*dtnstp c_0_01 = c_0_full (ix ,iyp,itnst)*dtnst + c_0_full (ix ,iyp,itnstp)*dtnstp c_0_10 = c_0_full (ixp,iy ,itnst)*dtnst + c_0_full (ixp,iy ,itnstp)*dtnstp c_0_11 = c_0_full (ixp,iyp,itnst)*dtnst + c_0_full (ixp,iyp,itnstp)*dtnstp c_d_00 = c_d_full (ix ,iy ,itnst)*dtnst + c_d_full (ix ,iy ,itnstp)*dtnstp c_d_01 = c_d_full (ix ,iyp,itnst)*dtnst + c_d_full (ix ,iyp,itnstp)*dtnstp c_d_10 = c_d_full (ixp,iy ,itnst)*dtnst + c_d_full (ixp,iy ,itnstp)*dtnstp c_d_11 = c_d_full (ixp,iyp,itnst)*dtnst + c_d_full (ixp,iyp,itnstp)*dtnstp w_0_00 = w_0_full (ix ,iy ,itnst)*dtnst + w_0_full (ix ,iy ,itnstp)*dtnstp w_0_01 = w_0_full (ix ,iyp,itnst)*dtnst + w_0_full (ix ,iyp,itnstp)*dtnstp w_0_10 = w_0_full (ixp,iy ,itnst)*dtnst + w_0_full (ixp,iy ,itnstp)*dtnstp w_0_11 = w_0_full (ixp,iyp,itnst)*dtnst + w_0_full (ixp,iyp,itnstp)*dtnstp w_d_00 = w_d_full (ix ,iy ,itnst)*dtnst + w_d_full (ix ,iy ,itnstp)*dtnstp w_d_01 = w_d_full (ix ,iyp,itnst)*dtnst + w_d_full (ix ,iyp,itnstp)*dtnstp w_d_10 = w_d_full (ixp,iy ,itnst)*dtnst + w_d_full (ixp,iy ,itnstp)*dtnstp w_d_11 = w_d_full (ixp,iyp,itnst)*dtnst + w_d_full (ixp,iyp,itnstp)*dtnstp ! Interpolate nst variables to obs location (water surface only) wgtavg = zero tr_tmp = zero dt_cool = zero z_c = zero dt_warm = zero z_w_tmp = zero c_0 = zero c_d = zero w_0 = zero w_d = 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 = z_c + w00*z_c_00 dt_warm = dt_warm + w00*dt_warm_00 z_w_tmp = z_w_tmp + w00*z_w_00 c_0 = c_0 + w00*c_0_00 c_d = c_d + w00*c_d_00 w_0 = w_0 + w00*w_0_00 w_d = w_d + w00*w_d_00 endif if(istyp01 == 0)then wgtavg = wgtavg + w01 tr_tmp = tr_tmp + w01*tref_01 dt_cool = dt_cool + w01*dt_cool_01 z_c = z_c + w01*z_c_01 dt_warm = dt_warm + w01*dt_warm_01 z_w_tmp = z_w_tmp + w01*z_w_01 c_0 = c_0 + w01*c_0_01 c_d = c_d + w01*c_d_01 w_0 = w_0 + w01*w_0_01 w_d = w_d + w01*w_d_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 = z_c + w10*z_c_10 dt_warm = dt_warm + w10*dt_warm_10 z_w_tmp = z_w_tmp + w10*z_w_10 c_0 = c_0 + w10*c_0_10 c_d = c_d + w10*c_d_10 w_0 = w_0 + w10*w_0_10 w_d = w_d + w10*w_d_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 = z_c + w11*z_c_11 dt_warm = dt_warm + w11*dt_warm_11 z_w_tmp = z_w_tmp + w11*z_w_11 c_0 = c_0 + w11*c_0_11 c_d = c_d + w11*c_d_11 w_0 = w_0 + w11*w_0_11 w_d = w_d + w11*w_d_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_cool = dt_cool/wgtavg z_c = z_c/wgtavg dt_warm = dt_warm/wgtavg c_0 = c_0/wgtavg c_d = c_d/wgtavg w_0 = w_0/wgtavg w_d = w_d/wgtavg 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,c_0,c_d,w_0,w_d,z_c,z_w,zob,tz_tr) end if end subroutine deter_nst_ !******************************************************************************************* subroutine cal_tztr_(dt_warm,c_0,c_d,w_0,w_d,zc,zw,z,tztr) ! ! abstract: calculate d(Tz)/d(Ts) with T-Profile info from NSST Model ! ! prgmmr: li, xu org: np23 date: 2011-04-08 ! input variables ! ! dt_warm : diurnal warming amount at the surface ! xz : DTL depth (M) ! c_0 : coefficint 1 to calculate d(Tc)/d(Ts) ! c_d : coefficint 2 to calculate d(Tc)/d(Ts) ! w_0 : coefficint 1 to calculate d(Tw)/d(Ts) ! w_d : coefficint 2 to calculate d(Tw)/d(Ts) ! ! output variables ! ! tztr : d(Tz)/d(Tr) use kinds, only: r_kind use constants, only: one,two,half,zero use gsi_nstcouplermod, only: fac_dtl,fac_tsl real(kind=r_kind), intent(in) :: dt_warm,c_0,c_d,w_0,w_d,zc,zw,z real(kind=r_kind), intent(out) :: tztr ! local variables real(kind=r_kind) :: c1,c2,c3 c1 = one-two*(fac_dtl*w_0-fac_tsl*c_0)-(fac_dtl*w_d-fac_tsl*c_d)*z c2 = one-two*(fac_dtl*w_0-fac_tsl*c_0)-fac_dtl*w_d*z c3 = one+fac_tsl*two*c_0+fac_dtl*c_d*z tztr = one if ( dt_warm > zero ) then if ( z <= zc .and. c1 /= zero ) then tztr = (one-fac_dtl*w_0+fac_tsl*c_0)/c1 elseif ( z > zc .and. z < zw .and. c2 /= zero ) then tztr = (one-fac_dtl*w_0+fac_tsl*c_0)/c2 endif elseif ( dt_warm == zero .and. c3 /= zero ) then if ( z <= zc ) then tztr = (one+fac_tsl*c_0)/c3 endif endif if ( tztr <= -1.0_r_kind .or. tztr > 4.0_r_kind ) then write(6,100) fac_dtl,fac_tsl,c1,c2,c3,dt_warm,c_0,c_d,w_0,w_d,zc,zw,z,tztr 100 format('CAL_TZTR compute ',2(i2,1x),12(g13.6,1x),' RESET tztr to 1.0') tztr = one endif end subroutine cal_tztr_ !******************************************************************************************* !******************************************************************************************* subroutine skindepth_(obstype,sd_rad) ! ! abstract: Get skin depth (instrument dependent). Ideally, a skin-depth model calculates the channel dependent sd ! ! program history log: ! 2011-04-08 li use kinds, only: r_kind implicit none character(10), intent(in) :: obstype real(kind=r_kind), intent(out) :: sd_rad sd_rad = 0.000015_r_kind if ( obstype == 'amsre' .or. obstype == 'amsr2' .or. obstype == 'gmi' ) then sd_rad = 0.03_r_kind elseif ( obstype == 'amsua' .or. obstype == 'amsub' .or. obstype == 'ssmis' .or. obstype == 'ssmi' .or. & obstype == 'mhs' .or. obstype == 'msu' .or. obstype == 'hsb' ) then sd_rad = 0.001_r_kind endif end subroutine skindepth_ !*******************************************************************************************