module gridio
!========================================================================
!$$$ Module documentation block
!
! This module contains various routines to ingest and update
! variables from Weather Research and Forecasting (WRF) model Advanced
! Research WRF (ARW) and Non-hydrostatic Mesoscale Model (NMM) dynamical
! cores which are required by the Ensemble Kalman Filter (ENKF) currently
! designed for operations within the National Centers for Environmental
! Prediction (NCEP) Global Forecasting System (GFS)
!
! prgmmr: Winterbottom org: ESRL/PSD1 date: 2011-11-30
!
! program history log:
!
! 2011-11-30 Winterbottom - Initial version.
! 2016-02-09 shlyaeva - update to read state and control variables;
! arw control now has Tv, specific humidity and
! surface pressure instead of Tp, mix ratio and
! dry surf pressure. nmm control now has Tv
! instead of Tsens.
! 2017-05-12 Y. Wang and X. Wang - add more state variables for radar DA,
! (Johnson et al. 2015 MWR; Wang and Wang
! 2017 MWR) POC: xuguang.wang@ou.edu
!
! attributes:
! language: f95
!
!$$$
!=========================================================================
! Define associated modules
use gridinfo, only: dimensions, npts, cross2dot, dot2cross
use kinds, only: r_double, r_kind, r_single, i_kind
use mpisetup, only: nproc
use netcdf_io
use params, only: nlevs, cliptracers, datapath, arw, nmm, datestring, &
pseudo_rh, nmm_restart
use mpeu_util, only: getindex
implicit none
!-------------------------------------------------------------------------
! Define all public subroutines within this module
private
public :: readgriddata
public :: writegriddata
!-------------------------------------------------------------------------
contains
! Generic WRF read routine, calls ARW-WRF or NMM-WRF
subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,reducedgrid,vargrid,qsat)
use constants, only: max_varname_length
implicit none
integer, intent(in) :: nanal1,nanal2, n2d, n3d, ndim, ntimes
character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d
character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d
integer, dimension(0:n3d), intent(in) :: levels
character(len=120), dimension(7), intent(in) :: fileprefixes
logical, intent(in) :: reducedgrid
real(r_single), dimension(npts,ndim,ntimes), intent(out) :: vargrid
real(r_double), dimension(npts,nlevs,ntimes), intent(out) :: qsat
if (arw) then
call readgriddata_arw(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,vargrid,qsat)
else if (nmm) then
call readgriddata_nmm(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,vargrid,qsat)
endif
end subroutine readgriddata
!========================================================================
! readgriddata_arw.f90: read WRF-ARW state or control vector
!-------------------------------------------------------------------------
subroutine readgriddata_arw(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,vargrid,qsat)
use constants
!======================================================================
! Define variables passed to subroutine
integer, intent(in) :: nanal1, nanal2, n2d, n3d,ndim, ntimes
character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d
character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d
integer, dimension(0:n3d), intent(in) :: levels
character(len=120), dimension(7), intent(in) :: fileprefixes
! Define variables returned by subroutine
real(r_single), dimension(npts,ndim,ntimes,nanal2-nanal1+1), intent(out) :: vargrid
real(r_double), dimension(npts,nlevs,ntimes,nanal2-nanal1+1), intent(out) :: qsat
! Define local variables
character(len=500) :: filename
character(len=7) :: charnanal
logical :: ice
real(r_single), dimension(:), allocatable :: znu, znw ! aeta1 and eta1
real(r_single), dimension(:), allocatable :: enkf_mu, enkf_mub
real(r_single), dimension(:,:), allocatable :: enkf_temp, enkf_virttemp
real(r_single), dimension(:,:), allocatable :: enkf_pressure
real(r_single), dimension(:), allocatable :: enkf_psfc
real(r_single), dimension(:,:), allocatable :: enkf_mixratio, enkf_spechumd
real(r_single), dimension(:), allocatable :: enkf_qintegral
real(r_single) :: ptop
! Define variables required for netcdf variable I/O
character(len=12) :: varstrname
! Define counting variables
integer :: i, k, nb, ne, nanal
integer :: u_ind, v_ind, tv_ind, q_ind, oz_ind, ql_ind, qr_ind, qi_ind, qg_ind, &
qs_ind, qnc_ind, qnr_ind, qni_ind, dbz_ind, w_ind
integer :: tsen_ind, prse_ind
integer :: ps_ind, sst_ind
!======================================================================
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
tv_ind = getindex(vars3d, 'tv') ! Tv (3D)
q_ind = getindex(vars3d, 'q') ! Q (3D)
oz_ind = getindex(vars3d, 'oz') ! Oz (3D)
tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D)
prse_ind = getindex(vars3d, 'prse') ! pressure
ql_ind = getindex(vars3d, 'ql') ! QL (3D)
qr_ind = getindex(vars3d, 'qr') ! QR (3D)
qi_ind = getindex(vars3d, 'qi') ! QI (3D)
qg_ind = getindex(vars3d, 'qg') ! QG (3D)
qs_ind = getindex(vars3d, 'qs') ! QS (3D)
qnc_ind = getindex(vars3d, 'qnc') ! QNC (3D)
qnr_ind = getindex(vars3d, 'qnr') ! QNR (3D)
qni_ind = getindex(vars3d, 'qni') ! QNI (3D)
dbz_ind = getindex(vars3d, 'dbz') ! DBZ (3D)
w_ind = getindex(vars3d, 'w') ! W (3D)
ps_ind = getindex(vars2d, 'ps') ! Ps (2D)
sst_ind = getindex(vars2d, 'sst') ! SST (2D)
! Initialize all constants required by routine
call init_constants(.true.)
if (ntimes > 1) then
write(6,*)'gridio/readgriddata: reading multiple backgrounds not yet supported'
call stop2(23)
endif
ne = 0
ensmemloop: do nanal=nanal1,nanal2
ne = ne + 1
backgroundloop: do nb=1,ntimes
! Define character string for ensemble member file
if (nanal > 0) then
write(charnanal,'(a3, i3.3)') 'mem', nanal
else
charnanal = 'ensmean'
endif
filename = trim(adjustl(datapath))//trim(adjustl(fileprefixes(nb)))//trim(charnanal)
!----------------------------------------------------------------------
! read u-component
if (u_ind > 0) then
varstrname = 'U'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(u_ind-1)+1:levels(u_ind),nb,ne),nlevs)
do k = levels(u_ind-1)+1, levels(u_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: u ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read v-component
if (v_ind > 0) then
varstrname = 'V'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(v_ind-1)+1:levels(v_ind),nb,ne),nlevs)
do k = levels(v_ind-1)+1, levels(v_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: v ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qcloud
if ( ql_ind > 0 ) then
varstrname = 'QCLOUD'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(ql_ind-1)+1:levels(ql_ind),nb,ne),nlevs)
do k = levels(ql_ind-1)+1, levels(ql_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: ql ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qrain
if ( qr_ind > 0 ) then
varstrname = 'QRAIN'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qr_ind-1)+1:levels(qr_ind),nb,ne),nlevs)
do k = levels(qr_ind-1)+1, levels(qr_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qr ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qice
if ( qi_ind > 0 ) then
varstrname = 'QICE'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qi_ind-1)+1:levels(qi_ind),nb,ne),nlevs)
do k = levels(qi_ind-1)+1, levels(qi_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qi ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qsnow
if ( qs_ind > 0 ) then
varstrname = 'QSNOW'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qs_ind-1)+1:levels(qs_ind),nb,ne),nlevs)
do k = levels(qs_ind-1)+1, levels(qs_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qs ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qgraup
if ( qg_ind > 0 ) then
varstrname = 'QGRAUP'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qg_ind-1)+1:levels(qg_ind),nb,ne),nlevs)
do k = levels(qg_ind-1)+1, levels(qg_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qg ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qncloud
if ( qnc_ind > 0 ) then
varstrname = 'QNCLOUD'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qnc_ind-1)+1:levels(qnc_ind),nb,ne),nlevs)
do k = levels(qnc_ind-1)+1, levels(qnc_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qnc ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qnice
if ( qi_ind > 0 ) then
varstrname = 'QNICE'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qni_ind-1)+1:levels(qni_ind),nb,ne),nlevs)
do k = levels(qni_ind-1)+1, levels(qni_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qni ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read qnrain
if ( qnr_ind > 0 ) then
varstrname = 'QNRAIN'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(qnr_ind-1)+1:levels(qnr_ind),nb,ne),nlevs)
do k = levels(qnr_ind-1)+1, levels(qnr_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: qnr ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read reflectivity
if ( dbz_ind > 0 ) then
varstrname = 'REFL_10CM'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(dbz_ind-1)+1:levels(dbz_ind),nb,ne),nlevs)
do k = levels(dbz_ind-1)+1, levels(dbz_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: dbz ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read w
if ( w_ind > 0 ) then
varstrname = 'W'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(w_ind-1)+1:levels(w_ind),nb,ne),nlevs)
do k = levels(w_ind-1)+1, levels(w_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: w ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! set ozone to zero for now (like in GSI?)
if (oz_ind > 0) then
vargrid(:,levels(oz_ind-1)+1:levels(oz_ind),nb,ne) = zero
endif
! set SST to zero for now
if (sst_ind > 0) then
vargrid(:,levels(n3d)+sst_ind,nb,ne) = zero
endif
ice = .false.
!----------------------------------------------------------------------
! Allocate memory for variables computed within routine
if(.not. allocated(enkf_temp)) allocate(enkf_temp(npts,nlevs))
if(.not. allocated(enkf_virttemp)) allocate(enkf_virttemp(npts,nlevs))
if(.not. allocated(enkf_psfc)) allocate(enkf_psfc(npts))
if(.not. allocated(enkf_qintegral)) allocate(enkf_qintegral(npts))
if(.not. allocated(enkf_pressure)) allocate(enkf_pressure(npts,nlevs))
if(.not. allocated(enkf_mixratio)) allocate(enkf_mixratio(npts,nlevs))
if(.not. allocated(enkf_spechumd)) allocate(enkf_spechumd(npts,nlevs))
!----------------------------------------------------------------------
! Ingest the perturbation potential temperature from the external file
varstrname= 'T'
call readwrfvar(filename, varstrname, enkf_temp, nlevs)
! Ingest the water vapor mixing ratio from the external file
varstrname = 'QVAPOR'
call readwrfvar(filename, varstrname, enkf_mixratio, nlevs)
! read pressure information
call readpressure_arw(filename, znu, znw, enkf_mu, enkf_mub, ptop)
! compute surface pressure
enkf_qintegral = one
do i = 1, npts
do k = 1, nlevs
enkf_qintegral(i) = enkf_qintegral(i) + &
(znw(k) - znw(k+1))*enkf_mixratio(i,k)
enddo
enddo
! compute dry surface pressure
enkf_psfc = r0_01 * (enkf_mu + enkf_mub + ptop)
! compute full surface pressure
enkf_psfc = (enkf_psfc - ptop) * enkf_qintegral + ptop
! compute specific humidity
enkf_spechumd = enkf_mixratio / (one + enkf_mixratio)
! compute pressure
do k = 1, nlevs
enkf_pressure(:,k) = r0_01 * (znu(k) * (100 * enkf_psfc - ptop) + ptop)
enddo
! compute sensible temperature
enkf_temp = (enkf_temp + 300.0) * (0.001 * enkf_pressure)**rd_over_cp_mass
! compute virtual temperature
enkf_virttemp = enkf_temp * (1. + fv*enkf_spechumd)
if (tsen_ind > 0) then
vargrid(:,levels(tsen_ind-1)+1:levels(tsen_ind),nb,ne) = enkf_temp
do k = levels(tsen_ind-1)+1, levels(tsen_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: tsen ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
if (q_ind > 0) then
vargrid(:,levels(q_ind-1)+1:levels(q_ind),nb,ne) = enkf_spechumd
do k = levels(q_ind-1)+1, levels(q_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: q ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
if (tv_ind > 0) then
vargrid(:,levels(tv_ind-1)+1:levels(tv_ind),nb,ne) = enkf_virttemp
do k = levels(tv_ind-1)+1, levels(tv_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: tv ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
if (ps_ind > 0) then
vargrid(:,levels(n3d)+ps_ind,nb,ne) = enkf_psfc
k = levels(n3d) + ps_ind
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: ps ', &
& minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
endif
if (prse_ind > 0) then
vargrid(:,levels(prse_ind-1)+1:levels(prse_ind)-1, nb,ne) = enkf_pressure
do k = levels(prse_ind-1)+1, levels(prse_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_ARW: prse ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
!----------------------------------------------------------------------
! Compute the saturation specific humidity
if (pseudo_rh) then
call genqsat1(enkf_spechumd,qsat(:,:,nb,ne),enkf_pressure,enkf_virttemp,ice, &
npts,nlevs)
else
qsat(:,:,nb,ne) = 1._r_double
endif
!======================================================================
! Deallocate memory
if(allocated(enkf_mu)) deallocate(enkf_mu)
if(allocated(enkf_mub)) deallocate(enkf_mub)
if(allocated(enkf_temp)) deallocate(enkf_temp)
if(allocated(enkf_psfc)) deallocate(enkf_psfc)
if(allocated(enkf_qintegral)) deallocate(enkf_qintegral)
if(allocated(enkf_virttemp)) deallocate(enkf_virttemp)
if(allocated(enkf_pressure)) deallocate(enkf_pressure)
if(allocated(enkf_mixratio)) deallocate(enkf_mixratio)
if(allocated(enkf_spechumd)) deallocate(enkf_spechumd)
end do backgroundloop ! loop over backgrounds to read in
end do ensmemloop ! loop over ens members to read in
return
end subroutine readgriddata_arw
!========================================================================
! readgriddata_nmm.f90: read WRF-NMM state or control vector
!-------------------------------------------------------------------------
subroutine readgriddata_nmm(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,vargrid,qsat)
use constants
!======================================================================
! Define variables passed to subroutine
integer, intent(in) :: nanal1, nanal2, n2d, n3d, ndim, ntimes
character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d
character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d
integer, dimension(0:n3d), intent(in) :: levels
character(len=120), dimension(7), intent(in) :: fileprefixes
! Define variables returned by subroutine
real(r_single), dimension(npts,ndim,ntimes,nanal2-nanal1+1), intent(out) :: vargrid
real(r_double), dimension(npts,nlevs,ntimes,nanal2-nanal1+1), intent(out) :: qsat
! Define variables computed within subroutine
logical :: ice
real :: pdtop, pt
real, dimension(:), allocatable :: aeta1, aeta2
real(r_single), dimension(:), allocatable :: enkf_psfc, enkf_pd
real(r_single), dimension(:,:), allocatable :: enkf_temp
real(r_single), dimension(:,:), allocatable :: enkf_pressure
real(r_single), dimension(:,:), allocatable :: enkf_spechumd
character(len=12) :: varstrname
character(len=500) :: filename
character(len=7) :: charnanal
! Define counting variables
integer(i_kind) :: nb, k, nanal, ne
integer(i_kind) :: u_ind, v_ind, tv_ind, q_ind, oz_ind
integer(i_kind) :: cw_ind, tsen_ind, prse_ind
integer(i_kind) :: ps_ind, sst_ind
!======================================================================
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
tv_ind = getindex(vars3d, 'tv') ! Tv (3D)
q_ind = getindex(vars3d, 'q') ! Q (3D)
oz_ind = getindex(vars3d, 'oz') ! Oz (3D)
cw_ind = getindex(vars3d, 'cw') ! CW (3D)
tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D)
prse_ind = getindex(vars3d, 'prse')
ps_ind = getindex(vars2d, 'ps') ! Ps (2D)
sst_ind = getindex(vars2d, 'sst')
! Initialize all constants required by routine
call init_constants(.true.)
!----------------------------------------------------------------------
if (ntimes > 1) then
write(6,*)'gridio/readgriddata: reading multiple backgrounds not yet supported'
call stop2(23)
endif
ne = 0
ensmemloop: do nanal=nanal1,nanal2
ne = ne + 1
backgroundloop: do nb=1,ntimes
! Define character string for ensemble member file
if (nanal > 0) then
write(charnanal,'(a3, i3.3)') 'mem', nanal
else
charnanal = 'ensmean'
endif
filename = trim(adjustl(datapath))//trim(adjustl(fileprefixes(nb)))//trim(charnanal)
!----------------------------------------------------------------------
! read u-component
if (u_ind > 0) then
varstrname = 'U'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(u_ind-1)+1:levels(u_ind),nb,ne), nlevs)
do k = levels(u_ind-1)+1, levels(u_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: u ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read v-component
if (v_ind > 0) then
varstrname = 'V'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(v_ind-1)+1:levels(v_ind),nb,ne), nlevs)
do k = levels(v_ind-1)+1, levels(v_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: v ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! read cwm
if (cw_ind > 0) then
varstrname = 'CWM'
call readwrfvar(filename, varstrname, &
vargrid(:,levels(cw_ind-1)+1:levels(cw_ind),nb,ne), nlevs)
do k = levels(cw_ind-1)+1, levels(cw_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: cw', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! set ozone to zero for now (like in GSI?)
if (oz_ind > 0) then
vargrid(:,levels(oz_ind-1)+1:levels(oz_ind),nb,ne) = zero
endif
! set SST to zero for now
if (sst_ind > 0) then
vargrid(:,levels(n3d)+sst_ind,nb,ne) = zero
endif
! Define all constants required by routine
ice = .false.
!----------------------------------------------------------------------
! Allocate memory for variables computed within routine
if(.not. allocated(enkf_psfc)) allocate(enkf_psfc(npts))
if(.not. allocated(enkf_temp)) allocate(enkf_temp(npts,nlevs))
if(.not. allocated(enkf_pressure)) allocate(enkf_pressure(npts,nlevs))
if(.not. allocated(enkf_spechumd)) allocate(enkf_spechumd(npts,nlevs))
!----------------------------------------------------------------------
! Ingest the (sensible) temperature from the external file
varstrname= 'T'
call readwrfvar(filename, varstrname, enkf_temp, nlevs)
if (tsen_ind > 0) then
vargrid(:,levels(tsen_ind-1)+1:levels(tsen_ind),nb,ne) = enkf_temp
do k = levels(tsen_ind-1)+1, levels(tsen_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: tsen ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! Ingest the specific humidity from the external file
varstrname = 'Q'
call readwrfvar(filename, varstrname, enkf_spechumd, nlevs)
if (q_ind > 0) then
vargrid(:,levels(q_ind-1)+1:levels(q_ind),nb,ne) = enkf_spechumd
do k = levels(q_ind-1)+1, levels(q_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: q ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! calculate virtual temperature
enkf_temp = enkf_temp * (one + fv*enkf_spechumd)
if (tv_ind > 0) then
vargrid(:,levels(tv_ind-1)+1:levels(tv_ind),nb,ne) = enkf_temp
do k = levels(tv_ind-1)+1, levels(tv_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: tv ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
! Ingest pressure info from the external file
call readpressure_nmm(filename, enkf_pd, aeta1, aeta2, pt, pdtop)
!----------------------------------------------------------------------
! calculate surface pressure
enkf_psfc = r0_01 * (enkf_pd + pdtop + pt)
if (ps_ind > 0) then
vargrid(:,levels(n3d)+ps_ind,nb,ne) = enkf_psfc
k = levels(n3d) + ps_ind
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: ps ', &
& minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
endif
! compute the pressure profile
do k = 1, nlevs
enkf_pressure(:,k) = r0_01 * (aeta1(k)*pdtop + aeta2(k)*enkf_pd + pt)
end do
if (prse_ind > 0) then
vargrid(:,levels(prse_ind-1)+1:levels(prse_ind)-1, nb,ne) = enkf_pressure
do k = levels(prse_ind-1)+1, levels(prse_ind)
if (nproc .eq. 0) &
write(6,*) 'READGRIDDATA_NMM: prse ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
!----------------------------------------------------------------------
! Compute the saturation specific humidity
if (pseudo_rh) then
call genqsat1(enkf_spechumd,qsat(:,:,nb,ne),enkf_pressure,enkf_temp,ice, &
npts,nlevs)
else
qsat(:,:,nb,ne) = 1._r_double
endif
! Deallocate memory for variables computed within routine
if(allocated(enkf_temp)) deallocate(enkf_temp)
if(allocated(enkf_psfc)) deallocate(enkf_psfc)
if(allocated(enkf_pressure)) deallocate(enkf_pressure)
if(allocated(enkf_spechumd)) deallocate(enkf_spechumd)
!======================================================================
end do backgroundloop ! loop over backgrounds to read in
end do ensmemloop ! loop over ens members to read in
! Return calculated values
return
end subroutine readgriddata_nmm
!========================================================================
! writegriddata.f90: write WRF-ARW or WRF-NMM analysis
!-------------------------------------------------------------------------
subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflate_flag)
use constants
use params, only: nbackgrounds, anlfileprefixes, fgfileprefixes
include 'netcdf.inc'
!----------------------------------------------------------------------
! Define variables passed to subroutine
integer, intent(in) :: nanal1,nanal2, n2d, n3d, ndim
character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d
character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d
integer, dimension(0:n3d), intent(in) :: levels
real(r_single), dimension(npts,ndim,nbackgrounds,nanal2-nanal1+1), intent(in) :: vargrid
logical, intent(in) :: no_inflate_flag
!Not used here, but added to make writegriddata(...) consistent with gridio_gfs.f90
!----------------------------------------------------------------------
! Define variables computed within subroutine
character(len=500) :: filename
character(len=3) :: charnanal
real :: clip
integer :: iyear,imonth,iday,ihour,dh1,ierr,iw3jdn
!----------------------------------------------------------------------
integer(i_kind) :: u_ind, v_ind, tv_ind, q_ind, ps_ind, ql_ind, qr_ind, qi_ind, qg_ind, &
qs_ind, qnc_ind, qnr_ind, qni_ind, dbz_ind
integer(i_kind) :: w_ind, cw_ind, ph_ind
!----------------------------------------------------------------------
! Define variables required by for extracting netcdf variable
! fields
character(len=19) :: DateStr
! Define variables required for netcdf variable I/O
character(len=12) :: varstrname
real, dimension(:,:,:), allocatable :: vargrid_native
!----------------------------------------------------------------------
! Define counting variables
integer :: k, nb, nanal, ne
real(r_single), dimension(:,:), allocatable :: enkf_t, enkf_q, enkf_field
real(r_single), dimension(:), allocatable :: enkf_psfc, pressure
real(r_single), dimension(:), allocatable :: enkf_mu, enkf_mub
real(r_single), dimension(:), allocatable :: znu, znw
real(r_single), dimension(:), allocatable :: qintegral
real(r_single) :: ptop
!----------------------------------------------------------------------
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
tv_ind = getindex(vars3d, 'tv') ! Tv (3D)
q_ind = getindex(vars3d, 'q') ! Q (3D)
cw_ind = getindex(vars3d, 'cw') ! CWM for WRF-NMM
w_ind = getindex(vars3d, 'w') ! W for WRF-ARW
ph_ind = getindex(vars3d, 'ph') ! PH for WRF-ARW
ql_ind = getindex(vars3d, 'ql') ! QL (3D) for WRF-ARW
qr_ind = getindex(vars3d, 'qr') ! QR (3D) for WRF-ARW
qi_ind = getindex(vars3d, 'qi') ! QI (3D) for WRF-ARW
qg_ind = getindex(vars3d, 'qg') ! QG (3D) for WRF-ARW
qs_ind = getindex(vars3d, 'qs') ! QS (3D) for WRF-ARW
qnc_ind = getindex(vars3d, 'qnc') ! QNC (3D) for WRF-ARW
qnr_ind = getindex(vars3d, 'qnr') ! QNR (3D) for WRF-ARW
qni_ind = getindex(vars3d, 'qni') ! QNI (3D) for WRF-ARW
dbz_ind = getindex(vars3d, 'dbz') ! DBZ (3D) for WRF-ARW
ps_ind = getindex(vars2d, 'ps') ! Ps (2D)
! Initialize constants required by routine
call init_constants(.true.)
!----------------------------------------------------------------------
if (nbackgrounds > 1) then
write(6,*)'gridio/writegriddata: writing multiple backgrounds not yet supported'
call stop2(23)
endif
ne = 0
ensmemloop: do nanal=nanal1,nanal2
ne = ne + 1
backgroundloop: do nb=1,nbackgrounds
!----------------------------------------------------------------------
! First guess file should be copied to analysis file at scripting
! level; only variables updated by EnKF are changed
write(charnanal,'(i3.3)') nanal
filename = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"mem"//charnanal
!----------------------------------------------------------------------
! Update u and v variables (same for NMM and ARW)
allocate(enkf_field(npts, nlevs))
if (u_ind > 0) then
varstrname = 'U'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(u_ind-1)+1:levels(u_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (v_ind > 0) then
varstrname = 'V'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(v_ind-1)+1:levels(v_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
! update CWM for WRF-NMM
if (nmm .and. cw_ind > 0) then
varstrname = 'CWM'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(cw_ind-1)+1:levels(cw_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
! update reflectivity and hydrometeor mixing ratios for WRF-ARW
if (arw .and. dbz_ind > 0) then
varstrname = 'REFL_10CM'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(dbz_ind-1)+1:levels(dbz_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. ql_ind > 0) then
varstrname = 'QCLOUD'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(ql_ind-1)+1:levels(ql_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qr_ind > 0) then
varstrname = 'QRAIN'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qr_ind-1)+1:levels(qr_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qi_ind > 0) then
varstrname = 'QICE'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qi_ind-1)+1:levels(qi_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qs_ind > 0) then
varstrname = 'QSNOW'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qs_ind-1)+1:levels(qs_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qg_ind > 0) then
varstrname = 'QGRAUP'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qg_ind-1)+1:levels(qg_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qnc_ind > 0) then
varstrname = 'QNCLOUD'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qnc_ind-1)+1:levels(qnc_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qni_ind > 0) then
varstrname = 'QNICE'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qni_ind-1)+1:levels(qni_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. qnr_ind > 0) then
varstrname = 'QNRAIN'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(qnr_ind-1)+1:levels(qnr_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
! update W and PH for WRF-ARW
if (arw .and. w_ind > 0) then
varstrname = 'W'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(w_ind-1)+1:levels(w_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
if (arw .and. ph_ind > 0) then
varstrname = 'PH'
call readwrfvar(filename, varstrname, enkf_field, nlevs)
enkf_field = enkf_field + vargrid(:,levels(ph_ind-1)+1:levels(ph_ind),nb,ne)
call writewrfvar(filename, varstrname, enkf_field, nlevs)
endif
deallocate(enkf_field)
allocate(enkf_t(npts, nlevs), enkf_q(npts,nlevs), enkf_psfc(npts))
if (nmm) then
! Update Tv and Q for NMM files (write out Tsen and Q)
if (tv_ind > 0 .or. q_ind > 0) then
! read background specific humidity and sensible temperature
varstrname = 'Q'
call readwrfvar(filename, varstrname, enkf_q, nlevs)
varstrname = 'T'
call readwrfvar(filename, varstrname, enkf_t, nlevs)
! compute background virtual temperature
enkf_t = enkf_t * (one + fv*enkf_q)
! add analysis increment to virtual temperature and specific humidity
if (tv_ind > 0) then
enkf_t = enkf_t + vargrid(:,levels(tv_ind-1)+1:levels(tv_ind),nb,ne)
endif
if (q_ind > 0) then
enkf_q = enkf_q + vargrid(:,levels(q_ind-1)+1:levels(q_ind),nb,ne)
endif
! clip Q if needed
if (cliptracers) then
clip = tiny(enkf_q(1,1))
where (enkf_q < clip) enkf_q = clip
end if
! compute analysis sensible temperature
enkf_t = enkf_t / (one + fv*enkf_q)
! write out analysis sensible temperature and specific humidity
if (tv_ind > 0) then
varstrname = 'T'
call writewrfvar(filename, varstrname, enkf_t, nlevs)
endif
if (q_ind > 0) then
varstrname = 'Q'
call writewrfvar(filename, varstrname, enkf_q, nlevs)
endif
endif
! update surface pressure for NMM
if (ps_ind > 0) then
varstrname = 'PD'
call readwrfvar(filename, varstrname, enkf_psfc, 1)
! add ps increment (mulitply by 100 since we're updating PD
enkf_psfc = enkf_psfc + 100.*vargrid(:,levels(n3d)+ps_ind,nb,ne)
call writewrfvar(filename, varstrname, enkf_psfc, 1)
endif
! for ARW, update Tv and Q, but write out Tp and mix ratio
elseif (arw) then
if (tv_ind > 0 .or. q_ind > 0 .or. ps_ind > 0) then
allocate(qintegral(npts), pressure(npts))
! read background potential temperature, mixing ratio
! and pressure information
varstrname = 'QVAPOR'
call readwrfvar(filename, varstrname, enkf_q, nlevs)
varstrname = 'T'
call readwrfvar(filename, varstrname, enkf_t, nlevs)
call readpressure_arw(filename, znu, znw, enkf_mu, enkf_mub, ptop)
! compute background dry surface pressure
enkf_psfc = r0_01*(enkf_mu + enkf_mub + ptop)
! compute background full surface pressure
qintegral = one
do k = 1, nlevs
qintegral(:) = qintegral(:) + (znw(k) - znw(k+1))*enkf_q(:,k)
enddo
enkf_psfc = (enkf_psfc - ptop) * qintegral + ptop
! compute background specific humidity
enkf_q = enkf_q / (one + enkf_q)
! compute background sensible temperature
do k = 1, nlevs
pressure = r0_01 * (znu(k)*(100*enkf_psfc-ptop)+ptop)
enkf_t(:,k) = (enkf_t(:,k) + 300.0) * &
(0.001 * pressure)**rd_over_cp_mass
enddo
! compute background virtual temperature
enkf_t = enkf_t * (one + fv*enkf_q)
! add analysis increment to virtual temperature, specific humidity
! and surface pressure
if (tv_ind > 0) then
enkf_t = enkf_t + vargrid(:,levels(tv_ind-1)+1:levels(tv_ind),nb,ne)
endif
if (q_ind > 0) then
enkf_q = enkf_q + vargrid(:,levels(q_ind-1)+1:levels(q_ind),nb,ne)
endif
if (ps_ind > 0) then
enkf_psfc = enkf_psfc + vargrid(:,levels(n3d)+ps_ind,nb,ne)
endif
! clip Q if needed
if (cliptracers) then
clip = tiny(enkf_q(1,1))
where (enkf_q < clip) enkf_q = clip
end if
! compute analysis sensible temperature
enkf_t = enkf_t / (one + fv*enkf_q)
! compute analysis mixing ratio
enkf_q = enkf_q / (one - enkf_q)
! compute analysis potential temperature
do k = 1, nlevs
pressure = r0_01 * (znu(k)*(100*enkf_psfc-ptop)+ptop)
enkf_t(:,k) = enkf_t(:,k) / &
(0.001 * pressure)**rd_over_cp_mass - 300.0
enddo
! compute analysis dry surface pressure
qintegral = one
do k = 1, nlevs
qintegral(:) = qintegral(:) + &
(znw(k) - znw(k+1))*enkf_q(:,k)
enddo
enkf_psfc = (enkf_psfc - ptop) / qintegral + ptop
! compute analysis mu
enkf_psfc = 100.*enkf_psfc - enkf_mub - ptop
! write out analysis virtual temperature, specific humidity
! and surface pressure
if (tv_ind > 0) then
varstrname = 'T'
call writewrfvar(filename, varstrname, enkf_t, nlevs)
endif
if (q_ind > 0) then
varstrname = 'QVAPOR'
call writewrfvar(filename, varstrname, enkf_q, nlevs)
endif
if (ps_ind > 0) then
varstrname = 'MU'
call writewrfvar(filename, varstrname, enkf_psfc, 1)
endif
endif
endif
!----------------------------------------------------------------------
! update NSTART_HOUR in NMM (HWRF) restart file.
read(datestring(1:4),'(i4)') iyear
read(datestring(5:6),'(i2)') imonth
read(datestring(7:8),'(i2)') iday
read(datestring(9:10),'(i2)') ihour
if (nmm .and. nmm_restart) then
varstrname = 'NSTART_HOUR'
if(.not. allocated(vargrid_native)) allocate(vargrid_native(1,1,1))
vargrid_native(1,1,1) = ihour
call writenetcdfdata(filename,vargrid_native,varstrname,1,1,1)
end if
!
! update START_DATE, SIMULATION_START_DATE, GMT, JULYR, JULDAY
! global attributes.
!
write(DateStr,'(i4,"-",i2.2,"-",i2.2,"-",i2.2,"_",i2.2,":",i2.2)') iyear,imonth,iday,ihour,0,0
ierr = NF_OPEN(trim(filename), NF_WRITE, dh1)
IF (ierr .NE. NF_NOERR) print *, 'OPEN ',NF_STRERROR(ierr)
ierr = NF_PUT_ATT_TEXT(dh1,NF_GLOBAL,'START_DATE',len(trim(DateStr)),DateStr)
IF (ierr .NE. NF_NOERR) print *,'PUT START_DATE', NF_STRERROR(ierr)
ierr = NF_PUT_ATT_TEXT(dh1,NF_GLOBAL,'SIMULATION_START_DATE',len(trim(DateStr)),DateStr)
IF (ierr .NE. NF_NOERR) print *,'PUT SIMULATION_START_DATE', NF_STRERROR(ierr)
ierr = NF_PUT_ATT_REAL(dh1,NF_GLOBAL,'GMT',NF_FLOAT,1,float(ihour))
IF (ierr .NE. NF_NOERR) print *,'PUT GMT', NF_STRERROR(ierr)
ierr = NF_PUT_ATT_INT(dh1,NF_GLOBAL,'JULYR',NF_INT,1,iyear)
IF (ierr .NE. NF_NOERR) print *,'PUT JULYR', NF_STRERROR(ierr)
ierr=NF_PUT_ATT_INT(dh1,NF_GLOBAL,'JULDAY',NF_INT,1,iw3jdn(iyear,imonth,iday)-iw3jdn(iyear,1,1)+1)
IF (ierr .NE. NF_NOERR) print *,'PUT JULDAY', NF_STRERROR(ierr)
ierr = NF_CLOSE(dh1)
IF (ierr .NE. NF_NOERR) print *, 'CLOSE ',NF_STRERROR(ierr)
!======================================================================
end do backgroundloop ! loop over backgrounds to read in
end do ensmemloop ! loop over ens members to read in
! Return calculated values
return
!======================================================================
end subroutine writegriddata
!======================================================================
! readwrfvar.f90: This subroutine reads a varname variable from WRF
! ARW or NMM netcdf file and returns the variable interpolated to
! unstaggered grid, in EnKF style (1D array for 2D field); all
! checks for grid staggering are contained within this subroutine
subroutine readwrfvar(filename, varname, grid, nlevs)
implicit none
character(len=500), intent(in) :: filename
character(len=12), intent(in) :: varname
integer(i_kind), intent(in) :: nlevs
real(r_single), dimension(npts,nlevs), intent(out) :: grid
! Define variables computed within subroutine
real, dimension(:,:,:), allocatable :: workgrid
real, dimension(:,:,:), allocatable :: vargrid_native
integer :: xdim, ydim, zdim
integer :: xdim_native, ydim_native, zdim_native
integer :: xdim_local, ydim_local, zdim_local
! Define variables requiredfor netcdf variable I/O
character(len=50) :: attstr
character(len=12) :: varstagger
character(len=12) :: varmemoryorder
! Define counting variables
integer :: i, j, k
integer :: counth
xdim = dimensions%xdim
ydim = dimensions%ydim
zdim = dimensions%zdim
! Define staggering attributes for variable grid
attstr = 'stagger'
call variableattribute_char(filename,varname,attstr, &
& varstagger)
xdim_native = xdim
ydim_native = ydim
zdim_native = zdim
! If variable grid is staggered assign array dimensions appropriately
if(varstagger(1:1) .eq. 'X') then
xdim_native = xdim + 1
else if(varstagger(1:1) .eq. 'Y') then
ydim_native = ydim + 1
else if(varstagger(1:1) .eq. 'Z') then
zdim_native = zdim + 1
end if ! if(varstagger(1:1) .eq. 'X')
! Define memory attributes for variable grid
attstr = 'MemoryOrder'
call variableattribute_char(filename,varname,attstr, &
& varmemoryorder)
! If variable is a 2-dimensional field, rescale variables appropriately
if(varmemoryorder(1:3) .eq. 'XY ') then
zdim_local = 1
zdim_native = 1
else
zdim_local = zdim
end if ! if(varmemoryorder(1:3) .eq. 'XY ')
! Define local variable dimensions
xdim_local = xdim
ydim_local = ydim
! Allocate memory for local variable arrays
if(.not. allocated(workgrid)) &
& allocate(workgrid(xdim_local,ydim_local,zdim_local))
if(.not. allocated(vargrid_native)) &
& allocate(vargrid_native(xdim_native,ydim_native,zdim_native))
! Ingest variable from external netcdf formatted file
call readnetcdfdata(filename,vargrid_native,varname, &
& xdim_native,ydim_native,zdim_native)
! Interpolate variable from staggered (i.e., E-) grid to
! unstaggered (i.e., A-) grid. If variable is staggered in
! vertical, intepolate from model layer interfaces
! (including surface and top) to model layer midpoints.
call cross2dot(vargrid_native,xdim_native,ydim_native, &
& zdim_native,xdim_local,ydim_local,zdim_local,workgrid)
!----------------------------------------------------------------------
! Loop through vertical coordinate
do k = 1, zdim_local
! Initialize counting variable
counth = 1
! Loop through meridional horizontal coordinate
do j = 1, ydim_local
! Loop through zonal horizontal coordinate
do i = 1, xdim_local
! Assign values to output variable array
grid(counth,k) = workgrid(i,j,k)
counth = counth + 1
end do ! do i = 1, xdim_local
end do ! do j = 1, ydim_local
end do ! do k = 1, zdim_local
!----------------------------------------------------------------------
! Deallocate memory for local variables
if(allocated(vargrid_native)) deallocate(vargrid_native)
if(allocated(workgrid)) deallocate(workgrid)
end subroutine readwrfvar
!======================================================================
! writewrfvar: write EnKF-style field in WRF netcdf file; variable is
! interpolated to the native variable grid; all checks for
! grid staggering are contained within this subroutine
subroutine writewrfvar(filename, varname, grid, nlevs)
implicit none
character(len=500), intent(in) :: filename
character(len=12), intent(in) :: varname
integer(i_kind), intent(in) :: nlevs
real(r_single), dimension(npts,nlevs), intent(in) :: grid
! Define variables computed within subroutine
real, dimension(:,:,:), allocatable :: workgrid
real, dimension(:,:,:), allocatable :: vargrid_native
integer :: xdim, ydim, zdim
integer :: xdim_native, ydim_native, zdim_native
integer :: xdim_local, ydim_local, zdim_local
! Define variables requiredfor netcdf variable I/O
character(len=50) :: attstr
character(len=12) :: varstagger
! Define counting variables
integer :: i, j, k
integer :: counth
xdim = dimensions%xdim
ydim = dimensions%ydim
zdim = dimensions%zdim
! Allocate memory for local variable
allocate(workgrid(xdim,ydim,zdim))
xdim_native = xdim
ydim_native = ydim
zdim_native = zdim
if (arw) then
attstr = 'stagger'
call variableattribute_char(filename,varname,attstr, &
& varstagger)
!----------------------------------------------------------------------
! If variable grid is staggered, assign array dimensions appropriately
if(varstagger(1:1) .eq. 'X') then
xdim_native = xdim + 1
else if(varstagger(1:1) .eq. 'Y') then
ydim_native = ydim + 1
else if(varstagger(1:1) .eq. 'Z') then
zdim_native = zdim + 1
end if ! if(varstagger(1:1) .eq. 'X')
endif
zdim_local = nlevs
if(nlevs == 1) then
zdim_native = 1
end if
! Define local variable dimensions
xdim_local = xdim
ydim_local = ydim
!----------------------------------------------------------------------
! Allocate memory local arrays (first check whether they are
! already allocated)
if (allocated(vargrid_native)) deallocate(vargrid_native)
allocate(vargrid_native(xdim_native,ydim_native,zdim_native))
!----------------------------------------------------------------------
! Loop through vertical coordinate
do k = 1, zdim_local
! Initialize counting variable
counth = 1
! Loop through meridional horizontal coordinate
do j = 1, ydim
! Loop through zonal horizontal coordinate
do i = 1, xdim
! Assign values to local array
workgrid(i,j,k) = grid(counth,k)
counth = counth + 1
end do ! do i = 1, xdim
end do ! do j = 1, ydim
end do ! k = 1, zdim_local
! Interpolate increments to native grid (i.e., from A-grid to
! C-grid; if necessary); on input, workgrid is increments on
! unstaggered grid; on output vargrid_native is increments on
! model-native (i.e., staggered grid); vargridin_native is
! unmodified first guess on native staggered grid
call dot2cross(xdim_local,ydim_local,zdim_local,xdim_native, &
ydim_native,zdim_native,workgrid,vargrid_native)
!----------------------------------------------------------------------
! Write analysis variable.
call writenetcdfdata(filename,vargrid_native,varname, &
xdim_native,ydim_native,zdim_native)
! Deallocate memory for local variables
if(allocated(vargrid_native)) deallocate(vargrid_native)
if(allocated(workgrid)) deallocate(workgrid)
end subroutine writewrfvar
!========================================================================
! read pressure information (pd, aeta1, aeta2, pl, pdtop from WRF-NMM file
! subroutine allocates space for pd, aeta1 and aeta2
subroutine readpressure_nmm(filename, pd, aeta1, aeta2, pt, pdtop)
implicit none
character(len=500), intent(in) :: filename
real(r_single), dimension(:), allocatable :: aeta1, aeta2
real(r_single), dimension(:), allocatable :: pd
real(r_single) :: pt, pdtop
real, dimension(:,:,:), allocatable :: wrfnmm
character(len=12) :: varstrname
integer :: zdim
zdim = dimensions%zdim
allocate(aeta1(zdim), aeta2(zdim), pd(npts))
! Ingest surface pressure from the external file
varstrname = 'PD'
call readwrfvar(filename,varstrname,pd,1)
! Ingest hybrid vertical coordinate from the external file
varstrname = 'AETA1'
allocate(wrfnmm(1,1,zdim))
call readnetcdfdata(filename,wrfnmm,varstrname,1,1,zdim)
aeta1 = wrfnmm(1,1,:)
varstrname = 'AETA2'
call readnetcdfdata(filename,wrfnmm,varstrname,1,1,zdim)
aeta2 = wrfnmm(1,1,:)
deallocate(wrfnmm)
allocate(wrfnmm(1,1,1))
! Ingest pressure at top of domain from the external file
varstrname = 'PT'
call readnetcdfdata(filename,wrfnmm,varstrname,1,1,1)
pt = wrfnmm(1,1,1)
! Ingest mass within pressure domain from the external file
varstrname = 'PDTOP'
call readnetcdfdata(filename,wrfnmm,varstrname,1,1,1)
pdtop = wrfnmm(1,1,1)
deallocate(wrfnmm)
end subroutine readpressure_nmm
!========================================================================
! read pressure information (aeta1, eta1, mu, mub, ptop from WRF-ARW file
! subroutine allocates space for znu, znw, mu, mub
subroutine readpressure_arw(filename, znu, znw, mu, mub, ptop)
implicit none
character(len=500), intent(in) :: filename
real(r_single), dimension(:), allocatable :: znu, znw ! aeta1 and eta1
real(r_single), dimension(:), allocatable :: mu, mub
real(r_single) :: ptop
real, dimension(:,:,:), allocatable :: wrfarw
character(len=12) :: varstrname
integer :: zdim
zdim = dimensions%zdim
allocate(znu(zdim), znw(zdim + 1), mu(npts), mub(npts))
! Ingest the model vertical (eta) levels from the external file
varstrname = 'ZNU'
allocate(wrfarw(1,1,zdim))
call readnetcdfdata(filename,wrfarw,varstrname,1,1,zdim)
znu = wrfarw(1,1,:)
deallocate(wrfarw)
! Ingest the model vertical (aeta) levels from the external file
varstrname = 'ZNW'
allocate(wrfarw(1,1,zdim+1))
call readnetcdfdata(filename,wrfarw,varstrname,1,1,zdim+1)
znw = wrfarw(1,1,:)
deallocate(wrfarw)
! Ingest the model top pressure level from the external file
varstrname = 'P_TOP'
allocate(wrfarw(1,1,1))
call readnetcdfdata(filename,wrfarw,varstrname,1,1,1)
ptop = wrfarw(1,1,1)
deallocate(wrfarw)
! Ingest the model perturbation dry air mass from the external
! file
varstrname = 'MU'
call readwrfvar(filename,varstrname,mu,1)
! Ingest the model base state dry air mass from the external file
varstrname = 'MUB'
call readwrfvar(filename,varstrname,mub,1)
end subroutine readpressure_arw
end module gridio