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 ! 2021-02-08 CAPS (C. Tong and J. Park) ! -- add variables for direct reflectivitiy DA capability ! (hydrometeors and 'w') ! -- add code to update 'delp' directly ! from analysis icnrements ! attributes: ! language: f95 ! !$$$ !========================================================================= ! Define associated modules use gridinfo, only: npts use constants, only:max_varname_length 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,l_fv3reg_filecombined use params, only: pseudo_rh, l_use_enkf_directZDA 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,readgriddata_pnc public :: writegriddata,writegriddata_pnc public :: writeincrement, writeincrement_pnc !------------------------------------------------------------------------- integer(i_kind) ,parameter:: ndynvarslist=6, ntracerslist=8 character(len=max_varname_length), parameter, dimension(ndynvarslist) :: & vardynvars = [character(len=max_varname_length) :: & 'u', 'v', 'T', 'W', 'DZ', 'delp'] character(len=max_varname_length), parameter, dimension(ntracerslist) :: & vartracers = [character(len=max_varname_length) :: & 'sphum','o3mr', 'liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc'] type type_fv3lamfile logical l_filecombined character(len=max_varname_length), dimension(2):: fv3lamfilename integer (i_kind), dimension(2):: fv3lam_fileid(2) contains procedure, pass(this) :: setupfile => type_bound_setupfile procedure, pass(this):: get_idfn => type_bound_getidfn end type type(type_fv3lamfile) :: fv3lamfile contains subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,fileprefixes,filesfcprefixes,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 character(len=120), dimension(7), intent(in) :: filesfcprefixes 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,fv3filename1 character(len=7) :: charnanal integer(i_kind) file_id,file_id1 real(r_single), dimension(:,:,:), allocatable ::workvar3d,uworkvar3d,& vworkvar3d,tvworkvar3d,tsenworkvar3d,& workprsi,qworkvar3d,wworkvar3d 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 :: w_ind, ql_ind, qi_ind, qr_ind, qs_ind, qg_ind, qnr_ind integer :: ps_ind, sst_ind integer :: tmp_ind,ifile 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) w_ind = getindex(vars3d, 'w') ! W (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 ql_ind = getindex(vars3d, 'ql') ! Q cloud water (3D) qi_ind = getindex(vars3d, 'qi') ! Q cloud ice (3D) qr_ind = getindex(vars3d, 'qr') ! Q rain water (3D) qs_ind = getindex(vars3d, 'qs') ! Q snow (3D) qg_ind = getindex(vars3d, 'qg') ! Q graupel (3D) qnr_ind = getindex(vars3d, 'qnr') ! N rain (3D) 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) if(l_fv3reg_filecombined) then fv3filename=trim(adjustl(filename))//"_dynvartracer" call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_nowrite,file_id),& myname_,'open: '//trim(adjustl(fv3filename)) ) call fv3lamfile%setupfile(fileid1=file_id,fv3fn1=trim(adjustl(fv3filename))) else fv3filename=trim(adjustl(filename))//"_dynvars" call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_nowrite,file_id),& myname_,'open: '//trim(adjustl(fv3filename)) ) fv3filename1=trim(adjustl(filename))//"_tracer" call nc_check( nf90_open(trim(adjustl(fv3filename1)),nf90_nowrite,file_id1),& myname_,'open: '//trim(adjustl(fv3filename1)) ) call fv3lamfile%setupfile(fileid1=file_id,fv3fn1=trim(adjustl(fv3filename)) , & fileid2=file_id1,fv3fn2=trim(adjustl(fv3filename1)) ) endif !---------------------------------------------------------------------- ! read u-component !---------------------------------------------------------------------- ! 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 (w_ind > 0) then allocate(wworkvar3d(nx_res,ny_res,nlevs)) varstrname = 'W' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,wworkvar3d) do k=1,nlevs nn = nn_tile0 do j=1,ny_res do i=1,nx_res nn=nn+1 vargrid(nn,levels(w_ind-1)+k,nb,ne)=wworkvar3d(i,j,k) enddo enddo enddo do k = levels(w_ind-1)+1, levels(w_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : w ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo deallocate(wworkvar3d) endif if (tv_ind > 0 .or. tsen_ind > 0) then allocate(tsenworkvar3d(nx_res,ny_res,nlevs)) varstrname = 'T' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d) varstrname = 'sphum' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 if (ql_ind > 0) then varstrname = 'liq_wat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(ql_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(ql_ind-1)+1, levels(ql_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : ql ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif if (qi_ind > 0) then varstrname = 'ice_wat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qi_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(qi_ind-1)+1, levels(qi_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : qi ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif if (qr_ind > 0) then varstrname = 'rainwat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qr_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(qr_ind-1)+1, levels(qr_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : qr ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif if (qs_ind > 0) then varstrname = 'snowwat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qs_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(qs_ind-1)+1, levels(qs_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : qs ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif if (qg_ind > 0) then varstrname = 'graupel' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qg_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(qg_ind-1)+1, levels(qg_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : qg ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif if (qnr_ind > 0) then varstrname = 'rain_nc' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qnr_ind-1)+k,nb,ne)=workvar3d(i,j,k) enddo enddo enddo do k = levels(qnr_ind-1)+1, levels(qnr_ind) if (nproc .eq. 0) & write(6,*) 'READFVregional : qnr ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo endif ! 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)) call fv3lamfile%get_idfn('delp',file_id,fv3filename) call read_fv3_restart_data3d('delp',fv3filename,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 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 if(l_fv3reg_filecombined) then call nc_check( nf90_close(file_id),& myname_,'close '//trim(filename) ) else do ifile=1,2 file_id=fv3lamfile%fv3lam_fileid(ifile) filename=fv3lamfile%fv3lamfilename(ifile) call nc_check( nf90_close(file_id),& myname_,'close '//trim(filename) ) enddo 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 FV3-Lam state or control vector !------------------------------------------------------------------------- !======================================================================== ! writegriddata.f90: write FV3-LAM 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,fv3filename1 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) :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind, qnr_ind integer(i_kind) file_id,file_id1 real(r_single), dimension(:,:), allocatable ::pswork real(r_single), dimension(:,:,:), allocatable ::workvar3d,workinc3d,workinc3d2,uworkvar3d,& vworkvar3d,tvworkvar3d,tsenworkvar3d,& workprsi,qworkvar3d,wworkvar3d real(r_single) :: clip !---------------------------------------------------------------------- ! 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,ifile,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 oz_ind = getindex(vars3d, 'oz') ! Oz (3D) 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 FV3 qi_ind = getindex(vars3d, 'qi') ! QI (3D) for FV3 qr_ind = getindex(vars3d, 'qr') ! QR (3D) for FV3 qs_ind = getindex(vars3d, 'qs') ! QS (3D) for FV3 qg_ind = getindex(vars3d, 'qg') ! QG (3D) for FV3 qnr_ind = getindex(vars3d, 'qnr') ! QNR (3D) for FV3 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) if(l_fv3reg_filecombined) then fv3filename=trim(adjustl(filename))//"_dynvartracer" call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_write,file_id),& myname_,'open: '//trim(adjustl(fv3filename)) ) call fv3lamfile%setupfile(fileid1=file_id,fv3fn1=trim(adjustl(fv3filename))) else fv3filename=trim(adjustl(filename))//"_dynvars" call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_write,file_id),& myname_,'open: '//trim(adjustl(fv3filename)) ) fv3filename1=trim(adjustl(filename))//"_tracer" call nc_check( nf90_open(trim(adjustl(fv3filename1)),nf90_write,file_id1),& myname_,'open: '//trim(adjustl(fv3filename1)) ) call fv3lamfile%setupfile(fileid1=file_id,fv3fn1=trim(adjustl(fv3filename)) , & fileid2=file_id1,fv3fn2=trim(adjustl(fv3filename1)) ) endif !---------------------------------------------------------------------- ! read u-component ! update CWM for WRF-NMM if (u_ind > 0) then varstrname = 'u' allocate(uworkvar3d(nx_res,ny_res+1,nlevs)) call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 (w_ind > 0) then varstrname = 'W' allocate(wworkvar3d(nx_res,ny_res,nlevs)) call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,wworkvar3d) 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(w_ind-1)+k,nb,ne) enddo enddo enddo wworkvar3d(1:nx_res,:,:)=wworkvar3d(1:nx_res,:,:)+workinc3d wworkvar3d(nx_res+1,:,:)=wworkvar3d(nx_res,:,:) call write_fv3_restart_data3d(varstrname,fv3filename,file_id,wworkvar3d) deallocate(wworkvar3d) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d) varstrname = 'sphum' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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 (ql_ind > 0) then varstrname = 'liq_wat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(ql_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d) endif if (qi_ind > 0) then varstrname = 'ice_wat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qi_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d) endif if (qr_ind > 0) then varstrname = 'rainwat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qr_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d) endif if (qs_ind > 0) then varstrname = 'snowwat' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qs_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d) endif if (qg_ind > 0) then varstrname = 'graupel' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qg_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if call write_fv3_restart_data3d(varstrname,fv3filename,file_id,workvar3d) endif if (qnr_ind > 0) then varstrname = 'rain_nc' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) 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(qnr_ind-1)+k,nb,ne) enddo enddo enddo workvar3d=workvar3d+workinc3d if ( l_use_enkf_directZDA .and. cliptracers ) then ! set cliptracers to remove negative hydrometers clip = tiny(workvar3d(1,1,1)) where (workvar3d < clip) workvar3d = clip end if 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 fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,filename,file_id,workvar3d) ! Pascal !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) end if if(l_fv3reg_filecombined) then call nc_check( nf90_close(file_id),& myname_,'close '//trim(filename) ) else do ifile=1,2 file_id=fv3lamfile%fv3lam_fileid(ifile) filename=fv3lamfile%fv3lamfilename(ifile) call nc_check( nf90_close(file_id),& myname_,'close '//trim(filename) ) enddo endif !---------------------------------------------------------------------- ! 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 subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) !Dummy subroutine declaration in place of the actual subroutine definition in !the GFS EnKF !to be implemented in the future use constants, only: max_varname_length use params, only: nbackgrounds implicit none integer, intent(in) :: nanal1,nanal2 character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d integer, intent(in) :: n2d,n3d,ndim integer, dimension(0:n3d), intent(in) :: levels real(r_single), dimension(npts,ndim,nbackgrounds,1), intent(inout) :: grdin logical, intent(in) :: no_inflate_flag end subroutine writeincrement subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) !Dummy subroutine declaration in place of the actual subroutine definition in !the GFS EnKF !to be implemented in the future use constants, only: max_varname_length use params, only: nbackgrounds implicit none character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d integer, intent(in) :: n2d,n3d,ndim integer, dimension(0:n3d), intent(in) :: levels real(r_single), dimension(npts,ndim,nbackgrounds,1), intent(inout) :: grdin logical, intent(in) :: no_inflate_flag end subroutine writeincrement_pnc subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & fileprefixes,filesfcprefixes,reducedgrid,grdin,qsat) !Dummy subroutine declaration in place of the actual subroutine definition in !the GFS EnKF !to be implemented in the future use constants, only: max_varname_length implicit none character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d integer, intent(in) :: n2d, n3d integer, dimension(0:n3d), intent(in) :: levels integer, intent(in) :: ndim, ntimes character(len=120), dimension(7), intent(in) :: fileprefixes character(len=120), dimension(7), intent(in) :: filesfcprefixes logical, intent(in) :: reducedgrid real(r_single), dimension(npts,ndim,ntimes,1), intent(out) :: grdin real(r_double), dimension(npts,nlevs,ntimes,1), intent(out) :: qsat end subroutine readgriddata_pnc subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) !Dummy subroutine declaration in place of the actual subroutine definition in !the GFS EnKF !to be implemented in the future use constants, only: max_varname_length use params, only: nbackgrounds implicit none character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d integer, intent(in) :: n2d,n3d,ndim integer, dimension(0:n3d), intent(in) :: levels real(r_single), dimension(npts,ndim,nbackgrounds,1), intent(inout) :: grdin logical, intent(in) :: no_inflate_flag end subroutine writegriddata_pnc subroutine type_bound_setupfile(this,fileid1,fv3fn1,fileid2,fv3fn2) class (type_fv3lamfile) :: this integer(i_kind) fileid1 integer(i_kind), optional :: fileid2 character(len=*)::fv3fn1 character(len=*),optional ::fv3fn2 if (present (fileid2)) then this%l_filecombined=.false. this%fv3lamfilename(1)=trim(fv3fn1) this%fv3lamfilename(2)=trim(fv3fn2) this%fv3lam_fileid(1)=fileid1 this%fv3lam_fileid(2)=fileid2 else this%l_filecombined=.true. this%fv3lamfilename(1)=fv3fn1 this%fv3lam_fileid(1)=fileid1 endif end subroutine type_bound_setupfile subroutine type_bound_getidfn(this,vnamloc,fileid,fv3fn) class (type_fv3lamfile) :: this integer(i_kind) fileid character(len=*)::fv3fn,vnamloc if (.not.this%l_filecombined) then if(ifindstrloc(vardynvars,vnamloc)> 0) then fv3fn=trim(this%fv3lamfilename(1)) fileid=this%fv3lam_fileid(1) else if(ifindstrloc(vartracers,vnamloc)> 0) then fv3fn=trim(this%fv3lamfilename(2)) fileid=this%fv3lam_fileid(2) else write(6,*)"the varname ",trim(vnamloc)," is not recognized in the ype_bound_getidfn, stop" call stop2(23) endif else fv3fn=trim(this%fv3lamfilename(1)) fileid=this%fv3lam_fileid(1) endif end subroutine type_bound_getidfn function ifindstrloc(str_array,strin) integer(i_kind) ifindstrloc character(len=max_varname_length),dimension(:) :: str_array character(len=*) :: strin integer(i_kind) i ifindstrloc=0 do i=1,size(str_array) if(trim(str_array(i)) == trim(strin)) then ifindstrloc=i exit endif enddo end function ifindstrloc end module gridio