From 45e3e6cc7020cc9495e78ff03f74db52faad6b7b Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Mon, 11 Dec 2023 20:52:54 -0600 Subject: [PATCH] Squash HAFS EnKF DA related updates developed mainly by @TingLei-NOAA. --- src/enkf/enkf_main.f90 | 8 +- src/enkf/gridinfo_fv3reg.f90 | 6 +- src/enkf/gridio_fv3reg.f90 | 609 +++++++++++++++++++++-------------- src/enkf/loadbal.f90 | 70 +++- src/enkf/params.f90 | 11 +- 5 files changed, 441 insertions(+), 263 deletions(-) diff --git a/src/enkf/enkf_main.f90 b/src/enkf/enkf_main.f90 index 1a38dd7525..eb4e89d90e 100644 --- a/src/enkf/enkf_main.f90 +++ b/src/enkf/enkf_main.f90 @@ -76,7 +76,7 @@ program enkf_main ! reads namelist parameters. use params, only : read_namelist,cleanup_namelist,letkf_flag,readin_localization,lupd_satbiasc,& numiter, nanals, lupd_obspace_serial, write_spread_diag, & - lobsdiag_forenkf, netcdf_diag, efsoi_cycling, ntasks_io + lobsdiag_forenkf, ldo_enscalc_option, netcdf_diag, efsoi_cycling, ntasks_io ! mpi functions and variables. use mpisetup, only: mpi_initialize, mpi_initialize_io, mpi_cleanup, nproc, & mpi_wtime @@ -154,7 +154,7 @@ program enkf_main ! read obs, initial screening. t1 = mpi_wtime() - call readobs() + if(ldo_enscalc_option==0) call readobs() t2 = mpi_wtime() if (nproc == 0) print *,'time in read_obs =',t2-t1,'on proc',nproc @@ -167,10 +167,12 @@ program enkf_main endif ! print innovation statistics for prior on root task. + if(ldo_enscalc_option==0) then if (nproc == 0) then print *,'innovation statistics for prior:' call print_innovstats(obfit_prior, obsprd_prior) end if + end if ! read state/control vector info from anavinfo call init_controlvec() @@ -204,6 +206,7 @@ program enkf_main t1 = mpi_wtime() ! state and bias correction coefficient update iteration. + if(ldo_enscalc_option ==0 ) then if(letkf_flag) then ! do ob space update using serial filter if desired if (lupd_obspace_serial) call enkf_update() @@ -257,6 +260,7 @@ program enkf_main call obsmod_cleanup() t1 = mpi_wtime() + end if ! ldo_enscalc_option call gather_chunks() t2 = mpi_wtime() if (nproc == 0) print *,'time in gather_chunks =',t2-t1,'on proc',nproc diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90 index 4eff63c003..778a73d89a 100644 --- a/src/enkf/gridinfo_fv3reg.f90 +++ b/src/enkf/gridinfo_fv3reg.f90 @@ -72,10 +72,10 @@ module gridinfo integer(i_kind),public :: npts integer(i_kind),public :: ntrunc ! supported variable names in anavinfo -character(len=max_varname_length),public, dimension(16) :: & +character(len=max_varname_length),public, dimension(17) :: & vars3d_supported = [character(len=max_varname_length) :: & - 'u', 'v', 'w', 't', 'q', 'oz', 'cw', 'tsen', 'prse', & - 'ql', 'qi', 'qr', 'qs', 'qg', 'qnr','dbz'] + 'u', 'v', 'w', 't', 'q', 'oz', 'cw', 'tsen', 'prse', 'delp', & + 'ql', 'qi', 'qr', 'qs', 'qg', 'qnr', 'dbz'] character(len=max_varname_length),public, dimension(3) :: & vars2d_supported = [character(len=max_varname_length) :: & 'ps', 'pst', 'sst'] diff --git a/src/enkf/gridio_fv3reg.f90 b/src/enkf/gridio_fv3reg.f90 index 068e6cba8b..203c87eacc 100644 --- a/src/enkf/gridio_fv3reg.f90 +++ b/src/enkf/gridio_fv3reg.f90 @@ -3,18 +3,14 @@ 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) + ! variables from FV3-LAM warmstart files ! ! prgmmr: Winterbottom org: ESRL/PSD1 date: 2011-11-30 ! ! program history log: - ! + ! ! 2011-11-30 Winterbottom - Initial version. ! ! 2019-01- Ting -- modified for fv3sar @@ -42,7 +38,7 @@ module gridio use params, only: nlevs, cliptracers, datapath, arw, nmm, datestring use params, only: nx_res,ny_res,nlevs,ntiles,l_fv3reg_filecombined,& fv3_io_layout_nx,fv3_io_layout_ny,nanals - use params, only: pseudo_rh + use params, only: pseudo_rh,ldo_enscalc_option 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 @@ -87,7 +83,7 @@ module gridio 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 constants, only:izero,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 @@ -119,18 +115,19 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ! 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(i_kind) :: nlevsp1 integer (i_kind):: i,j, k,nn,ntile,nn_tile0, nb,nanal,ne - integer(i_kind) :: u_ind, v_ind, tv_ind,tsen_ind, q_ind, oz_ind - integer(i_kind) :: dbz_ind - integer(i_kind) :: w_ind, ql_ind, qi_ind, qr_ind, qs_ind, qg_ind, qnr_ind - integer (i_kind):: ps_ind, sst_ind + integer(i_kind) :: u_ind=izero, v_ind=izero, tv_ind,tsen_ind=izero, & + q_ind=izero,delp_ind=izero, oz_ind=izero + integer(i_kind) :: dbz_ind=izero + integer(i_kind) :: w_ind=izero, ql_ind=izero, qi_ind=izero, qr_ind=izero,& + qs_ind=izero, qg_ind=izero, qnr_ind=izero + integer (i_kind):: ps_ind=izero, sst_ind=izero integer (i_kind):: tmp_ind,ifile logical (i_kind):: ice @@ -143,6 +140,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, w_ind = getindex(vars3d, 'w') ! W (3D) tv_ind = getindex(vars3d, 't') ! Tv (3D) q_ind = getindex(vars3d, 'q') ! Q (3D) + delp_ind = getindex(vars3d, 'delp') ! delp (3D) oz_ind = getindex(vars3d, 'oz') ! Oz (3D) tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) @@ -156,6 +154,17 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ps_ind = getindex(vars2d, 'ps') ! Ps (2D) sst_ind = getindex(vars2d, 'sst') ! SST (2D) + if(ps_ind > 0 .and. delp_ind > 0) then + write(6,*)'ps_ind and delp_ind should be >0 at the same time, stop' + call stop2(23) + endif + if ( ldo_enscalc_option > 0 ) then + if (pseudo_rh) then + write(6,*)'For ensemble mean or recenter tasks, pseudo_rh should be false' + call stop2(23) + endif + endif + ! Initialize all constants required by routine allocate(workvar3d(nx_res,ny_res,nlevs)) @@ -282,34 +291,36 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, enddo 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) + if(q_ind > 0) then varstrname = 'sphum' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d) + 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 (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 (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) + if (.not. (q_ind > 0) ) then + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d) endif if(tv_ind > 0) then do k=1,nlevs @@ -521,67 +532,90 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, !---------------------------------------------------------------------- ! Allocate memory for variables computed within routine - - if (ps_ind > 0) then + + + if ((ps_ind > 0.or. delp_ind>0) .or. pseudo_rh) then !if q_ind > 0 qsat will be use when pseudo_rh=.true. + 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) + call read_fv3_restart_data3d('delp',fv3filename,file_id,workvar3d) + workvar3d=workvar3d*0.01_r_kind 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) + workprsi(:,:,i)=workvar3d(:,:,i)+workprsi(:,:,i+1) enddo - - pswork(:,:)=workprsi(:,:,1) + pswork(:,:)=workprsi(:,:,1) + if(ps_ind > 0) then + 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 + else if (delp_ind >0 ) then + do k=1,nlevs + nn = nn_tile0 + do j=1,ny_res + do i=1,nx_res + nn=nn+1 + vargrid(nn,levels(delp_ind-1)+k,nb,ne)=workvar3d(i,j,k) + enddo + enddo + enddo - 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 + endif - - - do k=1,nlevs - do j=1,ny_res + 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. - if (pseudo_rh) then + if( pseudo_rh) then + if(tv_ind < 0 .and. tsen_ind > 0) then + do k=1,nlevs + do j=1,ny_res + do i=1,nx_res + tvworkvar3d(i,j,k)=tsenworkvar3d(i,j,k)*(one+fv*qworkvar3d(i,j,k)) + enddo + enddo + enddo + elseif (tv_ind > 0) then + write(6,*)'tvworkvar3d has been calculated' + else + write(6,*)'One and only one of tv_ind and tsen_ind should be > 0,stop ' + call stop2(23) + endif + + ice=.true. 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 + endif + if(.not. pseudo_rh ) qsatworkvar3d(:,:,:) = 1._r_double + 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) - 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) ) @@ -604,20 +638,15 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, return -end subroutine readgriddata + end subroutine readgriddata !======================================================================== - ! readgriddata_nmm.f90: read FV3-Lam state or control vector + ! writegriddata: write fv3-lam !------------------------------------------------------------------------- - - !======================================================================== - ! 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 + subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflate_flag) + use constants, only: izero,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 @@ -643,9 +672,10 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid character(len=7) :: charnanal !---------------------------------------------------------------------- - integer(i_kind) :: u_ind, v_ind, tv_ind, tsen_ind,q_ind, ps_ind,oz_ind,dbz_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) :: u_ind=izero, v_ind=izero, tv_ind=izero, tsen_ind=izero, q_ind=izero,& + delp_ind=izero, ps_ind=izero, oz_ind=izero, dbz_ind=izero + integer(i_kind) :: w_ind=izero, cw_ind=izero, ph_ind=izero + integer(i_kind) :: ql_ind=izero, qi_ind=izero, qr_ind=izero, qs_ind=izero, qg_ind=izero, qnr_ind=izero integer(i_kind) file_id,file_id1,file_id2 real(r_single), dimension(:,:), allocatable ::pswork @@ -663,15 +693,13 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid character(len=12) :: varstrname character(len=1) char_tile character(len=24),parameter :: myname_ = 'fv3: writegriddata' - + logical :: l_q_updated !---------------------------------------------------------------------- ! Define counting variables integer(i_kind) :: i,j,k,ifile,nn,ntile,nn_tile0, nb,ne,nanal - - write(6,*)"anlfileprefixes, fgfileprefixes are not used in the current implementation", & - anlfileprefixes, fgfileprefixes + anlfileprefixes, fgfileprefixes write(6,*)"the no_inflate_flag is not used in the currrent implementation ",no_inflate_flag !---------------------------------------------------------------------- nlevsp1=nlevs+1 @@ -682,6 +710,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid tsen_ind = getindex(vars3d, 'tsen') ! Tv (3D) q_ind = getindex(vars3d, 'q') ! Q (3D) oz_ind = getindex(vars3d, 'oz') ! Oz (3D) + delp_ind = getindex(vars3d, 'delp') ! delp (3D) w_ind = getindex(vars3d, 'w') ! W for WRF-ARW ql_ind = getindex(vars3d, 'ql') ! QL (3D) for FV3 @@ -818,7 +847,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid endif - if (tv_ind > 0.or.tsen_ind>0 ) then + if (tv_ind > 0.or.tsen_ind>0 ) then varstrname = 'T' if(tsen_ind>0) then @@ -850,9 +879,9 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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) + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d) !enforce lower positive bound (clip) to replace negative hydrometers if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip tvworkvar3d=tsenworkvar3d*(one+fv*qworkvar3d) @@ -875,6 +904,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid varstrname = 'T' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call write_fv3_restart_data3d(varstrname,fv3filename,file_id,tsenworkvar3d) + l_q_updated=.true. do k=1,nlevs if (nproc .eq. 0) & write(6,*) 'WRITEregional : T ', & @@ -882,20 +912,6 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid 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) @@ -903,6 +919,35 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid endif endif + if(q_ind>0) then + if(.not.l_q_updated ) then + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,qworkvar3d) + 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 + !enforce lower positive bound (clip) to replace negative q + if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip + endif + 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 + + if (oz_ind > 0) then varstrname = 'o3mr' @@ -1059,47 +1104,66 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid endif - if (ps_ind > 0) then + if (ps_ind > 0 .or. delp_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 - 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 + workvar3d=workvar3d*0.01_r_kind + if(ps_ind > 0) then + workprsi(:,:,nlevsp1)=eta1_ll(nlevsp1) !etal_ll is needed + do i=nlevs,1,-1 + workprsi(:,:,i)=workvar3d(:,:,i)+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 - workprsi=workprsi+workinc3d2 - else - workprsi(:,:,1)=workprsi(:,:,1)+pswork - do k=2,nlevsp1 - workprsi(:,:,k)=eta1_ll(k)+eta2_ll(k)*workprsi(:,:,1) + elseif (delp_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(delp_ind-1)+k,nb,ne) + enddo + enddo + enddo + workvar3d=(workvar3d+workinc3d)*100.0_r_kind + + do k=1,nlevs + if (nproc .eq. 0) & + write(6,*) 'WRITEregional : delp ', & + & k, minval(workvar3d(:,:,k)), maxval(workvar3d(:,:,k)) enddo - endif - do k=1,nlevs - workvar3d(:,:,k)=(workprsi(:,:,k)-workprsi(:,:,k+1))*100.0 - enddo + endif 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) ) @@ -1203,7 +1267,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & character(len=1) char_tile character(len=24) strsubid integer(i_kind) :: nlevsp1 - integer(i_kind) :: u_ind, v_ind, tv_ind,tsen_ind, q_ind, oz_ind + integer(i_kind) :: u_ind, v_ind, tv_ind,tsen_ind, delp_ind,q_ind, oz_ind integer(i_kind) :: w_ind, ql_ind, qi_ind, qr_ind, qs_ind, qg_ind, qnr_ind integer(i_kind) :: ps_ind, sst_ind integer(i_kind) :: tmp_ind,ifile @@ -1221,6 +1285,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & w_ind = getindex(vars3d, 'w') ! W (3D) tv_ind = getindex(vars3d, 't') ! Tv (3D) q_ind = getindex(vars3d, 'q') ! Q (3D) + delp_ind = getindex(vars3d, 'delp') ! delp (3D) oz_ind = getindex(vars3d, 'oz') ! Oz (3D) tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) @@ -1418,6 +1483,33 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & write(6,*) 'READFVregional : w ', & & k, minval(vargrid(:,k,nb,ne)), maxval(vargrid(:,k,nb,ne)) enddo + endif !iope == 0 + + endif + if (q_ind > 0) then + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) + do k=1,nlevs + call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), & + mpi_real4, qworkvar3d(:,:,k), recvcounts2d, displs2d,& + mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) + enddo + if(iope==0) then + 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 !iope == 0 endif @@ -1434,48 +1526,33 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & varstrname = 'T' call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,loctsenworkvar3d) - varstrname = 'sphum' - call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) do k=1,nlevs call mpi_gatherv(loctsenworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4,& tsenworkvar3d(:,:,k), recvcounts2d, displs2d,& mpi_real4, 0,iocomms(mem_pe(nproc)),iret) enddo - do k=1,nlevs - call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4, & - qworkvar3d(:,:,k), recvcounts2d, displs2d,& - mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) - enddo + if( .not. (q_ind>0) ) then + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) + do k=1,nlevs + call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4, & + qworkvar3d(:,:,k), recvcounts2d, displs2d,& + mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) + enddo + endif if(iope== 0) then - 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 + do k=1,nlevs + do j=1,ny_res + do i=1,nx_res + tvworkvar3d(i,j,k)=tsenworkvar3d(i,j,k)*(one+fv*qworkvar3d(i,j,k)) + enddo + enddo + enddo 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 + workvar3d=tvworkvar3d else! tsen_id >0 workvar3d=tsenworkvar3d endif @@ -1705,13 +1782,14 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & !---------------------------------------------------------------------- ! Allocate memory for variables computed within routine - if (ps_ind > 0) then + if (ps_ind > 0.or. delp_ind >0 .or. q_ind > 0) then !if q_ind > 0 qsat will be use when pseudo_rh=.true. if (iope ==0) then allocate(workprsi(nx_res,ny_res,nlevsp1)) allocate(pswork(nx_res,ny_res)) endif call fv3lamfile%get_idfn('delp',file_id,fv3filename) call read_fv3_restart_data3d('delp',fv3filename,file_id,locworkvar3d) + locworkvar3d=locworkvar3d*0.01_r_kind do k=1,nlevs call mpi_gatherv(locworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4,& workvar3d(:,:,k), recvcounts2d, displs2d,& @@ -1720,29 +1798,43 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & if(iope == 0 ) then 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) + workprsi(:,:,i)=workvar3d(:,:,i)+workprsi(:,:,i+1) enddo pswork(:,:)=workprsi(:,:,1) + if(ps_ind >0) then + 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 + if(nproc==0) then + write(6,*) 'READFVregional : ps ', & + & k, minval(vargrid(:,levels(n3d)+ps_ind,nb,ne)), maxval(vargrid(:,k,nb,ne)) - 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 + endif !nproc==0 + elseif (delp_ind >0) then - if(nproc==0) then - write(6,*) 'READFVregional : ps ', & - & k, minval(vargrid(:,levels(n3d)+ps_ind,nb,ne)), maxval(vargrid(:,k,nb,ne)) + do k=1,nlevs + nn = nn_tile0 + do j=1,ny_res + do i=1,nx_res + nn=nn+1 + vargrid(nn,levels(delp_ind-1)+k,nb,ne)=workvar3d(i,j,k) + enddo + enddo + if(nproc==0) then + write(6,*) 'READFVregional : delp ', & + & k, minval(vargrid(:,levels(n3d)+delp_ind,nb,ne)), maxval(vargrid(:,k,nb,ne)) - endif !nproc==0 - + endif !nproc==0 + enddo + endif - do k=1,nlevs do j=1,ny_res do i=1,nx_res @@ -1845,7 +1937,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat integer(i_kind) :: nlevsp1 integer(i_kind) :: i,j,ii, k,nn,ntile,nn_tile0, nb,nanal,ne,iret,ierror integer(i_kind) :: irec,isend - integer(i_kind) :: u_ind, v_ind, tv_ind,tsen_ind, q_ind, oz_ind + integer(i_kind) :: u_ind, v_ind, tv_ind,tsen_ind,delp_ind, q_ind, oz_ind integer(i_kind) :: w_ind, ql_ind, qi_ind, qr_ind, qs_ind, qg_ind, qnr_ind integer (i_kind):: ps_ind, sst_ind integer (i_kind):: tmp_ind,ifile @@ -1862,6 +1954,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat tv_ind = getindex(vars3d, 't') ! Tv (3D) tsen_ind = getindex(vars3d, 'tsen') ! Tv (3D) q_ind = getindex(vars3d, 'q') ! Q (3D) + delp_ind = getindex(vars3d, 'delp') ! delp (3D) oz_ind = getindex(vars3d, 'oz') ! Oz (3D) w_ind = getindex(vars3d, 'w') ! W for WRF-ARW @@ -1897,7 +1990,6 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat allocate(qworkvar3d(1,1,nlevs)) endif - i=ij_pe(1,iope) j=ij_pe(2,iope) nxloc=nxlocgroup(i,j) @@ -2091,6 +2183,35 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat enddo call write_fv3_restart_data3d(varstrname,fv3filename,file_id,locworkvar3d) + endif + if (q_ind > 0) then + varstrname = 'sphum' + + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) + do k=1,nlevs + call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4,& + qworkvar3d(:,:,k), recvcounts2d, displs2d,& + mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) + enddo + if(iope == 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(1:nx_res,:,:)+workinc3d + endif + do k=1,nlevs + call mpi_scatterv(qworkvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,& + locqworkvar3d(:,:,k),recvcounts2d(iope+1),mpi_real4,0,iocomms(mem_pe(nproc)),ierror) + enddo + call write_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) + endif if (tv_ind > 0.or.tsen_ind>0 ) then @@ -2149,14 +2270,16 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat tsenworkvar3d(:,:,k), recvcounts2d, displs2d,& mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) enddo - varstrname = 'sphum' - call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) - call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) - do k=1,nlevs - call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4, & - qworkvar3d(:,:,k), recvcounts2d, displs2d,& - mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) - enddo + if(.not. (q_ind > 0) ) then + varstrname = 'sphum' + call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) + call read_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) + do k=1,nlevs + call mpi_gatherv(locqworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4, & + qworkvar3d(:,:,k), recvcounts2d, displs2d,& + mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) + enddo + endif if(iope ==0 ) then if ( cliptracers ) where (qworkvar3d < clip) qworkvar3d = clip @@ -2192,25 +2315,6 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat enddo endif - - - if(q_ind>0) then - varstrname='sphum' - - call fv3lamfile%get_idfn(varstrname,file_id,fv3filename) - do k=1,nlevs - call mpi_scatterv(qworkvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,& - locqworkvar3d(:,:,k),recvcounts2d(iope+1),mpi_real4,0,iocomms(mem_pe(nproc)),ierror) - enddo - call write_fv3_restart_data3d(varstrname,fv3filename,file_id,locqworkvar3d) - if(nproc ==0 ) then - do k=1,nlevs - write(6,*) 'WRITEregional : sphum ', & - & k, minval(qworkvar3d(:,:,k)), maxval(qworkvar3d(:,:,k)) - enddo - endif - - endif !q_ind deallocate(tsenworkvar3d) if (allocated(loctsenworkvar3d)) deallocate(loctsenworkvar3d) deallocate(locqworkvar3d) @@ -2432,52 +2536,68 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,vargrid,no_inflat call write_fv3_restart_data3d(varstrname,fv3filename,file_id,locworkvar3d) endif - if (ps_ind > 0) then + if (ps_ind > 0.or.delp_ind > 0) then if(iope == 0) allocate(workprsi(nx_res,ny_res,nlevsp1)) if(iope == 0) 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,locworkvar3d) ! Pascal + locworkvar3d=locworkvar3d*0.01_r_kind do k=1,nlevs call mpi_gatherv(locworkvar3d(1:nxloc,1:nyloc,k), recvcounts2d(iope+1), mpi_real4, & workvar3d(:,:,k), recvcounts2d, displs2d,& mpi_real4, 0,iocomms(mem_pe(nproc)) ,iret) enddo if(iope == 0 ) then - 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 + if(ps_ind > 0 ) then + workprsi(:,:,nlevsp1)=eta1_ll(nlevsp1) !etal_ll is needed + do i=nlevs,1,-1 + workprsi(:,:,i)=workvar3d(:,:,i)+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) + 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 - 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) + 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_r_kind enddo - workprsi=workprsi+workinc3d2 - else - workprsi(:,:,1)=workprsi(:,:,1)+pswork - do k=2,nlevsp1 - workprsi(:,:,k)=eta1_ll(k)+eta2_ll(k)*workprsi(:,:,1) + else !delp_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(delp_ind-1)+k,nb,ne) enddo - endif - do k=1,nlevs - workvar3d(:,:,k)=(workprsi(:,:,k)-workprsi(:,:,k+1))*100.0_r_kind - enddo + enddo + enddo + workvar3d=(workvar3d+workinc3d) + endif + workvar3d=workvar3d*100.0_r_kind + endif !iope==0 do k=1,nlevs - call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,& - locworkvar3d(:,:,k),recvcounts2d(iope+1),mpi_real4,0,iocomms(mem_pe(nproc)),ierror) + call mpi_scatterv(workvar3d(:,:,k),recvcounts2d,displs2d,mpi_real4,& + locworkvar3d(:,:,k),recvcounts2d(iope+1),mpi_real4,0,iocomms(mem_pe(nproc)),ierror) enddo @@ -2579,5 +2699,4 @@ function ifindstrloc(str_array,strin) end function ifindstrloc - end module gridio diff --git a/src/enkf/loadbal.f90 b/src/enkf/loadbal.f90 index 3292c99ee4..7109cb2751 100644 --- a/src/enkf/loadbal.f90 +++ b/src/enkf/loadbal.f90 @@ -142,6 +142,7 @@ subroutine load_balance() ! Uses "Graham's rule", which simply ! stated, assigns each new work item to the task that currently has the ! smallest load. +use params, only: ldo_enscalc_option implicit none integer(i_kind), allocatable, dimension(:) :: rtmp,numobs !real(r_single), allocatable, dimension(:) :: buffer @@ -179,6 +180,7 @@ subroutine load_balance() numptsperproc = 0 np = 0 test_loadbal = .false. ! simple partition for testing +if(ldo_enscalc_option .ne. 0 ) test_loadbal=.true. do n=1,npts if (test_loadbal) then ! use simple partition (does not use estimated workload) for testing @@ -240,7 +242,8 @@ subroutine load_balance() end do ! for serial filter, partition obs for observation space update. -if (.not. letkf_flag .or. lupd_obspace_serial) then +!cltorg if (.not. letkf_flag .or. lupd_obspace_serial) then +if ((.not. letkf_flag .or. lupd_obspace_serial).and.ldo_enscalc_option==0) then allocate(numobsperproc(numproc)) allocate(iprocob(nobstot)) ! default is to partition obs simply, since @@ -346,6 +349,7 @@ subroutine scatter_chunks ! decomposition from load_balance use controlvec, only: ncdim, grdin use params, only: nbackgrounds, ntasks_io, nanals_per_iotask +use params, only: ldo_enscalc_option implicit none integer(i_kind), allocatable, dimension(:) :: scounts, displs, rcounts @@ -475,23 +479,67 @@ subroutine scatter_chunks allocate(ensmean_chunk_prior(npts_max,ncdim,nbackgrounds)) ensmean_chunk = 0_r_single -!==> compute mean, remove it from anal_chunk +if(ldo_enscalc_option ==0) then !regular enkf run +!==> compute ensemble of first guesses on each task, remove mean from anal. !$omp parallel do schedule(dynamic,1) private(nn,i,n,nb) -do nb=1,nbackgrounds - do nn=1,ncdim - do i=1,numptsperproc(nproc+1) - ensmean_chunk(i,nn,nb) = sum(anal_chunk(:,i,nn,nb))/float(nanals) - ensmean_chunk_prior(i,nn,nb) = ensmean_chunk(i,nn,nb) - ! remove mean from ensemble. - do nanal=1,nanals - anal_chunk(nanal,i,nn,nb) = anal_chunk(nanal,i,nn,nb)-ensmean_chunk(i,nn,nb) - anal_chunk_prior(nanal,i,nn,nb)=anal_chunk(nanal,i,nn,nb) + do nb=1,nbackgrounds + do nn=1,ncdim + do i=1,numptsperproc(nproc+1) + ensmean_chunk(i,nn,nb) = sum(anal_chunk(:,i,nn,nb))/float(nanals) + ensmean_chunk_prior(i,nn,nb) = ensmean_chunk(i,nn,nb) +! remove mean from ensemble. + do nanal=1,nanals + anal_chunk(nanal,i,nn,nb) = anal_chunk(nanal,i,nn,nb)-ensmean_chunk(i,nn,nb) + anal_chunk_prior(nanal,i,nn,nb)=anal_chunk(nanal,i,nn,nb) + end do end do end do + enddo ! loop over nbackgrounds +!$omp end parallel do + +else if (ldo_enscalc_option ==1 ) then !to calculate ensemble mean +!$omp parallel do schedule(dynamic,1) private(nn,i,n,nb) +do nb=1,nbackgrounds + do nn=1,ncdim + do i=1,numptsperproc(nproc+1) + ensmean_chunk(i,nn,nb) = sum(anal_chunk(:,i,nn,nb))/float(nanals) + ensmean_chunk_prior(i,nn,nb) = ensmean_chunk(i,nn,nb) + do nanal=2,nanals + anal_chunk(nanal,i,nn,nb) = anal_chunk(nanal,i,nn,nb)-ensmean_chunk(i,nn,nb) + anal_chunk_prior(nanal,i,nn,nb)=anal_chunk(nanal,i,nn,nb) + end do + anal_chunk(1,i,nn,nb) =ensmean_chunk(i,nn,nb)- anal_chunk(1,i,nn,nb) + anal_chunk_prior(1,i,nn,nb)=0.0_r_kind + end do + end do +end do +!$omp end parallel do + +else if (ldo_enscalc_option ==2 ) then !to recentter , the first member is the control run +!$omp parallel do schedule(dynamic,1) private(nn,i,n,nb) +do nb=1,nbackgrounds + do nn=1,ncdim + do i=1,numptsperproc(nproc+1) + ensmean_chunk(i,nn,nb) = sum(anal_chunk(2:nanals,i,nn,nb))/float(nanals-1) + ensmean_chunk_prior(i,nn,nb) = ensmean_chunk(i,nn,nb) +! remove mean from ensemble. + do nanal=2,nanals + anal_chunk(nanal,i,nn,nb) = (anal_chunk(1,i,nn,nb)-ensmean_chunk(i,nn,nb)) + end do + do nanal=2,nanals + anal_chunk_prior(nanal,i,nn,nb)=0.0_r_kind + enddo + anal_chunk(1,i,nn,nb) = ensmean_chunk(i,nn,nb)-anal_chunk(1,i,nn,nb) + anal_chunk_prior(1,i,nn,nb)=0.0_r_kind + end do end do end do !$omp end parallel do +else +write(6,*)'this ldo_enscalc_option = ',ldo_enscalc_option,' is not available now' + +endif end subroutine scatter_chunks diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 36b0c9c207..667cdc2f01 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -215,6 +215,9 @@ module params ! next two are no longer used, instead they are inferred from anavinfo logical,public :: massbal_adjust = .false. +integer (i_kind), public :: ldo_enscalc_option = 0 +! true, the program will be used to calculate ensemble mean (=1) or recenter(=2) + integer(i_kind),public :: nvars = -1 ! sort obs in LETKF in order of decreasing DFS @@ -236,7 +239,7 @@ module params logical,public :: netcdf_diag = .false. ! use fv3 cubed-sphere tiled restart files -logical,public :: fv3_native = .false. +logical,public :: fv3_native = .true. character(len=500),public :: fv3fixpath = ' ' integer(i_kind),public :: ntiles=6 integer(i_kind),public :: nx_res=0,ny_res=0 @@ -275,7 +278,7 @@ module params paoverpb_thresh,latbound,delat,pseudo_rh,numiter,biasvar,& lupd_satbiasc,cliptracers,simple_partition,adp_anglebc,angord,& newpc4pred,nmmb,nhr_anal,nhr_state, fhr_assim,nbackgrounds,nstatefields, & - save_inflation,nobsl_max,lobsdiag_forenkf,netcdf_diag,forecast_impact,& + save_inflation,nobsl_max,lobsdiag_forenkf,ldo_enscalc_option,netcdf_diag,forecast_impact,& letkf_flag,massbal_adjust,use_edges,emiss_bc,iseed_perturbed_obs,npefiles,& getkf,getkf_inflation,denkf,modelspace_vloc,dfs_sort,write_spread_diag,& covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff,letkf_bruteforce_search,& @@ -658,6 +661,10 @@ subroutine read_namelist() print *, 'must specify nx_res,ny_res and fv3fixpath when fv3_native is true' call stop2(19) endif + if(ldo_enscalc_option .ne.0 ) then + l_pres_add_saved=.false. + endif + if (letkf_flag .and. univaroz) then print *,'univaroz is not supported in LETKF!' call stop2(19)