diff --git a/fix b/fix index 384204259f..f0aa67e488 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 384204259f72b55ae8dd1b3e6959b5bc735ce904 +Subproject commit f0aa67e488cb343b681f078639bf3205b0d6b663 diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90 index 9e1f3310ed..a50355d4e8 100644 --- a/src/enkf/gridinfo_fv3reg.f90 +++ b/src/enkf/gridinfo_fv3reg.f90 @@ -44,7 +44,7 @@ module gridinfo use mpisetup, only: nproc, mpi_integer, mpi_real4, mpi_comm_world use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio, fgfileprefixes, & - fv3fixpath, nx_res,ny_res, ntiles + fv3fixpath, nx_res,ny_res, ntiles,l_fv3reg_filecombined use kinds, only: r_kind, i_kind, r_double, r_single use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length use constants, only: half @@ -91,7 +91,7 @@ subroutine getgridinfo(fileprefix, reducedgrid) real(r_kind), allocatable, dimension(:,:) :: pressimn,presslmn real(r_kind) kap,kapr,kap1 -integer(i_kind) file_id,var_id,dim_id,nlevsp1,nx_tile,ny_tile,ntile +integer(i_kind) :: file_id,var_id,dim_id,nlevsp1,nx_tile,ny_tile,ntile integer (i_kind):: nn_tile0 integer(i_kind) :: nlevsp1n real(r_single), allocatable, dimension(:,:) :: lat_tile,lon_tile,ps @@ -212,7 +212,11 @@ subroutine getgridinfo(fileprefix, reducedgrid) nn_tile0=(ntile-1)*nx_res*ny_res nn=nn_tile0 write(char_tile, '(i1)') ntile - filename = 'fv3sar_tile'//char_tile//"_ensmean_dynvartracer" + if(l_fv3reg_filecombined) then + filename = 'fv3sar_tile'//char_tile//"_ensmean_dynvartracer" + else + filename = 'fv3sar_tile'//char_tile//"_ensmean_dynvars" + endif !print *,trim(adjustl(filename)) call nc_check( nf90_open(trim(adjustl(filename)),nf90_nowrite,file_id),& myname_,'open: '//trim(adjustl(filename)) ) diff --git a/src/enkf/gridio_fv3reg.f90 b/src/enkf/gridio_fv3reg.f90 index 0dd03dfb82..930264bc26 100644 --- a/src/enkf/gridio_fv3reg.f90 +++ b/src/enkf/gridio_fv3reg.f90 @@ -31,11 +31,12 @@ module gridio !========================================================================= ! 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 + 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 @@ -52,7 +53,23 @@ module gridio public :: writeincrement, writeincrement_pnc !------------------------------------------------------------------------- - + + integer(i_kind) ,parameter:: ndynvarslist=6, ntracerslist=8 + character(len=max_varname_length), parameter :: vardynvars(ndynvarslist) =(/"u","v", & + "T","W","DZ","delp"/) + character(len=max_varname_length), parameter :: vartracers(ntracerslist) =(/'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 @@ -76,9 +93,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f ! Define local variables character(len=500) :: filename - character(len=:),allocatable :: fv3filename + character(len=:),allocatable :: fv3filename,fv3filename1 character(len=7) :: charnanal - integer(i_kind) file_id + integer(i_kind) file_id,file_id1 real(r_single), dimension(:,:,:), allocatable ::workvar3d,uworkvar3d,& vworkvar3d,tvworkvar3d,tsenworkvar3d,& workprsi,qworkvar3d,wworkvar3d @@ -98,7 +115,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f 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 + integer :: tmp_ind,ifile logical :: ice !====================================================================== @@ -141,372 +158,401 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f 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" + if (nanal > 0) then + write(charnanal,'(a3, i3.3)') 'mem', nanal + else + charnanal = 'ensmean' + endif - !---------------------------------------------------------------------- - ! read u-component - call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_nowrite,file_id),& - myname_,'open: '//trim(adjustl(fv3filename)) ) + do ntile=1,ntiles + nn_tile0=(ntile-1)*nx_res*ny_res + write(char_tile, '(i1)') ntile - !---------------------------------------------------------------------- - ! 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' + 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 - call read_fv3_restart_data3d(varstrname,fv3filename,file_id,uworkvar3d) + 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(u_ind-1)+k,nb,ne)=uworkvar3d(i,j,k) + 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 enddo - do k = levels(u_ind-1)+1, levels(u_ind) + do k = levels(v_ind-1)+1, levels(v_ind) if (nproc .eq. 0) & - write(6,*) 'READFVregional : u ', & - & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) + write(6,*) 'READFVregional : v ', & + & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo + deallocate(vworkvar3d) - 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 (w_ind > 0) then - allocate(wworkvar3d(nx_res,ny_res,nlevs)) - varstrname = 'W' - 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) + 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 - 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 + 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) - 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) + endif + if (tv_ind > 0.or.tsen_ind) 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 - if (q_ind > 0) then - varstrname = 'sphum' - do k=1,nlevs - nn = nn_tile0 + endif + if(tv_ind > 0) then + do k=1,nlevs 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 + do i=1,nx_res + workvar3d(i,j,k)=tsenworkvar3d(i,j,k)*(one+fv*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 + 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 - 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 + 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 - 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 - endif - if(allocated(tsenworkvar3d)) deallocate(tsenworkvar3d) - + 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) + + 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 - 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 + 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 + endif - if (ql_ind > 0) then - varstrname = 'liq_wat' - 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) + 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 - 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 + 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 + endif - if (qi_ind > 0) then - varstrname = 'ice_wat' - 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) + 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 - 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 + endif - if (qr_ind > 0) then - varstrname = 'rainwat' - 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) + 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 - 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 + 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 + endif - if (qs_ind > 0) then - varstrname = 'snowwat' - 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) + 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 - 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 + endif - if (qg_ind > 0) then - varstrname = 'graupel' - 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) + 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 - 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 + endif - if (qnr_ind > 0) then - varstrname = 'rain_nc' - 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) + 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 - 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 + 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 - - 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 + 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)) - 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) + !---------------------------------------------------------------------- + ! 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) + 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 - 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 + + 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 - 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) + 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 - 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 + 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 @@ -548,7 +594,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid !---------------------------------------------------------------------- ! Define variables computed within subroutine character(len=500) :: filename - character(len=:),allocatable :: fv3filename + character(len=:),allocatable :: fv3filename,fv3filename1 character(len=7) :: charnanal !---------------------------------------------------------------------- @@ -556,7 +602,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 + 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,& @@ -575,7 +621,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid !---------------------------------------------------------------------- ! Define counting variables - integer :: i,j,k,nn,ntile,nn_tile0, nb,ne,nanal + integer :: i,j,k,ifile,nn,ntile,nn_tile0, nb,ne,nanal @@ -611,10 +657,10 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid call stop2(23) endif ne = 0 - ensmemloop: do nanal=nanal1,nanal2 + ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 - backgroundloop: do nb=1,nbackgrounds + 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)) @@ -633,13 +679,26 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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" + 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 - call nc_check( nf90_open(trim(adjustl(fv3filename)),nf90_write,file_id),& - myname_,'open: '//trim(adjustl(fv3filename)) ) ! update CWM for WRF-NMM @@ -647,6 +706,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -668,6 +728,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -689,6 +750,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -719,6 +781,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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) @@ -735,25 +798,28 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid varstrname = 'T' allocate(tsenworkvar3d(nx_res,ny_res,nlevs)) - call read_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d) + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + 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 + 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 + 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) & @@ -767,6 +833,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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) & @@ -784,6 +851,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -801,7 +869,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -824,6 +892,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -846,6 +915,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -868,6 +938,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -890,6 +961,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -911,7 +983,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -935,6 +1007,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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 @@ -972,26 +1045,35 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid call write_fv3_restart_data3d(varstrname,filename,file_id,workvar3d) end if - call nc_check( nf90_close(file_id),& - myname_,'close '//trim(filename) ) + 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 + 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) + 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 + end do backgroundloop ! loop over backgrounds to read in + end do ensmemloop ! loop over ens members to read in ! Return calculated values @@ -1064,6 +1146,58 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ 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 diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index e46bfa195b..fc44da61c9 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -243,6 +243,7 @@ module params integer(i_kind),public :: ntiles=6 integer(i_kind),public :: nx_res=0,ny_res=0 logical,public ::l_pres_add_saved +logical,public ::l_fv3reg_filecombined =.true. !=.true., the dynvar and tracer files would be combined for enkf fv3_reg ! for parallel netCDF logical, public :: paranc = .false. @@ -286,7 +287,7 @@ module params lnsigcutoffrdrnh,lnsigcutoffrdrsh,lnsigcutoffrdrtr,& l_use_enkf_directZDA namelist /nam_wrf/arw,nmm,nmm_restart -namelist /nam_fv3/fv3fixpath,nx_res,ny_res,ntiles,l_pres_add_saved +namelist /nam_fv3/fv3fixpath,nx_res,ny_res,ntiles,l_pres_add_saved,l_fv3reg_filecombined namelist /satobs_enkf/sattypes_rad,dsis namelist /ozobs_enkf/sattypes_oz