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.
!
! 2019-01- Ting -- modified for fv3sar
! attributes:
! language: f95
!
!$$$
!=========================================================================
! Define associated modules
use gridinfo, only: npts
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
use params, only: nx_res,ny_res,nlevs,ntiles
use params, only: pseudo_rh
use mpeu_util, only: getindex
use read_fv3regional_restarts,only:read_fv3_restart_data1d,read_fv3_restart_data2d
use read_fv3regional_restarts,only:read_fv3_restart_data3d,read_fv3_restart_data4d
use netcdf_mod,only: nc_check
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:zero,one,half,fv, max_varname_length
use gridinfo,only: eta1_ll
use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr
use netcdf, only: nf90_inq_dimid,nf90_inq_varid
use netcdf, only: nf90_nowrite,nf90_write,nf90_inquire,nf90_inquire_dimension
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,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=:),allocatable :: fv3filename
character(len=7) :: charnanal
integer(i_kind) file_id
real(r_single), dimension(:,:,:), allocatable ::workvar3d,uworkvar3d,&
vworkvar3d,tvworkvar3d,tsenworkvar3d,&
workprsi,qworkvar3d
real(r_double),dimension(:,:,:),allocatable:: qsatworkvar3d
real(r_single), dimension(:,:), allocatable ::pswork
! Define variables required for netcdf variable I/O
character(len=12) :: varstrname
character(len=1) char_tile
character(len=24),parameter :: myname_ = 'fv3: getgriddata'
! Define counting variables
integer :: nlevsp1
integer :: i,j, k,nn,ntile,nn_tile0, nb,nanal,ne
integer :: u_ind, v_ind, tv_ind,tsen_ind, q_ind, oz_ind
integer :: ps_ind, sst_ind
integer :: tmp_ind
logical :: ice
!======================================================================
write (6,*)"The input fileprefix, reducedgrid are not used in the current implementation", &
fileprefixes, reducedgrid
nlevsp1=nlevs+1
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
tv_ind = getindex(vars3d, 't') ! 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
ps_ind = getindex(vars2d, 'ps') ! Ps (2D)
sst_ind = getindex(vars2d, 'sst') ! SST (2D)
! Initialize all constants required by routine
allocate(workvar3d(nx_res,ny_res,nlevs))
allocate(qworkvar3d(nx_res,ny_res,nlevs))
allocate(qsatworkvar3d(nx_res,ny_res,nlevs))
allocate(tvworkvar3d(nx_res,ny_res,nlevs))
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
do ntile=1,ntiles
nn_tile0=(ntile-1)*nx_res*ny_res
write(char_tile, '(i1)') ntile
filename = "fv3sar_tile"//char_tile//"_"//trim(charnanal)
fv3filename=trim(adjustl(filename))//"_dynvartracer"
!----------------------------------------------------------------------
! read u-component
call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_nowrite,file_id),&
myname_,'open: '//trim(adjustl(fv3filename)) )
!----------------------------------------------------------------------
! Update u and v variables (same for NMM and ARW)
if (u_ind > 0) then
allocate(uworkvar3d(nx_res,ny_res+1,nlevs))
varstrname = 'u'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,uworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(u_ind-1)+k,nb,ne)=uworkvar3d(i,j,k)
enddo
enddo
enddo
do k = levels(u_ind-1)+1, levels(u_ind)
if (nproc .eq. 0) &
write(6,*) 'READFVregional : u ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
deallocate(uworkvar3d)
endif
if (v_ind > 0) then
allocate(vworkvar3d(nx_res+1,ny_res,nlevs))
varstrname = 'v'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,vworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(v_ind-1)+k,nb,ne)=vworkvar3d(i,j,k)
enddo
enddo
enddo
do k = levels(v_ind-1)+1, levels(v_ind)
if (nproc .eq. 0) &
write(6,*) 'READFVregional : v ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
deallocate(vworkvar3d)
endif
if (tv_ind > 0.or.tsen_ind) then
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
varstrname = 'T'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
if (q_ind > 0) then
varstrname = 'sphum'
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(q_ind-1)+k,nb,ne)=qworkvar3d(i,j,k)
enddo
enddo
enddo
do k = levels(q_ind-1)+1, levels(q_ind)
if (nproc .eq. 0) &
write(6,*) 'READFVregional : q ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
if(tv_ind > 0) then
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
workvar3d(i,j,k)=tsenworkvar3d(i,j,k)*(one+fv*qworkvar3d(i,j,k))
enddo
enddo
enddo
tvworkvar3d=workvar3d
else! tsen_id >0
workvar3d=tsenworkvar3d
endif
tmp_ind=max(tv_ind,tsen_ind) !then can't be both >0
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(tmp_ind-1)+k,nb,ne)=workvar3d(i,j,k)
enddo
enddo
enddo
do k = levels(tmp_ind-1)+1, levels(tmp_ind)
if (nproc .eq. 0) then
write(6,*) 'READFVregional : t ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
endif
enddo
endif
if(allocated(tsenworkvar3d)) deallocate(tsenworkvar3d)
if (oz_ind > 0) then
varstrname = 'o3mr'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(oz_ind-1)+k,nb,ne)=workvar3d(i,j,k)
enddo
enddo
enddo
do k = levels(oz_ind-1)+1, levels(oz_ind)
if (nproc .eq. 0) &
write(6,*) 'READFVregional : oz ', &
& k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne))
enddo
endif
call nc_check( nf90_close(file_id),&
myname_,'close '//trim(fv3filename) )
! set SST to zero for now
if (sst_ind > 0) then
vargrid(:,levels(n3d)+sst_ind,nb,ne) = zero
endif
!----------------------------------------------------------------------
! Allocate memory for variables computed within routine
if (ps_ind > 0) then
allocate(workprsi(nx_res,ny_res,nlevsp1))
allocate(pswork(nx_res,ny_res))
fv3filename=trim(adjustl(filename))//"_dynvartracer"
call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_nowrite,file_id),&
myname_,'open: '//trim(adjustl(fv3filename)) )
call read_fv3_restart_data3d('delp',fv3filename,file_id,workvar3d)
!print *,'min/max delp',ntile,minval(delp),maxval(delp)
call nc_check( nf90_close(file_id),&
myname_,'close '//trim(fv3filename) )
workprsi(:,:,nlevsp1)=eta1_ll(nlevsp1) !etal_ll is needed
do i=nlevs,1,-1
workprsi(:,:,i)=workvar3d(:,:,i)*0.01_r_kind+workprsi(:,:,i+1)
enddo
pswork(:,:)=workprsi(:,:,1)
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
vargrid(nn,levels(n3d)+ps_ind, nb,ne) =pswork(i,j)
enddo
enddo
do k=1,nlevs
do j=1,ny_res
do i=1,nx_res
workvar3d(i,j,k)=(workprsi(i,j,k)+workprsi(i,j,k+1))*half
enddo
enddo
enddo
ice=.true. !tothink
if (pseudo_rh) then
call genqsat1(qworkvar3d,qsatworkvar3d,workvar3d,tvworkvar3d,ice, &
nx_res*ny_res,nlevs)
else
qsatworkvar3d(:,:,:) = 1._r_double
endif
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
qsat(nn,k,nb,ne)=qsatworkvar3d(i,j,k)
enddo
enddo
enddo
if(allocated(workprsi)) deallocate(workprsi)
if(allocated(pswork)) deallocate(pswork)
if(allocated(tvworkvar3d)) deallocate(tvworkvar3d)
if(allocated(qworkvar3d)) deallocate(qworkvar3d)
if(allocated(qsatworkvar3d)) deallocate(qsatworkvar3d)
endif
!======================================================================
! Deallocate memory
if(allocated(workvar3d)) deallocate(workvar3d)
end do ! ntile loop
end do backgroundloop ! loop over backgrounds to read in
end do ensmemloop ! loop over ens members to read in
return
end subroutine readgriddata
!========================================================================
! readgriddata_nmm.f90: read WRF-NMM state or control vector
!-------------------------------------------------------------------------
!========================================================================
! 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, only: zero, one,fv,half
use gridinfo,only: eta1_ll,eta2_ll
use params, only: nbackgrounds, anlfileprefixes, fgfileprefixes
use params, only: nx_res,ny_res,nlevs,ntiles,l_pres_add_saved
use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr
use netcdf, only: nf90_inq_dimid,nf90_inq_varid
use netcdf, only: nf90_write,nf90_write,nf90_inquire,nf90_inquire_dimension
use write_fv3regional_restarts,only:write_fv3_restart_data1d,write_fv3_restart_data2d
use write_fv3regional_restarts,only:write_fv3_restart_data3d,write_fv3_restart_data4d
include 'netcdf.inc'
!----------------------------------------------------------------------
! Define variables passed to subroutine
integer, intent(in) :: nanal1,nanal2, n2d, n3d, ndim
character(len=*), dimension(n2d), intent(in) :: vars2d
character(len=*), 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
!----------------------------------------------------------------------
! Define variables computed within subroutine
character(len=500) :: filename
character(len=:),allocatable :: fv3filename
character(len=7) :: charnanal
!----------------------------------------------------------------------
integer(i_kind) :: u_ind, v_ind, tv_ind, tsen_ind,q_ind, ps_ind,oz_ind
integer(i_kind) :: w_ind, cw_ind, ph_ind
integer(i_kind) file_id
real(r_single), dimension(:,:), allocatable ::pswork
real(r_single), dimension(:,:,:), allocatable ::workvar3d,workinc3d,workinc3d2,uworkvar3d,&
vworkvar3d,tvworkvar3d,tsenworkvar3d,&
workprsi,qworkvar3d
!----------------------------------------------------------------------
! Define variables required by for extracting netcdf variable
! fields
integer :: nlevsp1
! Define variables required for netcdf variable I/O
character(len=12) :: varstrname
character(len=1) char_tile
character(len=24),parameter :: myname_ = 'fv3: writegriddata'
!----------------------------------------------------------------------
! Define counting variables
integer :: i,j,k,nn,ntile,nn_tile0, nb,ne,nanal
write(6,*)"anlfileprefixes, fgfileprefixes are not used in the current implementation", &
anlfileprefixes, fgfileprefixes
write(6,*)"the no_inflate_flag is not used in the currrent implementation ",no_inflate_flag
!----------------------------------------------------------------------
nlevsp1=nlevs+1
u_ind = getindex(vars3d, 'u') !< indices in the state var arrays
v_ind = getindex(vars3d, 'v') ! U and V (3D)
tv_ind = getindex(vars3d, 't') ! Tv (3D)
tsen_ind = getindex(vars3d, 'tsen') ! 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
ps_ind = getindex(vars2d, 'ps') ! Ps (2D)
!----------------------------------------------------------------------
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
allocate(workinc3d(nx_res,ny_res,nlevs),workinc3d2(nx_res,ny_res,nlevsp1))
allocate(workvar3d(nx_res,ny_res,nlevs))
allocate(qworkvar3d(nx_res,ny_res,nlevs))
allocate(tvworkvar3d(nx_res,ny_res,nlevs))
!----------------------------------------------------------------------
! First guess file should be copied to analysis file at scripting
! level; only variables updated by EnKF are changed
write(charnanal,'(a3, i3.3)') 'mem', nanal
!----------------------------------------------------------------------
! Update u and v variables (same for NMM and ARW)
do ntile=1,ntiles
nn_tile0=(ntile-1)*nx_res*ny_res
write(char_tile, '(i1)') ntile
filename = "fv3sar_tile"//char_tile//"_"//trim(charnanal)
fv3filename=trim(adjustl(filename))//"_dynvartracer"
!----------------------------------------------------------------------
! read u-component
call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_write,file_id),&
myname_,'open: '//trim(adjustl(fv3filename)) )
! update CWM for WRF-NMM
if (u_ind > 0) then
varstrname = 'u'
allocate(uworkvar3d(nx_res,ny_res+1,nlevs))
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,uworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(u_ind-1)+k,nb,ne)
enddo
enddo
enddo
uworkvar3d(:,1:ny_res,:)=uworkvar3d(:,1:ny_res,:)+workinc3d
uworkvar3d(:,ny_res+1,:)=uworkvar3d(:,ny_res,:)
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,uworkvar3d)
deallocate(uworkvar3d)
endif
if (v_ind > 0) then
varstrname = 'v'
allocate(vworkvar3d(nx_res+1,ny_res,nlevs))
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,vworkvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(v_ind-1)+k,nb,ne)
enddo
enddo
enddo
vworkvar3d(1:nx_res,:,:)=vworkvar3d(1:nx_res,:,:)+workinc3d
vworkvar3d(nx_res+1,:,:)=vworkvar3d(nx_res,:,:)
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,vworkvar3d)
deallocate(vworkvar3d)
endif
if (tv_ind > 0.or.tsen_ind>0 ) then
varstrname = 'T'
if(tsen_ind>0) then
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tsen_ind-1)+k,nb,ne)
enddo
enddo
enddo
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
workvar3d=workvar3d+workinc3d
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
else ! tv_ind >0
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(tv_ind-1)+k,nb,ne)
enddo
enddo
enddo
varstrname = 'T'
allocate(tsenworkvar3d(nx_res,ny_res,nlevs))
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
varstrname = 'sphum'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
tvworkvar3d=tsenworkvar3d*(one+fv*qworkvar3d)
tvworkvar3d=tvworkvar3d+workinc3d
if(q_ind > 0) then
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(q_ind-1)+k,nb,ne)
enddo
enddo
enddo
qworkvar3d=qworkvar3d+workinc3d
endif
tsenworkvar3d=tvworkvar3d/(one+fv*qworkvar3d)
varstrname = 'T'
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d)
do k=1,nlevs
if (nproc .eq. 0) &
write(6,*) 'WRITEregional : T ', &
& k, minval(tsenworkvar3d(:,:,k)), maxval(tsenworkvar3d(:,:,k))
enddo
if(q_ind>0) then
varstrname='sphum'
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d)
do k=1,nlevs
if (nproc .eq. 0) &
write(6,*) 'WRITEregional : sphum ', &
& k, minval(qworkvar3d(:,:,k)), maxval(qworkvar3d(:,:,k))
enddo
endif
deallocate(tsenworkvar3d)
endif
endif
if (oz_ind > 0) then
varstrname = 'o3mr'
call read_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
do k=1,nlevs
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
workinc3d(i,j,k)=vargrid(nn,levels(oz_ind-1)+k,nb,ne)
enddo
enddo
enddo
workvar3d=workvar3d+workinc3d
call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d)
endif
if (ps_ind > 0) then
allocate(workprsi(nx_res,ny_res,nlevsp1))
allocate(pswork(nx_res,ny_res))
varstrname = 'delp'
call read_fv3_restart_data3d(varstrname,filename,file_id,workvar3d)
!print *,'min/max delp',ntile,minval(delp),maxval(delp)
workprsi(:,:,nlevsp1)=eta1_ll(nlevsp1) !etal_ll is needed
do i=nlevs,1,-1
workprsi(:,:,i)=workvar3d(:,:,i)*0.01_r_kind+workprsi(:,:,i+1)
enddo
nn = nn_tile0
do j=1,ny_res
do i=1,nx_res
nn=nn+1
pswork(i,j)=vargrid(nn,levels(n3d)+ps_ind,nb,ne)
enddo
enddo
if(l_pres_add_saved) then
do k=1,nlevs+1
do j=1,ny_res
do i=1,nx_res
workinc3d2(i,j,k)=eta2_ll(k)*pswork(i,j)
enddo
enddo
enddo
workprsi=workprsi+workinc3d2
else
workprsi(:,:,1)=workprsi(:,:,1)+pswork
do k=2,nlevsp1
workprsi(:,:,k)=eta1_ll(k)+eta2_ll(k)*workprsi(:,:,1)
enddo
endif
do k=1,nlevs
workvar3d(:,:,k)=(workprsi(:,:,k)-workprsi(:,:,k+1))*100.0
enddo
call write_fv3_restart_data3d(varstrname,filename,file_id,workvar3d)
endif
call nc_check( nf90_close(file_id),&
myname_,'close '//trim(filename) )
!----------------------------------------------------------------------
! update time stamp is to be considered NSTART_HOUR in NMM (HWRF) restart file.
!======================================================================
end do ! tiles
if(allocated(workinc3d)) deallocate(workinc3d)
if(allocated(workinc3d2)) deallocate(workinc3d2)
if(allocated(workprsi)) deallocate(workprsi)
if(allocated(pswork)) deallocate(pswork)
if(allocated(tvworkvar3d)) deallocate(tvworkvar3d)
if(allocated(qworkvar3d)) deallocate(qworkvar3d)
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
end module gridio