!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !MODULE: guess_grids --- Guess-related grid definitions
!
! !INTERFACE:
!
module guess_grids
! !USES:
use kinds, only: r_single,r_kind,i_kind
use constants, only: max_varname_length
use gridmod, only: regional
use gridmod, only: wrf_nmm_regional,nems_nmmb_regional,fv3_regional
use gridmod, only: eta1_ll
use gridmod, only: eta2_ll
use gridmod, only: aeta1_ll
use gridmod, only: aeta2_ll
use gridmod, only: pdtop_ll
use gridmod, only: pt_ll
use gsi_bundlemod, only : gsi_bundlegetpointer
! meteorological guess (beyond standard ones)
use gsi_metguess_mod, only: gsi_metguess_create_grids
use gsi_metguess_mod, only: gsi_metguess_destroy_grids
use gsi_metguess_mod, only: gsi_metguess_get
use gsi_metguess_mod, only: gsi_metguess_bundle
! chem trace gases
use gsi_chemguess_mod, only: gsi_chemguess_create_grids
use gsi_chemguess_mod, only: gsi_chemguess_destroy_grids
use gsi_chemguess_mod, only: gsi_chemguess_get
! derivatives
use derivsmod, only: create_ges_derivatives
use derivsmod, only: destroy_ges_derivatives
! tendencies
use tendsmod, only: create_ges_tendencies
use tendsmod, only: destroy_ges_tendencies
use mpeu_util, only: die,tell
implicit none
! !DESCRIPTION: module containing variables related to the guess fields
!
! !REVISION HISTORY:
!
! 2003-12-01 kleist
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue; added onlys
! 2005-06-01 treadon - add routine add_rtm_layers
! 2005-06-03 parrish - add horizontal derivatives of guess fields
! 2005-06-10 devenyi/treadon - initialize nfldsig and nfldsfc
! 2005-08-03 parrish - add array to hold roughness length
! 2005-09-29 kleist - add derivatives of terrain, move prsi allocation
! 2005-11-21 kleist - add tendency arrays
! 2005-11-29 derber - add ozmz remove psfcg
! 2005-11-30 derber - combine create_atm_grids and create_pcp_grids (and destroys)
! 2006-02-02 treadon - prefix prsi,prsl,lnprsl,prslk with "ges_"
! 2006-03-07 treadon - remove ges_prslk (no longer needed)
! 2006-04-14 treadon - add bias_tskin
! 2006-04-17 treadon - add ges_psfcavg and ges_prslavg
! 2006-04-21 kleist - add ges divt and agvt arrays
! 2006-07-28 derber - clean up add_rtm_layers routine
! 2006-07-28 derber - add ges_tsen (sensible temperature) array
! 2006-07-31 kleist - use ges_ps instead of ln(ps)
! 2006-09-20 cucurull - add ln(ges_prsi) array
! 2006-09-29 treadon - add flags to control 10m wind factor recompute
! 2007-05-30 h.liu - remove ozmz
! 2007-06-21 rancic - add pbl (ges_teta)
! 2006-12-01 todling - remove bias stuff; merging GMAO bias correction scheme
! 2006-12-15 todling - add _initialized parameters to control allocations
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2008-02-07 eliu - fixed the unit difference between prsitmp
! (kPa) and toa_pressure (hPa).
! 2009-08-19 guo - added sfc_grids_allocated_, ges_grids_allocated_,
! and gesfinfo_created_ to track the state of the data.
! for multi-pass observer.
! - merged destroy_sfc_grids() and destroy_sfct().
! 2008-08-25 hu - add array definitions for hydrometeor fields
! - add subroutine create_cld_grids and destroy_cld_grids
! 2010-04-16 hou - add array definitions ges_co2 (co2 mixing ratio) and
! control variable igfsco2
! 2010-04-22 todling - remove tracers,vtid,pdryini,xncld
! 2010-05-19 todling - add chem init and destroy (revamp Hou's implementation)
! 2010-08-31 cucurull - add logical use_compress
! 2010-09-15 pagowski - add cmaq
! 2010-12-20 cucurull - add integer nsig_ext
! 2011-01-05 cucurull - add real gpstop
! 2011-02-11 zhu - add ges_gust,ges_vis,ges_pblh
! 2011-03-13 li - add for nst FCST file
! 2011-04-29 todling - some of cloud fields move to wrf_guess_mod; some to met_guess
! 2011-05-01 todling - cwmr no longer in guess-grids; use metguess bundle now
! 2011-11-01 eliu - modified condition to allocate/deallocate arrays related to
! cloud water tendencies and derivatives
! 2011-12-27 kleist - add 4d guess array for saturation specific humidity
! 2012-01-11 Hu - add GSD PBL height
! 2013-02-22 Carley - Add NMMB to GSD PBL height calc
! 2013-10-19 todling - metguess now holds background
! all tendencies now in a bundle (see tendsmod)
! all derivaties now in a bundle (see derivsmod)
! 2015-01-15 Hu - Add coast_prox to hold coast proximity
! 2017-05-12 Y. Wang and X. Wang - add bottom and top levels of w and rho for
! radar DA later, POC: xuguang.wang@ou.edu
! 2017-10-10 wu - Add code for fv3_regional
!
! !AUTHOR:
! kleist org: np20 date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------
! set default to private
private
! set subroutines to public
public :: create_sfc_grids
public :: create_ges_grids
public :: destroy_ges_grids
public :: destroy_sfc_grids
public :: create_gesfinfo
public :: destroy_gesfinfo
public :: load_prsges
public :: load_geop_hgt
public :: load_gsdpbl_hgt
public :: add_rtm_layers
public :: load_fact10
public :: comp_fact10
public :: guess_grids_print
public :: guess_grids_stats
public :: create_metguess_grids
public :: destroy_metguess_grids
public :: create_chemges_grids
public :: destroy_chemges_grids
! set passed variables to public
public :: ntguessig,ges_prsi,ges_psfcavg,ges_prslavg
public :: isli2,ges_prsl,nfldsig
public :: ges_teta
public :: fact_tv,tropprs,sfct
public :: ntguessfc,ntguesnst,dsfct,ifilesig,veg_frac,soil_type,veg_type
public :: sno2,ifilesfc,ifilenst,sfc_rough,fact10,sno,isli,soil_temp,soil_moi,coast_prox
public :: nfldsfc,nfldnst,hrdifsig,ges_tsen,sfcmod_mm5,sfcmod_gfs,ifact10,hrdifsfc,hrdifnst
public :: geop_hgti,ges_lnprsi,ges_lnprsl,geop_hgtl,pt_ll,pbl_height,ges_geopi
public :: wgt_lcbas
public :: ges_qsat
public :: use_compress,nsig_ext,gpstop
public :: ges_initialized
public :: nfldsig_all,nfldsig_now,hrdifsig_all
public :: nfldsfc_all,nfldsfc_now,hrdifsfc_all
public :: nfldnst_all,nfldnst_now,hrdifnst_all
public :: extrap_intime
public :: ntguessig_ref
public :: ntguessfc_ref
public :: ntguesnst_ref
public :: ges_w_btlev
public :: ges_rho
logical:: sfcmod_gfs = .false. ! .true. = recompute 10m wind factor using gfs physics
logical:: sfcmod_mm5 = .false. ! .true. = recompute 10m wind factor using mm5 physics
logical:: use_compress = .false. ! true to turn on compressibility factor in geopotential heights
logical, save :: ges_initialized = .false.
integer(i_kind) ntguessig ! location of actual guess time for sigma fields
integer(i_kind) ntguessfc ! location of actual guess time for sfc fields
integer(i_kind) ntguesnst ! location of actual guess time for nst FCST fields
integer(i_kind), save:: ntguessig_ref ! replace ntguessig as the storage for its original value
integer(i_kind), save:: ntguessfc_ref ! replace ntguessfc as the storage for its original value
integer(i_kind), save:: ntguesnst_ref ! replace ntguesnst as the storage for its original value
integer(i_kind):: ifact10 = 0 ! 0 = use 10m wind factor from guess
integer(i_kind):: nsig_ext = 13 ! use 13 layers above model top to compute the bending angle for gpsro
! number of guess sigma/surface times are set in GSI_gridComp.rc
real(r_kind), allocatable, dimension(:), save:: hrdifsig_all ! a list of all times
real(r_kind), allocatable, dimension(:), save:: hrdifsfc_all ! a list of all times
real(r_kind), allocatable, dimension(:), save:: hrdifnst_all ! a list of all times
integer(i_kind), save:: nfldsig_all ! expected total count of time slots
integer(i_kind), save:: nfldsfc_all
integer(i_kind), save:: nfldnst_all
integer(i_kind), save:: nfldsig ! actual count of in-cache time slots
integer(i_kind), save:: nfldsfc
integer(i_kind), save:: nfldnst ! actual count of in-cache time slots for NST file
integer(i_kind), save:: nfldsig_now ! current count of filled time slots
integer(i_kind), save:: nfldsfc_now
integer(i_kind), save:: nfldnst_now
logical, save:: extrap_intime ! compute o-f interpolate within the time ranges of guess_grids,
! or also extrapolate outside the time ranges.
real(r_kind), allocatable, dimension(:):: hrdifsig ! times for cached sigma guess_grid
real(r_kind), allocatable, dimension(:):: hrdifsfc ! times for cached surface guess_grid
real(r_kind), allocatable, dimension(:):: hrdifnst ! times for cached nst guess_grid
integer(i_kind),allocatable, dimension(:)::ifilesfc ! array used to open the correct surface guess files
integer(i_kind),allocatable, dimension(:)::ifilesig ! array used to open the correct sigma guess files
integer(i_kind),allocatable, dimension(:)::ifilenst ! array used to open the correct nst guess files
integer(i_kind),allocatable,dimension(:,:,:):: isli ! snow/land/ice mask
integer(i_kind),allocatable,dimension(:,:,:):: isli_g ! isli on horiz/global grid
integer(i_kind),allocatable,dimension(:,:):: isli2 ! snow/land/ice mask at analysis time
real(r_kind),allocatable,dimension(:,:):: coast_prox ! coast proximity mask
real(r_kind),allocatable,dimension(:,:,:):: sno2 ! sno depth on subdomain
real(r_kind):: gpstop=30.0_r_kind ! maximum gpsro height used in km
! geometric height for ref, impact height for bnd
real(r_kind):: ges_psfcavg ! average guess surface pressure
real(r_kind),allocatable,dimension(:):: ges_prslavg ! average guess pressure profile
real(r_kind),allocatable,dimension(:,:):: tropprs ! guess tropopause pressure
real(r_kind),allocatable,dimension(:,:,:):: dsfct ! delta skin temperature
real(r_kind),allocatable,dimension(:,:,:):: fact10 ! 10 meter wind factor
real(r_kind),allocatable,dimension(:,:,:):: sfct ! guess skin temperature
real(r_kind),allocatable,dimension(:,:,:):: sno ! snow-ice mask
real(r_kind),allocatable,dimension(:,:,:):: veg_type ! vegetation type
real(r_kind),allocatable,dimension(:,:,:):: veg_frac ! vegetation fraction(0-1.0)
real(r_kind),allocatable,dimension(:,:,:):: sfc_rough ! sfc roughness length
real(r_kind),allocatable,dimension(:,:,:):: soil_type ! soil type
real(r_kind),allocatable,dimension(:,:,:):: soil_temp ! soil temperature of first layer
real(r_kind),allocatable,dimension(:,:,:):: soil_moi ! soil moisture of first layer
real(r_kind),allocatable,dimension(:,:,:,:):: geop_hgtl ! guess geopotential height at mid-layers
real(r_kind),allocatable,dimension(:,:,:,:):: geop_hgti ! guess geopotential height at level interfaces
real(r_kind),allocatable,dimension(:,:,:,:):: ges_geopi ! input guess geopotential height at level interfaces
real(r_kind),allocatable,dimension(:,:,:):: pbl_height ! GSD PBL height in hPa
! Guess Fields ...
real(r_kind),allocatable,dimension(:,:):: wgt_lcbas ! weight given to base height of lowest cloud seen
real(r_kind),allocatable,dimension(:,:,:,:):: ges_prsi ! interface pressure
real(r_kind),allocatable,dimension(:,:,:,:):: ges_prsl ! layer midpoint pressure
real(r_kind),allocatable,dimension(:,:,:,:):: ges_lnprsl! log(layer midpoint pressure)
real(r_kind),allocatable,dimension(:,:,:,:):: ges_lnprsi! log(interface pressure)
real(r_kind),allocatable,dimension(:,:,:,:):: ges_tsen ! sensible temperature
real(r_kind),allocatable,dimension(:,:,:,:):: ges_teta ! potential temperature
real(r_kind),allocatable,dimension(:,:,:):: fact_tv ! 1./(one+fv*ges_q) for virt to sen calc.
real(r_kind),allocatable,dimension(:,:,:,:):: ges_qsat ! 4d qsat array
real(r_kind),allocatable,dimension(:,:,:,:):: ges_w_btlev
real(r_kind),allocatable,dimension(:,:,:,:):: ges_rho
interface guess_grids_print
module procedure print1r8_
module procedure print2r8_
module procedure print3r8_
module procedure print4r8_
end interface
interface guess_grids_stats
module procedure guess_grids_stats3d_
module procedure guess_grids_stats2d_
end interface
logical,save:: sfc_grids_allocated_=.false.
logical,save:: ges_grids_allocated_=.false.
logical,save:: gesfinfo_created_=.false.
character(len=*),parameter::myname='guess_grids'
contains
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_sfc_grids --- Allocate memory for surface related grids
!
! !INTERFACE:
!
subroutine create_sfc_grids
! !USES:
use gridmod, only: lat2,lon2,nlat,nlon
use constants, only: zero
implicit none
! !DESCRIPTION: allocate memory for surface related grids
!
! !REVISION HISTORY:
! 2003-12-01 kleist
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue
! 2004-07-28 treadon - remove subroutine call list, pass variables via modules
! 2005-06-03 parrish - allocate and initialize sfct_lat and sfct_lon
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2008-12-5 todling - add time dimension to dsfct
! 2009-01-23 todling - zero out arrays
! 2012-03-06 akella - add call to initialize arrays for NST analysis
! 2017-08-31 li - move gsi_nstcoupler_init & gsi_nstcoupler_final to
! satthin.F90 and read_obs.F90 respectivaly
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! kleist org: w/nmc20 date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind) :: i,j,it,istatus
if(sfc_grids_allocated_) call die('create_sfc_grids','alread allocated')
sfc_grids_allocated_=.true.
allocate( isli_g(nlat,nlon,nfldsfc),&
isli2(lat2,lon2),sno2(lat2,lon2,nfldsfc),&
stat=istatus)
if (istatus/=0) write(6,*)'CREATE_SFC_GRIDS(1): allocate error, istatus=',&
istatus,lat2,lon2,nlat,nlon,nfldsfc
#ifndef HAVE_ESMF
allocate( isli(lat2,lon2,nfldsfc),fact10(lat2,lon2,nfldsfc),&
dsfct(lat2,lon2,nfldsfc),sfct(lat2,lon2,nfldsfc),sno(lat2,lon2,nfldsfc),&
veg_type(lat2,lon2,nfldsfc),veg_frac(lat2,lon2,nfldsfc),&
sfc_rough(lat2,lon2,nfldsfc),&
soil_type(lat2,lon2,nfldsfc),soil_temp(lat2,lon2,nfldsfc),&
soil_moi(lat2,lon2,nfldsfc), coast_prox(lat2,lon2),&
stat=istatus)
if (istatus/=0) write(6,*)'CREATE_SFC_GRIDS(2): allocate error, istatus=',&
istatus,lat2,lon2,nlat,nlon,nfldsfc
#endif /* HAVE_ESMF */
do it=1,nfldsfc
do j=1,nlon
do i=1,nlat
isli_g(i,j,it)=0
end do
end do
end do
#ifndef HAVE_ESMF
do it=1,nfldsfc
do j=1,lon2
do i=1,lat2
isli(i,j,it)=0
coast_prox(i,j)=zero
fact10(i,j,it)=zero
sfct(i,j,it)=zero
dsfct(i,j,it)=zero
sno(i,j,it)=zero
veg_type(i,j,it)=zero
veg_frac(i,j,it)=zero
soil_type(i,j,it)=zero
soil_temp(i,j,it)=zero
soil_moi(i,j,it)=zero
end do
end do
end do
#endif
do it=1,nfldsfc
do j=1,lon2
do i=1,lat2
sno2(i,j,it)=zero
end do
end do
end do
do j=1,lon2
do i=1,lat2
isli2(i,j)=0
end do
end do
return
end subroutine create_sfc_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_ges_grids --- Alloc grid for guess and bias
!
! !INTERFACE:
!
subroutine create_ges_grids(switch_on_derivatives,tendsflag)
! !USES:
use control_vectors, only : w_exist
use constants,only: zero,one
use gridmod, only: lat2,lon2,nsig
implicit none
! !INPUT PARAMETERS:
logical,intent(in ) :: switch_on_derivatives ! for for horizontal derivatives
logical,intent(in ) :: tendsflag ! for time tendencies
! !OUTPUT PARAMETERS:
! !DESCRIPTION: allocate grids to hold guess and bias correction fields
!
! !REVISION HISTORY:
! 2004-06-03 treadon, original code
! 2004-07-15 todling, protex-compliant prologue; added onlys
! 2004-07-28 treadon - remove subroutine call list, pass variables via modules
! 2005-06-03 parrish - allocate/initialize _lat,_lon derivatives for u,v,cwmr,oz
! 2005-06-08 treadon - pass switch_on_derivatives via argument list
! 2005-07-27 kleist - modified to include some shared arrays
! 2006-01-10 treadon - remove mype from calling list (not used)
! 2006-07-31 kleist - use ges_ps arrays instead of ln(ps)
! 2006-12-04 todling - remove bias initialization; rename routine
! 2006-12-15 todling - protection to allow initializing ges/tnd/drv at will
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2011-02-09 zhu - add ges_gust,ges_vis,ges_pblh
! 2012-05-14 todling - revisit cw check to check also on some hydrometeors
! 2013-10-19 todling - revisit initialization of certain vars wrt ESMF
! 2014-06-09 carley/zhu - add wgt_lcbas
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! treadon org: w/nmc20 date: 2004-06-03
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind) i,j,k,n,istatus
if(ges_grids_allocated_) call die('create_ges_grids','already allocated')
ges_grids_allocated_=.true.
if ( .not. ges_initialized ) then
#ifndef HAVE_ESMF
nfldsig_all=nfldsig
nfldsfc_all=nfldsfc
nfldnst_all=nfldnst
nfldsig_now=0 ! _now variables are not used if not for ESMF
nfldsfc_now=0
nfldnst_now=0
extrap_intime=.true.
#endif /* HAVE_ESMF */
! Allocate and zero guess grids
allocate ( ges_prsi(lat2,lon2,nsig+1,nfldsig),ges_prsl(lat2,lon2,nsig,nfldsig),&
ges_lnprsl(lat2,lon2,nsig,nfldsig),ges_lnprsi(lat2,lon2,nsig+1,nfldsig),&
ges_tsen(lat2,lon2,nsig,nfldsig),&
ges_teta(lat2,lon2,nsig,nfldsig),&
ges_rho(lat2,lon2,nsig,nfldsig), &
geop_hgtl(lat2,lon2,nsig,nfldsig), &
geop_hgti(lat2,lon2,nsig+1,nfldsig),ges_prslavg(nsig),&
ges_geopi(lat2,lon2,nsig+1,nfldsig),&
tropprs(lat2,lon2),fact_tv(lat2,lon2,nsig),&
pbl_height(lat2,lon2,nfldsig),wgt_lcbas(lat2,lon2), &
ges_qsat(lat2,lon2,nsig,nfldsig),stat=istatus)
if(w_exist)then
allocate(ges_w_btlev(lat2,lon2,2,nfldsig),stat=istatus)
endif
if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_prsi,..): allocate error, istatus=',&
istatus,lat2,lon2,nsig,nfldsig
ges_initialized = .true.
! Default for ges_psfcavg
ges_psfcavg=zero
do i=1,nsig
ges_prslavg(i)=zero
end do
do j=1,lon2
do i=1,lat2
tropprs(i,j)=zero
end do
end do
do k=1,nsig
do j=1,lon2
do i=1,lat2
fact_tv(i,j,k)=one
end do
end do
end do
do n=1,nfldsig
do j=1,lon2
do i=1,lat2
pbl_height(i,j,n)=zero
end do
end do
end do
do n=1,nfldsig
do k=1,nsig
do j=1,lon2
do i=1,lat2
ges_prsl(i,j,k,n)=zero
ges_lnprsl(i,j,k,n)=zero
ges_rho(i,j,k,n)=zero
ges_qsat(i,j,k,n)=zero
ges_tsen(i,j,k,n)=zero
ges_teta(i,j,k,n)=zero
geop_hgtl(i,j,k,n)=zero
end do
end do
end do
do k=1,nsig+1
do j=1,lon2
do i=1,lat2
ges_prsi(i,j,k,n)=zero
ges_lnprsi(i,j,k,n)=zero
geop_hgti(i,j,k,n)=zero
ges_geopi(i,j,k,n)=zero
end do
end do
end do
end do
do j=1,lon2
do i=1,lat2
wgt_lcbas(i,j)=0.01_r_kind
end do
end do
if(w_exist) then
ges_w_btlev=zero
endif
end if ! ges_initialized
! If tendencies option on, allocate/initialize _ten arrays to zero
call create_ges_tendencies(tendsflag)
call create_ges_derivatives(switch_on_derivatives,nfldsig)
return
end subroutine create_ges_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_metguess_grids --- initialize meterological guess
!
! !INTERFACE:
!
subroutine create_metguess_grids(mype,istatus)
! !USES:
use gridmod, only: lat2,lon2,nsig
implicit none
! !INPUT PARAMETERS:
integer(i_kind), intent(in) :: mype
! !OUTPUT PARAMETERS:
integer(i_kind), intent(out) :: istatus
! !DESCRIPTION: initialize meteorological background fields beyond
! the standard ones - wired-in this module.
!
! !REVISION HISTORY:
! 2011-04-29 todling
! 2013-10-30 todling - update interface
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR:
! todling org: w/nmc20 date: 2011-04-29
!
!EOP
!-------------------------------------------------------------------------
character(len=*),parameter::myname_=myname//'*create_metguess_grids'
integer(i_kind) :: nmguess ! number of meteorol. fields (namelist)
character(len=max_varname_length),allocatable:: mguess(:) ! names of meterol. fields
istatus=0
! When proper connection to ESMF is complete,
! the following will not be needed here
! ------------------------------------------
call gsi_metguess_get('dim',nmguess,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble getting number of met-guess fields'
return
endif
if(nmguess==0) return
if (nmguess>0) then
allocate (mguess(nmguess))
call gsi_metguess_get('gsinames',mguess,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble getting name of met-guess fields'
return
endif
! Allocate memory for guess fields
! --------------------------------
call gsi_metguess_create_grids(lat2,lon2,nsig,nfldsig,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble allocating mem for met-guess'
return
endif
endif
end subroutine create_metguess_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_metguess_grids --- destroy meterological background
!
! !INTERFACE:
!
subroutine destroy_metguess_grids(mype,istatus)
! !USES:
implicit none
! !INPUT PARAMETERS:
integer(i_kind),intent(in)::mype
! !OUTPUT PARAMETERS:
integer(i_kind),intent(out)::istatus
! !DESCRIPTION: destroy meterological background
!
! !REVISION HISTORY:
! 2011-04-29 todling
! 2013-10-30 todling - update interface
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR:
! todling org: w/nmc20 date: 2011-04-29
!
!EOP
character(len=*),parameter::myname_=myname//'destroy_metguess_grids'
istatus=0
call gsi_metguess_destroy_grids(istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble deallocating mem for met-guess'
return
endif
end subroutine destroy_metguess_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_chemges_grids --- initialize chem component
!
! !INTERFACE:
!
subroutine create_chemges_grids(mype,istatus)
! !USES:
use gridmod, only: lat2,lon2,nsig
implicit none
! !INPUT PARAMETERS:
integer(i_kind), intent(in) :: mype
! !OUTPUT PARAMETERS:
integer(i_kind), intent(out) :: istatus
! !DESCRIPTION: initialize chem background
!
! !REVISION HISTORY:
! 2010-05-19 todling
! 2013-10-30 todling - update interface
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR:
! todling org: w/nmc20 date: 2010-05-19
!
!EOP
!-------------------------------------------------------------------------
character(len=*),parameter::myname_=myname//'*create_chemges_grids'
integer(i_kind) :: ntgases ! number of tracer gases (namelist)
character(len=max_varname_length),allocatable:: tgases(:) ! names of tracer gases
istatus=0
! When proper connection to ESMF is complete,
! the following will not be needed here
! ------------------------------------------
call gsi_chemguess_get('dim',ntgases,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble getting number of chem/gases'
return
endif
if(ntgases==0) return
if (ntgases>0) then
allocate (tgases(ntgases))
call gsi_chemguess_get('gsinames',tgases,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble getting name of chem/gases'
return
endif
! Allocate memory for guess files for trace gases
! ------------------------------------------------
call gsi_chemguess_create_grids(lat2,lon2,nsig,nfldsig,istatus)
if(istatus/=0) then
if(mype==0) write(6,*) myname_, ': trouble allocating mem for chem/gases'
return
endif
endif
end subroutine create_chemges_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_chemges_grids --- destroy chem background
!
! !INTERFACE:
!
subroutine destroy_chemges_grids(istatus)
! !USES:
implicit none
! !INPUT PARAMETERS:
! !OUTPUT PARAMETERS:
integer(i_kind),intent(out)::istatus
! !DESCRIPTION: destroy chem background
!
! !REVISION HISTORY:
! 2010-05-19 todling
! 2013-10-30 todling - update interface
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR:
! todling org: w/nmc20 date: 2010-05-19
!
!EOP
istatus=0
call gsi_chemguess_destroy_grids(istatus)
end subroutine destroy_chemges_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_ges_grids --- Dealloc guess and bias fields
!
! !INTERFACE:
!
subroutine destroy_ges_grids
! !USES:
use control_vectors, only : w_exist
implicit none
! !INPUT PARAMETERS:
! !DESCRIPTION: deallocate guess and bias grids
!
! !REVISION HISTORY:
! 2003-12-01 kleist
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue; added onlys
! 2005-06-03 parrish - deallocate _lat,_lon arrays for u,v,cwmr,oz
! 2005-06-08 treadon - check flag to see if need to deallocate derivatives
! 2005-07-27 kleist - modified to include some shared arrays
! 2006-07-31 kleist - use ges_ps arrays instead of ln(ps)
! 2006-12-04 todling - remove bias destroy; rename routine
! 2006-12-15 todling - using internal switches to deallc(tnds/drvs)
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2012-05-14 todling - revist cw check to check also on some hyrometeors
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! kleist org: w/nmc20 date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind):: istatus
call destroy_ges_derivatives
call destroy_ges_tendencies
!
deallocate(ges_prsi,ges_prsl,ges_lnprsl,ges_lnprsi,&
ges_tsen,ges_teta,geop_hgtl,geop_hgti,ges_geopi,ges_prslavg,ges_rho,&
tropprs,fact_tv,pbl_height,wgt_lcbas,ges_qsat,stat=istatus)
if(w_exist) deallocate(ges_w_btlev,stat=istatus)
if (istatus/=0) &
write(6,*)'DESTROY_GES_GRIDS(ges_prsi,..): deallocate error, istatus=',&
istatus
return
end subroutine destroy_ges_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_sfc_grids --- Deallocate surface fields
!
! !INTERFACE:
!
subroutine destroy_sfc_grids
! !USES:
implicit none
! !DESCRIPTION: deallocate surface related grids
!
! !REVISION HISTORY:
! 2003-12-01 kleist
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue
! 2005-06-03 parrish - deallocate sfct_lat and sfct_lon
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2008-06-30 derber - remove sfct deallocate to allow earlier call
! 2009-01-17 todling - move isli2,sno2 into destroy_sfct
! 2010-03-15 todling - esmf protection
! 2012-03-06 akella - add call to destroy NST analysis arrays
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! kleist org: w/nmc20 date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind):: istatus
if(.not.sfc_grids_allocated_) call die('destroy_sfc_grids_','not allocated')
sfc_grids_allocated_=.false.
deallocate(isli_g,isli2,sno2,stat=istatus)
if (istatus/=0) &
write(6,*)'DESTROY_SFC_GRIDS: deallocate error, istatus=',&
istatus
#ifndef HAVE_ESMF
if(allocated(isli))deallocate(isli)
if(allocated(fact10))deallocate(fact10)
if(allocated(sfct))deallocate(sfct)
if(allocated(sno))deallocate(sno)
if(allocated(veg_type))deallocate(veg_type)
if(allocated(veg_frac))deallocate(veg_frac)
if(allocated(soil_type))deallocate(soil_type)
if(allocated(sfc_rough))deallocate(sfc_rough)
if(allocated(soil_temp))deallocate(soil_temp)
if(allocated(soil_moi))deallocate(soil_moi)
if(allocated(dsfct))deallocate(dsfct)
if(allocated(coast_prox))deallocate(coast_prox)
#endif
return
end subroutine destroy_sfc_grids
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_gesfinfo --- Allocate guess-files information arrays
!
! !INTERFACE:
!
subroutine create_gesfinfo
! !USES:
implicit none
! !DESCRIPTION: allocate guess-files information arrays
!
! !REVISION HISTORY:
! 2009-01-08 todling
!
! !REMARKS:
! language: f90
! machine:ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! todling org: w/nmc2 date: 2009-01-08
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind):: istatus
if(gesfinfo_created_) call die('create_gesfinfo','already created')
gesfinfo_created_=.true.
#ifndef HAVE_ESMF
nfldsig_all=nfldsig
nfldsfc_all=nfldsfc
nfldnst_all=nfldnst
nfldsig_now=0 ! _now variables are not used if not for ESMF
nfldsfc_now=0
nfldnst_now=0
extrap_intime=.true.
allocate(hrdifsfc(nfldsfc),ifilesfc(nfldsfc), &
hrdifnst(nfldnst),ifilenst(nfldnst), &
hrdifsig(nfldsig),ifilesig(nfldsig), &
hrdifsfc_all(nfldsfc_all), &
hrdifnst_all(nfldnst_all), &
hrdifsig_all(nfldsig_all), &
stat=istatus)
if (istatus/=0) &
write(6,*)'CREATE_GESFINFO(hrdifsfc,..): allocate error, istatus=',&
istatus
#endif /* HAVE_ESMF */
return
end subroutine create_gesfinfo
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_gesfinfo --- Deallocate guess-files information
!
! !INTERFACE:
!
subroutine destroy_gesfinfo
! !USES:
implicit none
! !DESCRIPTION: deallocate guess-files information
!
! !REVISION HISTORY:
! 2009-01-08 todling
!
! !REMARKS:
! language: f90
! machine:ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! todling org: w/nmc2 date: 2009-01-08
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind):: istatus
if(.not.gesfinfo_created_) call die('destroy_gesfinfo','not created')
gesfinfo_created_=.false.
#ifndef HAVE_ESMF
deallocate(hrdifsfc,ifilesfc,hrdifnst,ifilenst,hrdifsig,ifilesig, &
hrdifsfc_all,hrdifnst_all,hrdifsig_all,stat=istatus)
if (istatus/=0) &
write(6,*)'DESTROY_GESFINFO: deallocate error, istatus=',&
istatus
nfldsfc_all=0
nfldnst_all=0
nfldsig_all=0
nfldsfc =0
nfldnst =0
nfldsig =0
#endif /* HAVE_ESMF */
return
end subroutine destroy_gesfinfo
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_prsges --- Populate guess pressure arrays
!
! !INTERFACE:
!
subroutine load_prsges
! !USES:
use constants,only: zero,one,rd_over_cp,one_tenth,half,ten,rd,r1000
use gridmod, only: lat2,lon2,nsig,ak5,bk5,ck5,tref5,idvc5,&
regional,wrf_nmm_regional,nems_nmmb_regional,wrf_mass_regional,&
cmaq_regional,pt_ll,aeta2_ll,fv3_regional,&
aeta1_ll,eta2_ll,pdtop_ll,eta1_ll,twodvar_regional,idsl5
use obsmod, only: dtype,ndat
implicit none
! !DESCRIPTION: populate guess pressure arrays
!
! !REVISION HISTORY:
! 2003-10-15 kleist
! 2004-03-22 parrish, regional capability added
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue; added onlys
! 2004-07-28 treadon - remove subroutine call list, pass variables via modules
! 2005-05-24 pondeca - add regional surface analysis option
! 2006-04-14 treadon - unify global calculations to use ak5,bk5
! 2006-04-17 treadon - add ges_psfcavg and ges_prslavg for regional
! 2006-07-31 kleist - use ges_ps instead of ln(ps)
! 2007-05-08 kleist - add fully generalized coordinate for pressure calculation
! 2011-07-07 todling - add cap for log(pressure) calculation
! 2017-03-23 Hu - add code to use hybrid vertical coodinate in WRF MASS
! core
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! kleist org: w/nmc20 date: 2003-10-15
!
!EOP
!-------------------------------------------------------------------------
! Declare local parameter
character(len=*),parameter::myname_=myname//'*load_prsges'
real(r_kind),parameter:: r1013=1013.0_r_kind
! Declare local variables
real(r_kind) kap1,kapr,trk
real(r_kind),dimension(:,:) ,pointer::ges_ps=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL()
real(r_kind) pinc(lat2,lon2)
integer(i_kind) i,j,k,ii,jj,itv,ips,kp
logical ihaveprs(nfldsig)
kap1=rd_over_cp+one
kapr=one/rd_over_cp
ihaveprs=.false.
do jj=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'ps' ,ges_ps,ips)
if(ips/=0) call die(myname_,': ps not available in guess, abort',ips)
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv,itv)
if(idvc5==3) then
if(itv/=0) call die(myname_,': tv must be present when idvc5=3, abort',itv)
endif
!!!!!!!!!!!! load delp to ges_prsi in read_fv3_netcdf_guess !!!!!!!!!!!!!!!!!
if (fv3_regional ) then
do j=1,lon2
do i=1,lat2
pinc(i,j)=(ges_ps(i,j)-ges_prsi(i,j,1,jj))
enddo
enddo
do k=1,nsig+1
do j=1,lon2
do i=1,lat2
ges_prsi(i,j,k,jj)=ges_prsi(i,j,k,jj)+eta2_ll(k)*pinc(i,j)
enddo
enddo
enddo
endif
do k=1,nsig+1
do j=1,lon2
do i=1,lat2
if(regional) then
if (wrf_nmm_regional.or.nems_nmmb_regional) &
ges_prsi(i,j,k,jj)=one_tenth* &
(eta1_ll(k)*pdtop_ll + &
eta2_ll(k)*(ten*ges_ps(i,j)-pdtop_ll-pt_ll) + &
pt_ll)
if (twodvar_regional .or. cmaq_regional ) &
ges_prsi(i,j,k,jj)=one_tenth*(eta1_ll(k)*(ten*ges_ps(i,j)-pt_ll) + pt_ll)
if (wrf_mass_regional) &
ges_prsi(i,j,k,jj)=one_tenth*(eta1_ll(k)*(ten*ges_ps(i,j)-pt_ll) + &
eta2_ll(k) + pt_ll)
else
if (idvc5==1 .or. idvc5==2) then
ges_prsi(i,j,k,jj)=ak5(k)+(bk5(k)*ges_ps(i,j))
else if (idvc5==3) then
if (k==1) then
ges_prsi(i,j,k,jj)=ges_ps(i,j)
else if (k==nsig+1) then
ges_prsi(i,j,k,jj)=zero
else
trk=(half*(ges_tv(i,j,k-1)+ges_tv(i,j,k))/tref5(k))**kapr
ges_prsi(i,j,k,jj)=ak5(k)+(bk5(k)*ges_ps(i,j))+(ck5(k)*trk)
end if
end if
endif
ges_prsi(i,j,k,jj)=max(ges_prsi(i,j,k,jj),zero)
ges_lnprsi(i,j,k,jj)=log(max(ges_prsi(i,j,k,jj),0.0001_r_kind))
end do
end do
end do
ihaveprs(jj)=.true.
end do
if(regional) then
if (wrf_nmm_regional.or.nems_nmmb_regional) then
! load using aeta coefficients
do jj=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'ps' ,ges_ps ,ips)
do k=1,nsig
do j=1,lon2
do i=1,lat2
ges_prsl(i,j,k,jj)=one_tenth* &
(aeta1_ll(k)*pdtop_ll + &
aeta2_ll(k)*(ten*ges_ps(i,j)-pdtop_ll-pt_ll) + &
pt_ll)
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
end if ! end if wrf_nmm regional block
if (fv3_regional) then
do jj=1,nfldsig
do k=1,nsig
kp=k+1
do j=1,lon2
do i=1,lat2
ges_prsl(i,j,k,jj)=(ges_prsi(i,j,k,jj)+ges_prsi(i,j,kp,jj))*half
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
end if ! end if fv3 regional
if (twodvar_regional .or. cmaq_regional) then
! load using aeta coefficients
do jj=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'ps' ,ges_ps ,ips)
do k=1,nsig
do j=1,lon2
do i=1,lat2
ges_prsl(i,j,k,jj)=one_tenth*(aeta1_ll(k)*(ten*ges_ps(i,j)-pt_ll)+pt_ll)
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
end if ! end if twodvar_regional, cmaq_regional block
if (wrf_mass_regional) then
! load using aeta coefficients
do jj=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'ps' ,ges_ps ,ips)
do k=1,nsig
do j=1,lon2
do i=1,lat2
ges_prsl(i,j,k,jj)=one_tenth*(aeta1_ll(k)*(ten*ges_ps(i,j)-pt_ll)+&
aeta2_ll(k) + pt_ll)
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
end if ! end if wrf_mass regional block
else
! load mid-layer pressure by using phillips vertical interpolation
if (idsl5/=2) then
do jj=1,nfldsig
if(.not.ihaveprs(jj)) then
call tell(myname,'3d pressure has not been calculated somehow',99)
exit ! won't die ...
endif
do j=1,lon2
do i=1,lat2
do k=1,nsig
ges_prsl(i,j,k,jj)=((ges_prsi(i,j,k,jj)**kap1-ges_prsi(i,j,k+1,jj)**kap1)/&
(kap1*(ges_prsi(i,j,k,jj)-ges_prsi(i,j,k+1,jj))))**kapr
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
! load mid-layer pressure by simple averaging
else
do jj=1,nfldsig
if(.not.ihaveprs(jj)) then
call tell(myname,'3d pressure has not been calculated somehow',99)
exit ! won't die ...
endif
do j=1,lon2
do i=1,lat2
do k=1,nsig
ges_prsl(i,j,k,jj)=(ges_prsi(i,j,k,jj)+ges_prsi(i,j,k+1,jj))*half
ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
end do
end do
end do
end do
endif
end if ! end regional/global block
! Compute density for dBZ assimilation purposes - multiply by 1000 to convert
! to Pa
do ii=1,ndat
if ( index(dtype(ii), 'dbz') /= 0 )then
do jj=1,nfldsig
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv,itv)
if(idvc5==3) then
if(itv/=0) call die(myname_,': tv must be present when idvc5=3,abort',itv)
endif
do j=1,lon2
do i=1,lat2
do k=1,nsig
ges_rho(i,j,k,jj)=(ges_prsl(i,j,k,jj)/(ges_tv(i,j,k)*rd))*r1000
end do
end do
end do
end do
end if
end do
! For regional applications only, load variables containing mean
! surface pressure and pressure profile at the layer midpoints
if (regional) then
ges_psfcavg = r1013
if (wrf_nmm_regional.or.nems_nmmb_regional) then
do k=1,nsig
ges_prslavg(k)=aeta1_ll(k)*pdtop_ll+aeta2_ll(k)*(r1013-pdtop_ll-pt_ll)+pt_ll
end do
elseif (fv3_regional) then
do k=1,nsig
ges_prslavg(k)=aeta1_ll(k)*ten+r1013*aeta2_ll(k)
end do
elseif (wrf_mass_regional) then
do k=1,nsig
ges_prslavg(k)=aeta1_ll(k)*(r1013-pt_ll)+aeta2_ll(k) + pt_ll
end do
else
do k=1,nsig
ges_prslavg(k)=aeta1_ll(k)*(r1013-pt_ll)+pt_ll
end do
endif
endif
return
end subroutine load_prsges
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_geop_hgt --- Populate guess geopotential height
!
! !INTERFACE:
!
subroutine load_geop_hgt
! !USES:
use constants, only: one,eps, rd, grav, half, t0c, fv
use constants, only: cpf_a0, cpf_a1, cpf_a2, cpf_b0, cpf_b1, cpf_c0, cpf_c1, cpf_d, cpf_e
use constants, only: psv_a, psv_b, psv_c, psv_d
use constants, only: ef_alpha, ef_beta, ef_gamma
use gridmod, only: lat2, lon2, nsig, twodvar_regional
implicit none
! !INPUT PARAMETERS:
! !DESCRIPTION: populate guess geopotential height
!
! !REVISION HISTORY:
! 2003-10-15 treadon
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-compliant prologue
! 2004-10-28 treadon - replace "tiny" with "tiny_r_kind"
! 2004-12-15 treadon - replace use of Paul van Delst's Geopotential
! function with simple integration of hydrostatic
! equation (done to be consistent with Lidia
! Cucurull's GPS work)
! 2005-05-24 pondeca - add regional surface analysis option
! 2010-08-27 cucurull - add option to compute and use compressibility factors in geopot heights
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! treadon org: w/nmc20 date: 2003-10-15
!
!EOP
!-------------------------------------------------------------------------
character(len=*),parameter::myname_=myname//'*load_geop_hgt'
real(r_kind),parameter:: thousand = 1000.0_r_kind
integer(i_kind) i,j,k,jj,ier,istatus
real(r_kind) h,dz,rdog
real(r_kind),dimension(nsig+1):: height
real(r_kind) cmpr, x_v, rl_hm, fact, pw, tmp_K, tmp_C, prs_sv, prs_a, ehn_fct, prs_v
real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_q=>NULL()
real(r_kind),dimension(:,: ),pointer::ges_z=>NULL()
if (twodvar_regional) return
rdog = rd/grav
if (use_compress) then
! Compute compressibility factor (Picard et al 2008) and geopotential heights at midpoint
! of each layer
do jj=1,nfldsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'z' ,ges_z ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'q' ,ges_q ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv ,istatus)
ier=ier+istatus
if(ier/=0) exit
do j=1,lon2
do i=1,lat2
k = 1
fact = one + fv * ges_q(i,j,k)
pw = eps + ges_q(i,j,k)*( one - eps )
tmp_K = ges_tv(i,j,k) / fact
tmp_C = tmp_K - t0c
prs_sv = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K) ! Pvap sat, eq A1.1 (Pa)
prs_a = thousand * exp(half*(log(ges_prsi(i,j,k,jj)) + log(ges_prsl(i,j,k,jj)))) ! (Pa)
ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
prs_v = ges_q(i,j,k) * prs_a / pw ! vapor pressure (Pa)
rl_hm = prs_v / prs_sv ! relative humidity
x_v = rl_hm * ehn_fct * prs_sv / prs_a ! molar fraction of water vapor (eq. A1.3)
! Compressibility factor (eq A1.4 from Picard et al 2008)
cmpr = one - (prs_a/tmp_K) * (cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
+ (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
+ (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
h = rdog * ges_tv(i,j,k)
dz = h * cmpr * log(ges_prsi(i,j,k,jj)/ges_prsl(i,j,k,jj))
height(k) = ges_z(i,j) + dz
do k=2,nsig
fact = one + fv * half * (ges_q(i,j,k-1)+ges_q(i,j,k))
pw = eps + half * (ges_q(i,j,k-1)+ges_q(i,j,k)) * (one - eps)
tmp_K = half * (ges_tv(i,j,k-1)+ges_tv(i,j,k)) / fact
tmp_C = tmp_K - t0c
prs_sv = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K) ! eq A1.1 (Pa)
prs_a = thousand * exp(half*(log(ges_prsl(i,j,k-1,jj))+log(ges_prsl(i,j,k,jj)))) ! (Pa)
ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
prs_v = half*(ges_q(i,j,k-1)+ges_q(i,j,k) ) * prs_a / pw ! (Pa)
rl_hm = prs_v / prs_sv ! relative humidity
x_v = rl_hm * ehn_fct * prs_sv / prs_a ! molar fraction of water vapor (eq. A1.3)
cmpr = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
+ (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
+ (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
h = rdog * half * (ges_tv(i,j,k-1)+ges_tv(i,j,k))
dz = h * cmpr * log(ges_prsl(i,j,k-1,jj)/ges_prsl(i,j,k,jj))
height(k) = height(k-1) + dz
end do
do k=1,nsig
geop_hgtl(i,j,k,jj)=height(k) - ges_z(i,j)
end do
enddo
enddo
enddo
if(ier/=0) return
! Compute compressibility factor (Picard et al 2008) and geopotential heights at interface
! between layers
do jj=1,nfldsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'z' ,ges_z ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'q' ,ges_q ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv ,istatus)
ier=ier+istatus
if(ier/=0) exit
do j=1,lon2
do i=1,lat2
k=1
height(k) = ges_z(i,j)
do k=2,nsig
fact = one + fv * ges_q(i,j,k-1)
pw = eps + ges_q(i,j,k-1)*(one - eps)
tmp_K = ges_tv(i,j,k-1) / fact
tmp_C = tmp_K - t0c
prs_sv = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K) ! eq A1.1 (Pa)
prs_a = thousand * exp(half*(log(ges_prsi(i,j,k-1,jj))+log(ges_prsi(i,j,k,jj))))
ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
prs_v = ges_q(i,j,k-1) * prs_a / pw ! vapor pressure (Pa)
rl_hm = prs_v / prs_sv ! relative humidity
x_v = rl_hm * ehn_fct * prs_sv / prs_a ! molar fraction of water vapor (eq. A1.3)
cmpr = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
+ (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
+ (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
h = rdog * ges_tv(i,j,k-1)
dz = h * cmpr * log(ges_prsi(i,j,k-1,jj)/ges_prsi(i,j,k,jj))
height(k) = height(k-1) + dz
enddo
k=nsig+1
fact = one + fv* ges_q(i,j,k-1)
pw = eps + ges_q(i,j,k-1)*(one - eps)
tmp_K = ges_tv(i,j,k-1) / fact
tmp_C = tmp_K - t0c
prs_sv = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K) ! eq A1.1 (Pa)
prs_a = thousand * exp(half*(log(ges_prsi(i,j,k-1,jj))+log(ges_prsl(i,j,k-1,jj)))) ! (Pa)
ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
prs_v = ges_q(i,j,k-1) * prs_a / pw
rl_hm = prs_v / prs_sv ! relative humidity
x_v = rl_hm * ehn_fct * prs_sv / prs_a ! molar fraction of water vapor (eq. A1.3)
cmpr = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
+ (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
+ (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
h = rdog * ges_tv(i,j,k-1)
dz = h * cmpr * log(ges_prsi(i,j,k-1,jj)/ges_prsl(i,j,k-1,jj))
height(k) = height(k-1) + dz
do k=1,nsig+1
geop_hgti(i,j,k,jj)=height(k) - ges_z(i,j)
end do
enddo
enddo
enddo
if(ier/=0) return
else
! Compute geopotential height at midpoint of each layer
do jj=1,nfldsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'z' ,ges_z ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv ,istatus)
ier=ier+istatus
if(ier/=0) exit
do j=1,lon2
do i=1,lat2
k = 1
h = rdog * ges_tv(i,j,k)
dz = h * log(ges_prsi(i,j,k,jj)/ges_prsl(i,j,k,jj))
height(k) = ges_z(i,j) + dz
do k=2,nsig
h = rdog * half * (ges_tv(i,j,k-1)+ges_tv(i,j,k))
dz = h * log(ges_prsl(i,j,k-1,jj)/ges_prsl(i,j,k,jj))
height(k) = height(k-1) + dz
end do
do k=1,nsig
geop_hgtl(i,j,k,jj)=height(k) - ges_z(i,j)
end do
end do
end do
end do
if(ier/=0) return
! Compute geopotential height at interface between layers
do jj=1,nfldsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'z' ,ges_z ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv ,istatus)
ier=ier+istatus
if(ier/=0) call die(myname_,'not all fields available, ier=',ier)
do j=1,lon2
do i=1,lat2
k=1
height(k) = ges_z(i,j)
do k=2,nsig
h = rdog * ges_tv(i,j,k-1)
dz = h * log(ges_prsi(i,j,k-1,jj)/ges_prsi(i,j,k,jj))
height(k) = height(k-1) + dz
end do
k=nsig+1
h = rdog * ges_tv(i,j,k-1)
dz = h * log(ges_prsi(i,j,k-1,jj)/ges_prsl(i,j,k-1,jj))
height(k) = height(k-1) + dz
do k=1,nsig+1
geop_hgti(i,j,k,jj)=height(k) - ges_z(i,j)
end do
end do
end do
end do
endif
return
end subroutine load_geop_hgt
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_gsdpbl_hgt --- Populate PBL height
!
! !INTERFACE:
!
subroutine load_gsdpbl_hgt(mype)
! !USES:
use constants, only: one,rd_over_cp_mass,r1000,ten,zero,two
use gridmod, only: lat2, lon2, nsig,wrf_mass_regional, &
twodvar_regional,nems_nmmb_regional,fv3_regional
implicit none
! !INPUT PARAMETERS:
! !DESCRIPTION: populate guess geopotential height in millibars
!
!
! !REVISION HISTORY:
! 2011-06-06 Ming Hu
! 2013-02-22 Jacob Carley - Added NMMB
!
! !REMARKS:
! language: f90
! machine: JET
!
! !AUTHOR:
!
!EOP
!-------------------------------------------------------------------------
character(len=*),parameter::myname_=myname//'*load_gsdpbl_hgt'
integer(i_kind) , intent(in ) :: mype
integer(i_kind) i,j,k,jj,ier,istatus
real(r_kind),dimension(nsig):: thetav
real(r_kind),dimension(nsig):: pbk
real(r_kind) :: thsfc, d
real(r_kind),dimension(:,: ),pointer::ges_ps_01=>NULL()
real(r_kind),dimension(:,: ),pointer::ges_ps=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL()
if (twodvar_regional) return
if (fv3_regional) then
if(mype==0)write(6,*)'not setup for fv3_regional in load_gsdpbl_hgt'
return
endif
! Compute geopotential height at midpoint of each layer
do jj=1,nfldsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(1) ,'ps' ,ges_ps_01 ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'ps' ,ges_ps ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(jj),'tv' ,ges_tv ,istatus)
ier=ier+istatus
if(ier/=0) call die(myname_,'not all fields available, ier=',ier)
do j=1,lon2
do i=1,lat2
do k=1,nsig
if (wrf_mass_regional) pbk(k) = aeta1_ll(k)*(ges_ps_01(i,j)*ten-pt_ll)+aeta2_ll(k)+pt_ll
if (nems_nmmb_regional) then
pbk(k) = aeta1_ll(k)*pdtop_ll + aeta2_ll(k)*(ten*ges_ps(i,j) &
-pdtop_ll-pt_ll) + pt_ll
end if
thetav(k) = ges_tv(i,j,k)*(r1000/pbk(k))**rd_over_cp_mass
end do
pbl_height(i,j,jj) = zero
thsfc = thetav(1)
k=1
DO while (abs(pbl_height(i,j,jj)) < 0.0001_r_kind)
if( thetav(k) > thsfc + 1.0_r_kind ) then
pbl_height(i,j,jj) = float(k) - (thetav(k) - (thsfc + 1.0_r_kind))/ &
max((thetav(k)-thetav(k-1)),0.01_r_kind)
endif
k=k+1
ENDDO
if(abs(pbl_height(i,j,jj)) < 0.0001_r_kind) pbl_height(i,j,jj)=two
k=int(pbl_height(i,j,jj))
if( k < 1 .or. k > nsig-1) then
write(6,*) ' Error in PBL height calculation ',mype,i,j,pbl_height(i,j,jj)
endif
d=pbl_height(i,j,jj) - k
pbl_height(i,j,jj) = pbk(k) * (one-d) + pbk(k+1) * d
end do
end do
end do
return
end subroutine load_gsdpbl_hgt
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: add_rtm_layers --- Add pressure layers for RTM use
!
! !INTERFACE:
!
subroutine add_rtm_layers(prsitmp,prsltmp,prsitmp_ext,prsltmp_ext,klevel)
! !USES:
use constants, only: half,ten,one_tenth
use gridmod, only: nsig,msig,nlayers
use crtm_module, only: toa_pressure
implicit none
! !INPUT PARAMETERS:
integer(i_kind),dimension(msig) ,intent( out) :: klevel
real(r_kind) ,dimension(nsig+1),intent(in ) :: prsitmp
real(r_kind) ,dimension(nsig) ,intent(in ) :: prsltmp
real(r_kind) ,dimension(msig+1),intent( out) :: prsitmp_ext
real(r_kind) ,dimension(msig) ,intent( out) :: prsltmp_ext
! !DESCRIPTION: Add pressure layers for use in RTM
!
! !REVISION HISTORY:
! 2005-06-01 treadon
! 2006-05-10 derber modify how levels are added above model top
! 2013-03-27 rancic fix for toa units: crtm(hPa); prsitmp(kPa)
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! treadon org: w/nmc20 date: 2005-06-01
!
!EOP
!-------------------------------------------------------------------------
! Declare local variables
integer(i_kind) k,kk,l
real(r_kind) dprs,toa_prs_kpa
! Convert toa_pressure to kPa
! ---------------------------
toa_prs_kpa = toa_pressure*one_tenth
! Check if model top pressure above rtm top pressure, where prsitmp
! is in kPa and toa_pressure is in hPa.
if (prsitmp(nsig) < toa_prs_kpa)then
write(6,*)'ADD_RTM_LAYERS: model top pressure(hPa)=', &
ten*prsitmp(nsig),&
' above rtm top pressure(hPa)=',toa_pressure
call stop2(35)
end if
! Linear in pressure sub-divsions
kk=0
do k = 1,nsig
if (nlayers(k)<=1) then
kk = kk + 1
prsltmp_ext(kk) = prsltmp(k)
prsitmp_ext(kk) = prsitmp(k)
klevel(kk) = k
else
if (k/=nsig) then
dprs = (prsitmp(k+1)-prsitmp(k))/nlayers(k)
else
dprs = (toa_prs_kpa -prsitmp(k))/nlayers(k)
end if
prsitmp_ext(kk+1) = prsitmp(k)
do l=1,nlayers(k)
kk=kk + 1
prsitmp_ext(kk+1) = prsitmp(k) + dprs*l
prsltmp_ext(kk) = half*(prsitmp_ext(kk+1)+prsitmp_ext(kk))
klevel(kk) = k
end do
endif
end do
! Set top of atmosphere pressure
prsitmp_ext(msig+1) = toa_prs_kpa
end subroutine add_rtm_layers
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_fact10 --- Compute 10m wind factor
!
! !INTERFACE:
!
subroutine load_fact10
! !USES:
use gridmod, only: lat2,lon2
implicit none
! !INPUT PARAMETERS:
! !DESCRIPTION: compute 10m wind factor
!
! !REVISION HISTORY:
! 2006-09-26 treadon
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! treadon org: w/nmc20 date: 2006-09-26
!
!EOP
!-------------------------------------------------------------------------
! Declare local variables
character(len=*),parameter::myname_=myname//'*load_fact10'
logical iqtflg
integer(i_kind):: i,j,it,itt,nt,regime,ier,istatus
integer(i_kind),dimension(nfldsfc):: indx
real(r_kind):: u10ges,v10ges,t2ges,q2ges
real(r_kind),dimension(:,: ),pointer::ges_ps=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_u=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_v=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_q=>NULL()
nt=0
indx=1
do i=1,nfldsfc
if(abs(hrdifsfc(i)-hrdifsig(i))<0.001_r_kind) then
nt=nt+1
indx(nt) = i
endif
end do
if (sfcmod_gfs) then
do it=1,nt
itt=indx(it)
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'ps',ges_ps,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'u' ,ges_u ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'v' ,ges_v ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'q' ,ges_q ,istatus)
ier=ier+istatus
if(ier/=0) call die(myname_,'not all fields available, ier=',ier)
do j=1,lon2
do i=1,lat2
call compute_fact10(ges_u(i,j,1),ges_v(i,j,1),&
ges_tsen(i,j,1,itt),ges_q(i,j,1),&
ges_ps(i,j),ges_prsi(i,j,1,itt), &
ges_prsi(i,j,2,itt),sfct(i,j,itt), &
sfc_rough(i,j,itt),isli(i,j,itt),fact10(i,j,itt))
end do
end do
end do
endif
if (sfcmod_mm5) then
iqtflg=.true.
do it=1,nt
itt=indx(it)
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'ps',ges_ps,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'u' ,ges_u ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'v' ,ges_v ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'tv',ges_tv,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itt),'q' ,ges_q ,istatus)
ier=ier+istatus
if(ier/=0) call die(myname_,'not all fields available, ier=',ier)
do j=1,lon2
do i=1,lat2
call SFC_WTQ_FWD (&
ges_ps(i,j),&
sfct(i,j,itt),&
ges_lnprsl(i,j,1,itt),&
ges_tv(i,j,1),&
ges_q(i,j,1),&
ges_u(i,j,1),&
ges_v(i,j,1),&
ges_lnprsl(i,j,2,itt),&
ges_tv(i,j,2),&
ges_q(i,j,2),&
geop_hgtl(i,j,1,itt),&
sfc_rough(i,j,itt),&
isli(i,j,itt),&
fact10(i,j,itt),&
u10ges,v10ges,t2ges,q2ges,regime,iqtflg)
end do
end do
end do
endif
return
end subroutine load_fact10
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: comp_fact10 --- compute 10m wind factor
!
! !INTERFACE:
!
subroutine comp_fact10(dlat,dlon,dtime,skint,sfcrough,islimsk,mype,factw)
! !USES:
use gridmod, only: nlat,nlon,&
lon1,istart,jstart
use constants, only: zero,one
implicit none
! !INPUT PARAMETERS:
real(r_kind) ,intent(in ) :: dlat,dlon,dtime,skint,sfcrough
real(r_kind) ,intent(inout) :: factw
integer(i_kind),intent(in ) :: mype,islimsk
! !DESCRIPTION: compute 10m wind factor
!
! !REVISION HISTORY:
! 2006-09-26 treadon
! 2008-12-05 todling - use dsfct(:,:,ntguessfc) for calculation
! 2015-03-10 todling - assign missing ps pointers
!
! !REMARKS:
! language: f90
! machine: ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
! treadon org: w/nmc20 date: 2006-09-26
!
!EOP
!-------------------------------------------------------------------------
! Declare local parameters
character(len=*),parameter::myname_=myname//'*comp_fact10'
! Declare local variables
logical iqtflg
integer(i_kind) ix,ix1,ixp,iy,iy1,iyp,regime
integer(i_kind) itsig,itsigp,j,m1,islimsk2,ier,istatus
real(r_kind) w00,w01,w10,w11
real(r_kind) delx,dely,delx1,dely1,dtsig,dtsigp
real(r_kind):: u10ges,v10ges,t2ges,q2ges
real(r_kind):: pgesin,ugesin,vgesin,qgesin,tgesin,prsigesin1
real(r_kind):: prsigesin2,lnpgesin1,lnpgesin2,tgesin2,qgesin2,geopgesin,ts
real(r_kind),dimension(:,: ),pointer::ges_ps_itsig=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_u_itsig=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_v_itsig=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_tv_itsig=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_q_itsig=>NULL()
real(r_kind),dimension(:,: ),pointer::ges_ps_itsigp=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_u_itsigp=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_v_itsigp=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_tv_itsigp=>NULL()
real(r_kind),dimension(:,:,:),pointer::ges_q_itsigp=>NULL()
islimsk2=islimsk
if(islimsk2 > 2)islimsk2=islimsk2-3
m1=mype+1
! Set spatial interpolation indices and weights
ix1=dlat
ix1=max(1,min(ix1,nlat))
delx=dlat-ix1
delx=max(zero,min(delx,one))
ix=ix1-istart(m1)+2
ixp=ix+1
if(ix1==nlat) then
ixp=ix
end if
delx1=one-delx
iy1=dlon
dely=dlon-iy1
iy=iy1-jstart(m1)+2
if(iy<1) then
iy1=iy1+nlon
iy=iy1-jstart(m1)+2
end if
if(iy>lon1+1) then
iy1=iy1-nlon
iy=iy1-jstart(m1)+2
end if
iyp=iy+1
dely1=one-dely
w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely
! Get time interpolation factors for sigma files
if(dtime > hrdifsig(1) .and. dtime < hrdifsig(nfldsig))then
do j=1,nfldsig-1
if(dtime > hrdifsig(j) .and. dtime <= hrdifsig(j+1))then
itsig=j
itsigp=j+1
dtsig=((hrdifsig(j+1)-dtime)/(hrdifsig(j+1)-hrdifsig(j)))
end if
end do
else if(dtime <=hrdifsig(1))then
itsig=1
itsigp=1
dtsig=one
else
itsig=nfldsig
itsigp=nfldsig
dtsig=one
end if
dtsigp=one-dtsig
ier=0
call gsi_bundlegetpointer(gsi_metguess_bundle(itsig) ,'ps' ,ges_ps_itsig ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'ps' ,ges_ps_itsigp,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsig) ,'u' ,ges_u_itsig ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'u' ,ges_u_itsigp,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsig) ,'v' ,ges_v_itsig ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'v' ,ges_v_itsigp,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsig) ,'tv',ges_tv_itsig ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'tv',ges_tv_itsigp,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsig) ,'q' ,ges_q_itsig ,istatus)
ier=ier+istatus
call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'q' ,ges_q_itsigp,istatus)
ier=ier+istatus
if(ier/=0) return
ts =(dsfct(ix,iy ,ntguessfc)*w00 + dsfct(ixp,iy ,ntguessfc)*w10 + &
dsfct(ix,iyp,ntguessfc)*w01 + dsfct(ixp,iyp,ntguessfc)*w11) + skint
pgesin=(ges_ps_itsig (ix,iy )*w00+ges_ps_itsig (ixp,iy )*w10+ &
ges_ps_itsig (ix,iyp)*w01+ges_ps_itsig (ixp,iyp)*w11)*dtsig + &
(ges_ps_itsigp(ix,iy )*w00+ges_ps_itsigp(ixp,iy )*w10+ &
ges_ps_itsigp(ix,iyp)*w01+ges_ps_itsigp(ixp,iyp)*w11)*dtsigp
ugesin=(ges_u_itsig (ix,iy ,1)*w00+ges_u_itsig (ixp,iy ,1)*w10+ &
ges_u_itsig (ix,iyp,1)*w01+ges_u_itsig (ixp,iyp,1)*w11)*dtsig + &
(ges_u_itsigp(ix,iy ,1)*w00+ges_u_itsigp(ixp,iy ,1)*w10+ &
ges_u_itsigp(ix,iyp,1)*w01+ges_u_itsigp(ixp,iyp,1)*w11)*dtsigp
vgesin=(ges_v_itsig (ix,iy ,1)*w00+ges_v_itsig (ixp,iy ,1)*w10+ &
ges_v_itsig (ix,iyp,1)*w01+ges_v_itsig (ixp,iyp,1)*w11)*dtsig + &
(ges_v_itsigp(ix,iy ,1)*w00+ges_v_itsigp(ixp,iy ,1)*w10+ &
ges_v_itsigp(ix,iyp,1)*w01+ges_v_itsigp(ixp,iyp,1)*w11)*dtsigp
qgesin=(ges_q_itsig (ix,iy ,1)*w00+ges_q_itsig (ixp,iy ,1)*w10+ &
ges_q_itsig (ix,iyp,1)*w01+ges_q_itsig (ixp,iyp,1)*w11)*dtsig + &
(ges_q_itsigp(ix,iy ,1)*w00+ges_q_itsigp(ixp,iy ,1)*w10+ &
ges_q_itsigp(ix,iyp,1)*w01+ges_q_itsigp(ixp,iyp,1)*w11)*dtsigp
if (sfcmod_gfs) then
tgesin =(ges_tsen(ix ,iy ,1,itsig )*w00+ &
ges_tsen(ixp,iy ,1,itsig )*w10+ &
ges_tsen(ix ,iyp,1,itsig )*w01+ &
ges_tsen(ixp,iyp,1,itsig )*w11)*dtsig + &
(ges_tsen(ix ,iy ,1,itsigp)*w00+ &
ges_tsen(ixp,iy ,1,itsigp)*w10+ &
ges_tsen(ix ,iyp,1,itsigp)*w01+ &
ges_tsen(ixp,iyp,1,itsigp)*w11)*dtsigp
prsigesin1=(ges_prsi(ix ,iy ,1,itsig )*w00+ &
ges_prsi(ixp,iy ,1,itsig )*w10+ &
ges_prsi(ix ,iyp,1,itsig )*w01+ &
ges_prsi(ixp,iyp,1,itsig )*w11)*dtsig + &
(ges_prsi(ix ,iy ,1,itsigp)*w00+ &
ges_prsi(ixp,iy ,1,itsigp)*w10+ &
ges_prsi(ix ,iyp,1,itsigp)*w01+ &
ges_prsi(ixp,iyp,1,itsigp)*w11)*dtsigp
prsigesin2=(ges_prsi(ix ,iy ,2,itsig )*w00+ &
ges_prsi(ixp,iy ,2,itsig )*w10+ &
ges_prsi(ix ,iyp,2,itsig )*w01+ &
ges_prsi(ixp,iyp,2,itsig )*w11)*dtsig + &
(ges_prsi(ix ,iy ,2,itsigp)*w00+ &
ges_prsi(ixp,iy ,2,itsigp)*w10+ &
ges_prsi(ix ,iyp,2,itsigp)*w01+ &
ges_prsi(ixp,iyp,2,itsigp)*w11)*dtsigp
call compute_fact10(ugesin,vgesin,tgesin,qgesin,pgesin, &
prsigesin1,prsigesin2,ts,sfcrough,islimsk,factw)
else if (sfcmod_mm5)then
iqtflg=.true.
lnpgesin1 =(ges_lnprsl(ix ,iy ,1,itsig )*w00+ &
ges_lnprsl(ixp,iy ,1,itsig )*w10+ &
ges_lnprsl(ix ,iyp,1,itsig )*w01+ &
ges_lnprsl(ixp,iyp,1,itsig )*w11)*dtsig + &
(ges_lnprsl(ix ,iy ,1,itsigp)*w00+ &
ges_lnprsl(ixp,iy ,1,itsigp)*w10+ &
ges_lnprsl(ix ,iyp,1,itsigp)*w01+ &
ges_lnprsl(ixp,iyp,1,itsigp)*w11)*dtsigp
lnpgesin2 =(ges_lnprsl(ix ,iy ,2,itsig )*w00+ &
ges_lnprsl(ixp,iy ,2,itsig )*w10+ &
ges_lnprsl(ix ,iyp,2,itsig )*w01+ &
ges_lnprsl(ixp,iyp,2,itsig )*w11)*dtsig + &
(ges_lnprsl(ix ,iy ,2,itsigp)*w00+ &
ges_lnprsl(ixp,iy ,2,itsigp)*w10+ &
ges_lnprsl(ix ,iyp,2,itsigp)*w01+ &
ges_lnprsl(ixp,iyp,2,itsigp)*w11)*dtsigp
tgesin =(ges_tv_itsig (ix ,iy ,1)*w00+ &
ges_tv_itsig (ixp,iy ,1)*w10+ &
ges_tv_itsig (ix ,iyp,1)*w01+ &
ges_tv_itsig (ixp,iyp,1)*w11)*dtsig + &
(ges_tv_itsigp(ix ,iy ,1)*w00+ &
ges_tv_itsigp(ixp,iy ,1)*w10+ &
ges_tv_itsigp(ix ,iyp,1)*w01+ &
ges_tv_itsigp(ixp,iyp,1)*w11)*dtsigp
tgesin2 =(ges_tv_itsig (ix ,iy ,2)*w00+ &
ges_tv_itsig (ixp,iy ,2)*w10+ &
ges_tv_itsig (ix ,iyp,2)*w01+ &
ges_tv_itsig (ixp,iyp,2)*w11)*dtsig + &
(ges_tv_itsigp(ix ,iy ,2)*w00+ &
ges_tv_itsigp(ixp,iy ,2)*w10+ &
ges_tv_itsigp(ix ,iyp,2)*w01+ &
ges_tv_itsigp(ixp,iyp,2)*w11)*dtsigp
qgesin2 =(ges_q_itsig (ix ,iy ,2)*w00+ &
ges_q_itsig (ixp,iy ,2)*w10+ &
ges_q_itsig (ix ,iyp,2)*w01+ &
ges_q_itsig (ixp,iyp,2)*w11)*dtsig + &
(ges_q_itsigp(ix ,iy ,2)*w00+ &
ges_q_itsigp(ixp,iy ,2)*w10+ &
ges_q_itsigp(ix ,iyp,2)*w01+ &
ges_q_itsigp(ixp,iyp,2)*w11)*dtsigp
geopgesin =(geop_hgtl(ix ,iy ,1,itsig )*w00+ &
geop_hgtl(ixp,iy ,1,itsig )*w10+ &
geop_hgtl(ix ,iyp,1,itsig )*w01+ &
geop_hgtl(ixp,iyp,1,itsig )*w11)*dtsig + &
(geop_hgtl(ix ,iy ,1,itsigp)*w00+ &
geop_hgtl(ixp,iy ,1,itsigp)*w10+ &
geop_hgtl(ix ,iyp,1,itsigp)*w01+ &
geop_hgtl(ixp,iyp,1,itsigp)*w11)*dtsigp
call SFC_WTQ_FWD (pgesin,ts,lnpgesin1,tgesin,qgesin,ugesin,vgesin, &
lnpgesin2,tgesin2,qgesin2,geopgesin,sfcrough,islimsk, &
factw,u10ges,v10ges,t2ges,q2ges,regime,iqtflg)
endif
return
end subroutine comp_fact10
!-------------------------------------------------------------------------
subroutine guess_grids_stats3d_(name,a,mype)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: guess_grids_stats3d_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! a
! mype - mpi task id
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: zero
use mpimod, only: ierror,mpi_rtype,mpi_sum,mpi_comm_world
use gridmod, only: lon1,lat1,nsig
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:,:,:), intent(in ) :: a
integer(i_kind) , intent(in ) :: mype
! local variables
integer(i_kind) :: i,j,k
real(r_kind),dimension(nsig+1):: work_a,work_a1
real(r_kind),dimension(nsig):: amz ! global mean profile of a
real(r_kind) :: rms
! start
! Calculate global means for a
! Calculate sums for a to estimate variance.
work_a = zero
do k = 1,nsig
do j = 2,lon1+1
do i = 2,lat1+1
work_a(k) = work_a(k) + a(i,j,k)
end do
end do
end do
work_a(nsig+1)=float(lon1*lat1)
call mpi_allreduce(work_a,work_a1,nsig+1,mpi_rtype,mpi_sum,&
mpi_comm_world,ierror)
amz=zero
do k=1,nsig
if (work_a1(nsig+1)>zero) amz(k)=work_a1(k)/work_a1(nsig+1)
rms=sqrt(amz(k)**2/work_a1(nsig+1))
if (mype==0) write(*,100) trim(name),k,amz(k),rms
enddo
100 format(a,': Level, Global mean, RMS = ',i3,1P2E16.8)
end subroutine guess_grids_stats3d_
!-------------------------------------------------------------------------
subroutine guess_grids_stats2d_(name,a,mype)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: guess_grids_stats2d_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! a
! mype - mpi task id
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: zero
use mpimod, only: ierror,mpi_rtype,mpi_sum,mpi_comm_world
use gridmod, only: lon1,lat1
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:,:), intent(in ) :: a
integer(i_kind) , intent(in ) :: mype
! local variables
integer(i_kind) :: i,j
real(r_kind),dimension(2):: work_a,work_a1
real(r_kind) :: amz, rms
! start
! Calculate global means for a
! Calculate sums for a to estimate variance.
work_a = zero
do j = 2,lon1+1
do i = 2,lat1+1
work_a(1) = work_a(1) + a(i,j)
end do
end do
work_a(2)=float(lon1*lat1)
call mpi_allreduce(work_a,work_a1,2,mpi_rtype,mpi_sum,&
mpi_comm_world,ierror)
amz=zero
if (work_a1(2)>zero) amz=work_a1(1)/work_a1(2)
rms=sqrt(amz**2/work_a1(2))
if (mype==0) write(*,100) trim(name),amz,rms
100 format(a,': Global mean, RMS = ',1P2E16.8)
end subroutine guess_grids_stats2d_
!-------------------------------------------------------------------------
subroutine pstats_(a,amiss,avg,rms)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: pstats_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! amiss
! a
!
! output argument log:
! avg,rms
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: zero
implicit none
real(r_kind),dimension(:,:), intent(in ) :: a ! array var
real(r_kind) , intent(in ) :: amiss ! undef
real(r_kind) , intent( out) :: avg,rms
! local variables
integer(i_kind) :: i,j
integer(i_kind) :: allcnt,cnt
! start
allcnt=0
cnt=0
avg=zero
rms=zero
do i=1,size(a,1)
do j=1,size(a,2)
if(a(i,j)/=amiss) then
cnt=cnt+1
avg=avg+a(i,j)
endif
allcnt = allcnt+1
end do
end do
avg=avg/max(1,cnt)
rms=sqrt(avg*avg/max(1,cnt))
end subroutine pstats_
!-------------------------------------------------------------------------
subroutine print1r8_(name,fld,undef)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: print1r8_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! fld
! undef
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:), intent(in ) :: fld
real(r_kind) , intent(in ) :: undef
!
write(6,100) trim(name),minval(fld),maxval(fld),sum(fld),undef
100 format(a,': range,sum = ',1P3E16.4)
end subroutine print1r8_
!-------------------------------------------------------------------------
subroutine print2r8_(name,fld,undef)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: print2r8_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! fld
! undef
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:,:), intent(in ) :: fld
real(r_kind) , intent(in ) :: undef
!
real(r_kind) avg,rms
write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
call pstats_(fld,UNDEF,avg,rms)
write(6,99) trim(name),avg,rms
100 format(a,': range,sum = ',1P3E16.4)
99 format(a,': avg, rms = ',1P2E16.4)
end subroutine print2r8_
!-------------------------------------------------------------------------
subroutine print3r8_(name,fld,undef,allk)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: print3r8_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! fld
! undef
! allk
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:,:,:), intent(in ) :: fld
real(r_kind) , intent(in ) :: undef
logical,optional , intent(in ) :: allk
!
logical prntlevs
integer(i_kind) k
real(r_kind) avg,rms
if(present(allk)) prntlevs=allk
if(prntlevs) then
do k=1,size(fld,3)
write(6,101) trim(name),k,minval(fld(:,:,k)),maxval(fld(:,:,k)),sum(fld(:,:,k))
call pstats_(fld(:,:,k),UNDEF,avg,rms)
write(6,99) trim(name),avg,rms
end do
else
write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
end if
101 format(a,': time or lev,range,sum = ',i3,1P3E16.4)
100 format(a,': range,sum = ',1P3E16.4)
99 format(a,': avg, rms = ',1P2E16.4)
end subroutine print3r8_
!-------------------------------------------------------------------------
subroutine print4r8_(name,fld,undef,allk)
!-------------------------------------------------------------------------
!$$$ subprogram documentation block
! . . . .
! subprogram: print4r8_
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! name
! fld
! undef
! allk
!
! output argument log:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
character(len=*) , intent(in ) :: name
real(r_kind),dimension(:,:,:,:), intent(in ) :: fld
real(r_kind) , intent(in ) :: undef
logical,optional , intent(in ) :: allk
!
logical prntlevs
integer(i_kind) k,it
real(r_kind) avg,rms
if(present(allk)) prntlevs=allk
if(prntlevs) then
do it=1,size(fld,4)
do k=1,size(fld,3)
write(6,101) trim(name),it,k,minval(fld(:,:,k,it)),maxval(fld(:,:,k,it)),sum(fld(:,:,k,it))
call pstats_(fld(:,:,k,it),UNDEF,avg,rms)
write(6,99) trim(name),avg,rms
end do
end do
else
write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
end if
101 format(a,': time,lev,range,sum = ',i3,i3,1P3E16.4)
100 format(a,': range,sum = ',1P3E16.4)
99 format(a,': avg, rms = ',1P2E16.4)
end subroutine print4r8_
end module guess_grids