diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index ddd9ded945..b8b62e3b72 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -53,7 +53,7 @@ module controlvec use params, only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, & nanals, pseudo_rh, use_qsatensmean, nlons, nlats,& nanals_per_iotask, ntasks_io, nanal1, nanal2, & - fgsfcfileprefixes, paranc, write_fv3_incr + fgsfcfileprefixes, paranc, write_fv3_incr, write_ensmean use kinds, only: r_kind, i_kind, r_double, r_single use mpeu_util, only: gettablesize, gettable, getindex use constants, only: max_varname_length @@ -291,13 +291,14 @@ subroutine write_control(no_inflate_flag) real(r_double) :: t1,t2 integer(i_kind) :: nb, nvar, ne integer(i_kind) :: q_ind, ierr -real(r_single), allocatable, dimension(:,:) :: grdin_mean, grdin_mean_tmp +real(r_single), allocatable, dimension(:,:) :: grdin_mean_tmp +real(r_single), allocatable, dimension(:,:,:,:) :: grdin_mean if (nproc <= ntasks_io-1) then allocate(grdin_mean_tmp(npts,ncdim)) if (nproc == 0) then - allocate(grdin_mean(npts,ncdim)) + allocate(grdin_mean(npts,ncdim,nbackgrounds,1)) grdin_mean = 0_r_single t1 = mpi_wtime() endif @@ -311,27 +312,24 @@ subroutine write_control(no_inflate_flag) do ne=1,nanals_per_iotask call mpi_reduce(grdin(:,:,nb,ne), grdin_mean_tmp, npts*ncdim, mpi_real4,& mpi_sum,0,mpi_comm_io,ierr) - if (nproc == 0) grdin_mean = grdin_mean + grdin_mean_tmp + if (nproc == 0) grdin_mean(:,:,nb,1) = grdin_mean(:,:,nb,1) + grdin_mean_tmp enddo ! print out ens mean increment info if (nproc == 0) then - grdin_mean = grdin_mean/real(nanals) + grdin_mean(:,:,nb,1) = grdin_mean(:,:,nb,1)/real(nanals) do nvar=1,nc3d write(6,100) trim(cvars3d(nvar)), & - minval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar))), & - maxval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar))) + minval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar),nb,1)), & + maxval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar),nb,1)) enddo do nvar=1,nc2d write(6,100) trim(cvars2d(nvar)), & - minval(grdin_mean(:,clevels(nc3d) + nvar)), & - maxval(grdin_mean(:,clevels(nc3d) + nvar)) + minval(grdin_mean(:,clevels(nc3d) + nvar,nb,1)), & + maxval(grdin_mean(:,clevels(nc3d) + nvar,nb,1)) enddo endif enddo 100 format('ens. mean anal. increment min/max ',a,2x,g19.12,2x,g19.12) - if (nproc == 0) then - deallocate(grdin_mean) - endif deallocate(grdin_mean_tmp) q_ind = getindex(cvars3d, 'q') @@ -353,6 +351,14 @@ subroutine write_control(no_inflate_flag) enddo enddo endif + if (nproc == 0 .and. write_ensmean) then + ! write_ensmean implies use_qsatensmean + do nb=1,nbackgrounds + ! re-scale normalized spfh with sat. sphf of ensmean first guess + grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1) = & + grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1)*qsatmean(:,:,nb) + enddo + endif end if if (.not. paranc) then if (write_fv3_incr) then @@ -361,6 +367,15 @@ subroutine write_control(no_inflate_flag) call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) end if if (nproc == 0) then + if (write_ensmean) then + ! also write out ens mean on root task. + if (write_fv3_incr) then + call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) + else + call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) + end if + endif + deallocate(grdin_mean) t2 = mpi_wtime() print *,'time in write_control on root',t2-t1,'secs' endif @@ -375,6 +390,15 @@ subroutine write_control(no_inflate_flag) call writegriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) end if if (nproc == 0) then + ! also write out ens mean on root task + if (write_ensmean) then ! FIXME use parallel IO to write ensmean + if (write_fv3_incr) then + call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) + else + call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) + end if + endif + deallocate(grdin_mean) t2 = mpi_wtime() print *,'time in write_control on root',t2-t1,'secs' endif diff --git a/src/enkf/enkf_obsmod.f90 b/src/enkf/enkf_obsmod.f90 index 72e196d61e..d9ab492e73 100644 --- a/src/enkf/enkf_obsmod.f90 +++ b/src/enkf/enkf_obsmod.f90 @@ -105,7 +105,7 @@ module enkf_obsmod use kinds, only : r_kind, r_double, i_kind, r_single use constants, only: zero, one, deg2rad, rad2deg, rd, cp, pi use params, only: & - datestring,datapath,sprd_tol,nanals,saterrfact, & + letkf_flag,nobsl_max,datestring,datapath,sprd_tol,nanals,saterrfact, & lnsigcutoffnh, lnsigcutoffsh, lnsigcutofftr, corrlengthnh,& corrlengthtr, corrlengthsh, obtimelnh, obtimeltr, obtimelsh,& lnsigcutoffsatnh, lnsigcutoffsatsh, lnsigcutoffsattr,& @@ -204,6 +204,13 @@ subroutine readobs() if (nproc == 0) then print *,'max time in mpireadobs = ',tdiffmax print *,'total number of obs ',nobstot + print *,'min/max obtime ',minval(obtime),maxval(obtime) +endif +! if nobsl_max set for LETKF, and the total number of obs < nobsl_max, +! reset nobsl_max to -1 +if (letkf_flag .and. nobsl_max > 0 .and. nobstot < nobsl_max) then + if (nproc == 0) print *,'resetting nobsl_max to -1' + nobsl_max=-1 endif allocate(obfit_prior(nobstot)) ! screen out some obs by setting ob error to a very large number diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 58ceff78a6..8a6bd12507 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -333,14 +333,19 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & if (iope==0) then do k=1,nlevs - krev = nlevs-k+1 - ! layer pressure from phillips vertical interolation + ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90) + if (prse_ind > 0) then + ug(:) = pressi(:,k) + call copytogrdin(ug,pslg(:,k)) + ! Jacobian for gps in pressure is saved in different units in GSI; need to + ! multiply pressure by 0.1 + grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) + endif + ! layer pressure from phillips vertical interolation (used for qsat + ! calculation) ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/& (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr call copytogrdin(ug,pslg(:,k)) - ! Jacobian for gps in pressure is saved in different units in GSI; need to - ! multiply pressure by 0.1 - if (prse_ind > 0) grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) end do if (pseudo_rh) then call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs) @@ -952,15 +957,19 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ! compute saturation q. do k=1,nlevs - ! layer pressure from phillips vertical interolation + ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90) + if (prse_ind > 0) then + ug(:) = pressi(:,k) + call copytogrdin(ug,pslg(:,k)) + ! Jacobian for gps in pressure is saved in different units in GSI; need to + ! multiply pressure by 0.1 + grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) + endif + ! layer pressure from phillips vertical interolation (used for qsat + ! calculation) ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/& (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr - call copytogrdin(ug,pslg(:,k)) - ! Jacobian for gps in pressure is saved in different units in GSI; need to - ! multiply pressure by 0.1 - if (prse_ind > 0) grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) - end do if (pseudo_rh) then call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs) @@ -1115,7 +1124,9 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ ! need to distribute grdin to all PEs in this subcommunicator ! bring all the subdomains back to the main PE call mpi_barrier(iocomms(mem_pe(nproc)), iret) - call mpi_bcast(grdin,npts*ndim*nbackgrounds, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) + do nb=1,nbackgrounds + call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) + enddo ! loop through times and do the read ne = 1 @@ -1838,7 +1849,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n write_attribute, quantize_data, has_var, has_attr use constants, only: grav use params, only: nbackgrounds,anlfileprefixes,fgfileprefixes,reducedgrid,& - nccompress + nccompress,write_ensmean implicit none integer, intent(in) :: nanal1,nanal2 @@ -1906,12 +1917,17 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n write(charnanal,'(i3.3)') nanal backgroundloop: do nb=1,nbackgrounds - if(no_inflate_flag) then - filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"nimem"//charnanal + if (nanal == 0 .and. write_ensmean) then + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean" + filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"ensmean" else - filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"mem"//charnanal - end if - filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal + if(no_inflate_flag) then + filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"nimem"//charnanal + else + filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"mem"//charnanal + end if + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal + endif if (use_gfs_nemsio) then clip = tiny(vg(1)) @@ -3289,7 +3305,7 @@ end subroutine writegriddata subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& - datestring,nhr_anal + datestring,nhr_anal,write_ensmean use constants, only: grav use mpi use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -3354,12 +3370,17 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, write(charnanal,'(i3.3)') nanal backgroundloop: do nb=1,nbackgrounds - if(no_inflate_flag) then - filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal + if (nanal == 0 .and. write_ensmean) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"ensmean" + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean" else - filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal - end if - filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal + if(no_inflate_flag) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal + else + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal + end if + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal + endif ! create the output netCDF increment file call nccheck_incr(nf90_create(path=trim(filenameout), cmode=nf90_netcdf4, ncid=ncid_out)) @@ -3773,7 +3794,9 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ! need to distribute grdin to all PEs in this subcommunicator ! bring all the subdomains back to the main PE call mpi_barrier(iocomms(mem_pe(nproc)), iret) - call mpi_bcast(grdin,npts*ndim*nbackgrounds, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) + do nb=1,nbackgrounds + call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) + enddo ! loop through times and do the read ne = 1 diff --git a/src/enkf/letkf.f90 b/src/enkf/letkf.f90 index 0f04739392..9b74cecd75 100644 --- a/src/enkf/letkf.f90 +++ b/src/enkf/letkf.f90 @@ -107,13 +107,13 @@ module letkf numobspersat, biaspreds, corrlengthsq,& probgrosserr, prpgerr, obtype, obpress,& lnsigl, anal_ob, anal_ob_modens, obloclat, obloclon, stattype -use constants, only: pi, one, zero, rad2deg, deg2rad +use constants, only: pi, one, zero, rad2deg, deg2rad, rearth use params, only: sprd_tol, datapath, nanals, iseed_perturbed_obs,& iassim_order,sortinc,deterministic,nlevs,& zhuberleft,zhuberright,varqc,lupd_satbiasc,huber,letkf_novlocal,& lupd_obspace_serial,corrlengthnh,corrlengthtr,corrlengthsh,& getkf,getkf_inflation,denkf,nbackgrounds,nobsl_max,& - neigv,vlocal_evecs,dfs_sort + neigv,vlocal_evecs,dfs_sort,mincorrlength_fact use gridinfo, only: nlevs_pres,lonsgrd,latsgrd,logp,npts,gridloc use kdtree2_module, only: kdtree2, kdtree2_create, kdtree2_destroy, & kdtree2_result, kdtree2_n_nearest, kdtree2_r_nearest @@ -135,13 +135,16 @@ subroutine letkf_update() integer(i_kind) nob,nf,nanal,nens,& i,nlev,nrej,npt,nn,nnmax,ierr integer(i_kind) nobsl, ngrd1, nobsl2, nthreads, nb, & - nobslocal_min,nobslocal_max, & - nobslocal_minall,nobslocal_maxall + nobslocal_mean,nobslocal_min,nobslocal_max, & + nobslocal_meanall,nobslocal_minall,nobslocal_maxall +real(r_single) robslocal_mean,robslocal_min,robslocal_max,re, & + robslocal_meanall,robslocal_minall,robslocal_maxall,& + coslatslocal_meanall, coslatslocal_mean, coslat integer(i_kind),allocatable,dimension(:) :: oindex -real(r_single) :: deglat, dist, corrsq, oberrfact, trpa, trpa_raw +real(r_single) :: deglat, dist, corrsq, trpa, trpa_raw, maxdfs real(r_double) :: t1,t2,t3,t4,t5,tbegin,tend,tmin,tmax,tmean real(r_kind) r_nanals,r_nanalsm1 -real(r_kind) normdepart, pnge, width +real(r_kind) normdepart, pnge, width, mincorrlength_factsq real(r_kind),dimension(nobstot):: oberrvaruse real(r_kind) vdist real(r_kind) corrlength @@ -152,7 +155,7 @@ subroutine letkf_update() real(r_single),allocatable,dimension(:,:,:) :: ens_tmp real(r_single),allocatable,dimension(:,:) :: wts_ensperts,pa real(r_single),allocatable,dimension(:) :: dfs,wts_ensmean -real(r_kind),allocatable,dimension(:) :: rdiag,rloc +real(r_kind),allocatable,dimension(:) :: rdiag,rloc,robs_local,coslats_local real(r_single),allocatable,dimension(:) :: dep ! kdtree stuff type(kdtree2_result),dimension(:),allocatable :: sresults @@ -160,6 +163,7 @@ subroutine letkf_update() real(r_kind) eps eps = epsilon(0.0_r_single) ! real(4) machine precision +re = rearth/1.e3_r_single !$omp parallel nthreads = omp_get_num_threads() @@ -170,6 +174,7 @@ subroutine letkf_update() ! define a few frequently used parameters r_nanals=one/float(nanals) r_nanalsm1=one/float(nanals-1) +mincorrlength_factsq = mincorrlength_fact**2 kdobs=associated(kdtree_obs2) if (.not. kdobs .and. nproc .eq. 0) then @@ -258,17 +263,34 @@ subroutine letkf_update() t4 = zero t5 = zero tbegin = mpi_wtime() + nobslocal_max = -999 nobslocal_min = nobstot +nobslocal_mean = 0 +allocate(robs_local(npts_max)) +robs_local = 0 +if (nobsl_max > 0) then + allocate(coslats_local(npts_max)) + coslats_local = 0 +endif ! Update ensemble on model grid. ! Loop for each horizontal grid points on this task. -!$omp parallel do schedule(dynamic) private(npt,nob,nobsl, & -!$omp nobsl2,oberrfact,ngrd1,corrlength,ens_tmp, & -!$omp nf,vdist,obens,indxassim,indxob, & +!$omp parallel do schedule(dynamic) default(none) private(npt,nob,nobsl, & +!$omp nobsl2,ngrd1,corrlength,ens_tmp,coslat, & +!$omp nf,vdist,obens,indxassim,indxob,maxdfs, & !$omp nn,hxens,wts_ensmean,dfs,rdiag,dep,rloc,i, & -!$omp oindex,deglat,dist,corrsq,nb,sresults, & -!$omp wts_ensperts,pa,trpa,trpa_raw) & +!$omp oindex,deglat,dist,corrsq,nb,nlev,nanal,sresults, & +!$omp wts_ensperts,pa,trpa,trpa_raw) shared(anal_ob, & +!$omp anal_ob_modens,anal_chunk,obsprd_post,obsprd_prior, & +!$omp oberrvar,oberrvaruse,nobsl_max,grdloc_chunk, & +!$omp obloc,corrlengthnh,corrlengthsh,corrlengthtr,& +!$omp vlocal_evecs,vlocal,oblnp,lnp_chunk,lnsigl,corrlengthsq,& +!$omp getkf,denkf,getkf_inflation,ensmean_chunk,ob,ensmean_ob, & +!$omp nproc,numptsperproc,nnmax,r_nanalsm1,kdtree_obs2,kdobs, & +!$omp mincorrlength_factsq,robs_local,coslats_local, & +!$omp lupd_obspace_serial,eps,dfs_sort,nanals,index_pres,& +!$omp neigv,nlevs,lonsgrd,latsgrd,nobstot,nens,ncdim,nbackgrounds,indxproc,rad2deg) & !$omp reduction(+:t1,t2,t3,t4,t5) & !$omp reduction(max:nobslocal_max) & !$omp reduction(min:nobslocal_min) @@ -279,6 +301,7 @@ subroutine letkf_update() ! find obs close to this grid point (using kdtree) ngrd1=indxproc(nproc+1,npt) deglat = latsgrd(ngrd1)*rad2deg + coslat = cos(latsgrd(ngrd1)) corrlength=latval(deglat,corrlengthnh,corrlengthtr,corrlengthsh) corrsq = corrlength**2 allocate(sresults(nobstot)) @@ -304,29 +327,26 @@ subroutine letkf_update() allocate(indxob(nobstot)) ! calculate integrated 1-DFS for each ob in local volume nobsl = 0 + maxdfs = -9.9e31 do nob=1,nobstot rloc(nob) = sum((obloc(:,nob)-grdloc_chunk(:,npt))**2,1) - dist = sqrt(rloc(nob)/corrlengthsq(nob)) + dist = sqrt(rloc(nob)/corrsq) if (dist < 1.0 - eps .and. & oberrvaruse(nob) < 1.e10_r_single) then nobsl = nobsl + 1 indxob(nobsl) = nob - oberrfact = taper(dist) - if (lupd_obspace_serial) then - ! use updated ensemble in ob space to estimate DFS - !dfs(nobsl) = obsprd_post(nob)/obsprd_prior(nob) - ! weight by distance to analysis point - dfs(nobsl) = oberrfact*obsprd_post(nob)/obsprd_prior(nob) - else - ! estimate DFS assuming each ob assimilated independently, one - ! at a time. - ! 1-DFS = HP_aH^T/HP_bH^T = R/(HP_bH^T + R) - dfs(nobsl) = (oberrvaruse(nob)/oberrfact)/((oberrvar(nob)/oberrfact)+obsprd_prior(nob)) - endif + ! use updated ensemble in ob space to compute DFS + ! DFS = Tr(R**-1*HPaHT) = dy_a/dy_o see eqn 4 in Liu et al 2009 + ! https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.511 + dfs(nobsl) = obsprd_post(nob)/oberrvaruse(nob) + ! use spread reduction instead. + !dfs(nobsl) = obsprd_post(nob)/obsprd_prior(nob) + !if (dfs(nobsl) > maxdfs) maxdfs = dfs(nobsl) endif enddo - ! sort on 1-DFS + ! sort on max(DFS)-DFS allocate(indxassim(nobsl)) + dfs = maxdfs-dfs call quicksort(nobsl,dfs(1:nobsl),indxassim) nobsl2 = min(nobsl_max,nobsl) do nob=1,nobsl2 @@ -345,7 +365,6 @@ subroutine letkf_update() else ! brute force search call find_localobs(grdloc_chunk(:,npt),obloc,corrsq,nobstot,nobsl_max,sresults,nobsl) - nobsl_max = nobsl endif !if (nproc == 0 .and. npt == 1) then ! do nob=1,nobsl @@ -372,6 +391,12 @@ subroutine letkf_update() if (allocated(ens_tmp)) deallocate(ens_tmp) cycle grdloop endif + if (nobsl_max > 0) then + robs_local(npt) = sqrt(sresults(nobsl)%dis) + coslats_local(npt) = coslat + else + robs_local(npt) = nobsl + endif ! Loop through vertical levels (nnmax=1 if no vertical localization) verloop: do nn=1,nnmax @@ -390,7 +415,20 @@ subroutine letkf_update() else vdist = zero endif - dist = sqrt(sresults(nob)%dis/corrlengthsq(nf)+vdist*vdist) + if (nobsl_max > 0 .and. corrlength < 0) then + ! if corrlength<0, set R localization scale to be max distance to find nobsl_max obs + ! (unless max distance is > abs(corrlength) or < abs(corrlength)/10) + if (sresults(nobsl)%dis > corrsq) then + dist = sqrt(sresults(nob)%dis/corrsq+vdist*vdist) + else if (sresults(nobsl)%dis < corrsq*mincorrlength_factsq) then + dist = sqrt(sresults(nob)%dis/(corrsq/mincorrlength_factsq)+vdist*vdist) + else + dist = sqrt(sresults(nob)%dis/sresults(nobsl)%dis+vdist*vdist) + endif + else + ! set R localization scale to specificed distance + dist = sqrt(sresults(nob)%dis/corrsq+vdist*vdist) + endif if (dist >= one) cycle rloc(nobsl2)=taper(dist) oindex(nobsl2)=nf @@ -490,7 +528,9 @@ subroutine letkf_update() ! make sure posterior perturbations still have zero mean. ! (roundoff errors can accumulate) -!$omp parallel do schedule(dynamic) private(npt,nb,i) +!$omp parallel do schedule(dynamic) default(none) private(npt,nb,i) & +!$omp shared(anal_chunk,r_nanals,nanals,& +!$omp npts_max,nbackgrounds,ncdim) do npt=1,npts_max do nb=1,nbackgrounds do i=1,ncdim @@ -529,10 +569,35 @@ subroutine letkf_update() call mpi_reduce(t5,tmin,1,mpi_real8,mpi_min,0,mpi_comm_world,ierr) call mpi_reduce(t5,tmax,1,mpi_real8,mpi_max,0,mpi_comm_world,ierr) if (nproc .eq. 0) print *,',min/max/mean t5 = ',tmin,tmax,tmean + +if (nobsl_max > 0) then + ! compute and print min/max/mean search radius to find nobsl_max + robslocal_mean = sum(robs_local*coslats_local)/numptsperproc(nproc+1) + coslatslocal_mean = sum(coslats_local)/numptsperproc(nproc+1) + robslocal_min = minval(robs_local(1:numptsperproc(nproc+1))) + robslocal_max = maxval(robs_local(1:numptsperproc(nproc+1))) + call mpi_reduce(robslocal_max,robslocal_maxall,1,mpi_real4,mpi_max,0,mpi_comm_world,ierr) + call mpi_reduce(robslocal_min,robslocal_minall,1,mpi_real4,mpi_min,0,mpi_comm_world,ierr) + call mpi_reduce(robslocal_mean,robslocal_meanall,1,mpi_real4,mpi_sum,0,mpi_comm_world,ierr) + call mpi_reduce(coslatslocal_mean,coslatslocal_meanall,1,mpi_real4,mpi_sum,0,mpi_comm_world,ierr) + if (nproc == 0) print *,'min/max/mean distance searched for local obs',re*robslocal_minall,re*robslocal_maxall,re*robslocal_meanall/coslatslocal_meanall + deallocate(coslats_local) +else + ! compute and print min/max/mean number of obs found within search radius + nobslocal_mean = nint(sum(robs_local)/numptsperproc(nproc+1)) + nobslocal_min = minval(robs_local(1:numptsperproc(nproc+1))) + nobslocal_max = maxval(robs_local(1:numptsperproc(nproc+1))) + call mpi_reduce(nobslocal_max,nobslocal_maxall,1,mpi_integer,mpi_max,0,mpi_comm_world,ierr) + call mpi_reduce(nobslocal_min,nobslocal_minall,1,mpi_integer,mpi_min,0,mpi_comm_world,ierr) + call mpi_reduce(nobslocal_mean,nobslocal_meanall,1,mpi_integer,mpi_sum,0,mpi_comm_world,ierr) + if (nproc == 0) print *,'min/max/mean number of obs in local volume',nobslocal_minall,nobslocal_maxall,nint(nobslocal_meanall/float(numproc)) +endif call mpi_reduce(nobslocal_max,nobslocal_maxall,1,mpi_integer,mpi_max,0,mpi_comm_world,ierr) call mpi_reduce(nobslocal_min,nobslocal_minall,1,mpi_integer,mpi_max,0,mpi_comm_world,ierr) if (nproc == 0) print *,'min/max number of obs in local volume',nobslocal_minall,nobslocal_maxall + if (nrej > 0 .and. nproc == 0) print *, nrej,' obs rejected by varqc' +deallocate(robs_local) if (allocated(ens_tmp)) deallocate(ens_tmp) diff --git a/src/enkf/loadbal.f90 b/src/enkf/loadbal.f90 index 595ef20277..ff99bc0736 100644 --- a/src/enkf/loadbal.f90 +++ b/src/enkf/loadbal.f90 @@ -383,7 +383,7 @@ subroutine scatter_chunks ! allocate array to hold pieces of state vector on each proc. allocate(anal_chunk(nanals,npts_max,ncdim,nbackgrounds)) -if (nproc == 0) print *,'anal_chunk size = ',size(anal_chunk) +if (nproc == 0) print *,'anal_chunk size = ',size(anal_chunk,kind=8) allocate(anal_chunk_prior(nanals,npts_max,ncdim,nbackgrounds)) allocate(ensmean_chunk(npts_max,ncdim,nbackgrounds)) diff --git a/src/enkf/mpi_readobs.f90 b/src/enkf/mpi_readobs.f90 index bca5d7c715..e48f5a1804 100644 --- a/src/enkf/mpi_readobs.f90 +++ b/src/enkf/mpi_readobs.f90 @@ -250,7 +250,18 @@ subroutine mpi_getobs(obspath, datestring, nobs_conv, nobs_oz, nobs_sat, nobs_to ! exchange obs prior ensemble members across all tasks to fully populate shared ! memory array pointer on each node. if (nproc_shm == 0) then - call mpi_allreduce(mpi_in_place,anal_ob,nanals*nobs_tot,mpi_real4,mpi_sum,mpi_comm_shmemroot,ierr) + if (real(nanals)*real(nobs_tot) < 2**32/2. - 1) then + call mpi_allreduce(mpi_in_place,anal_ob,nanals*nobs_tot,mpi_real4,mpi_sum,mpi_comm_shmemroot,ierr) + else + ! count won't fit in 32-bit integer and mpi_allreduce doesn't handle + ! 64 bit counts. Split up into smaller chunks. + mem_ob = 0. + do na=1,nanals + mem_ob(:) = anal_ob(na,:) + call mpi_allreduce(mpi_in_place,mem_ob,nobs_tot,mpi_real4,mpi_sum,mpi_comm_shmemroot,ierr) + anal_ob(na,:) = mem_ob(:) + enddo + endif !print *,nproc,'min/max anal_ob',minval(anal_ob),maxval(anal_ob) if (neigv > 0) then mem_ob_modens = 0. diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 6a3a977f56..e46bfa195b 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -115,6 +115,9 @@ module params real(r_single),public :: covinflatemax,covinflatemin,smoothparm,biasvar real(r_single),public :: corrlengthnh,corrlengthtr,corrlengthsh real(r_single),public :: obtimelnh,obtimeltr,obtimelsh +! factor for minimum allowed horiz cov length scale +! to apply for LETKF when corrlengthnh,tr,sh < 0 and nobsl_max > 0 +real(r_single),public :: mincorrlength_fact = 0.1 real(r_single),public :: zhuberleft,zhuberright real(r_single),public :: lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,& lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,& @@ -129,7 +132,10 @@ module params real(r_single),public :: tar_minlat,tar_maxlat,tar_minlon,tar_maxlon real(r_single),public :: covl_minfact, covl_efold -real(r_single),public :: covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff +real(r_single),public :: covinflatenh=0 +real(r_single),public :: covinflatetr=0 +real(r_single),public :: covinflatesh=0 +real(r_single),public :: lnsigcovinfcutoff ! if npefiles=0, diag files are read (concatenated pe* files written by gsi) ! if npefiles>0, npefiles+1 pe* files read directly ! the pe* files are assumed to be located in /gsitmp_mem### @@ -245,10 +251,12 @@ module params ! for writing increments logical,public :: write_fv3_incr = .false. character(len=12),dimension(10),public :: incvars_to_zero='NONE' !just picking 10 arbitrarily +! write ensemble mean analysis (or analysis increment) +logical,public :: write_ensmean = .false. namelist /nam_enkf/datestring,datapath,iassim_order,nvars,& covinflatemax,covinflatemin,deterministic,sortinc,& - corrlengthnh,corrlengthtr,corrlengthsh,& + mincorrlength_fact,corrlengthnh,corrlengthtr,corrlengthsh,& varqc,huber,nlons,nlats,smoothparm,use_qsatensmean,& readin_localization, zhuberleft,zhuberright,& obtimelnh,obtimeltr,obtimelsh,reducedgrid,& @@ -273,7 +281,7 @@ module params eft,wmoist,adrate,andataname,& gdatehr,datehr,& tar_minlat,tar_maxlat,tar_minlon,tar_maxlon,tar_minlev,tar_maxlev,& - fv3_native, paranc, nccompress, write_fv3_incr,incvars_to_zero, & + fv3_native, paranc, nccompress, write_fv3_incr,incvars_to_zero,write_ensmean, & corrlengthrdrnh,corrlengthrdrsh,corrlengthrdrtr,& lnsigcutoffrdrnh,lnsigcutoffrdrsh,lnsigcutoffrdrtr,& l_use_enkf_directZDA @@ -664,6 +672,10 @@ subroutine read_namelist() letkf_flag) then print *,'warning: no time localization in LETKF!' endif + if ((write_ensmean .and. pseudo_rh) .and. .not. use_qsatensmean) then + print *,'write_ensmean=T requires use_qsatensmean=T when pseudo_rh=T' + call stop2(19) + endif print *, trim(adjustl(datapath)) diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index 8c0f7ea958..e1977298a6 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -335,8 +335,8 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) errorlimit2=errorlimit2_obs - if (obtype == 'gps' ) then - if (GPS_Type(i)==1) errorlimit2=errorlimit2_bnd + if (obtype == 'gps') then + if (GPS_Type(i)==1) errorlimit2=errorlimit2_bnd endif ! for q, normalize by qsatges @@ -664,7 +664,7 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & do i = 1, nobs nobdiag = nobdiag + 1 ! special handling for error limits for GPS bend angle - if (obtype == 'gps' ) then + if (obtype == 'gps') then if (GPS_Type(i)==1) errorlimit2=errorlimit2_bnd endif diff --git a/src/enkf/readozobs.f90 b/src/enkf/readozobs.f90 index efd79b1c8c..efebb00855 100644 --- a/src/enkf/readozobs.f90 +++ b/src/enkf/readozobs.f90 @@ -926,7 +926,7 @@ end subroutine write_ozobs_data_bin subroutine write_ozobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & x_fit, x_sprd, x_used, id, gesid) use netcdf, only: nf90_inq_dimid, nf90_open, nf90_close, NF90_NETCDF4, & - nf90_inquire_dimension, NF90_WRITE, nf90_create, nf90_def_dim + nf90_inquire_dimension, NF90_WRITE, NF90_NOWRITE, nf90_create, nf90_def_dim use ncdw_climsg, only: nclayer_check use constants, only: r_missing @@ -972,7 +972,7 @@ subroutine write_ozobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & inquire(file=obsfile,exist=fexist) if (.not. fexist) cycle peloop - call nclayer_check(nf90_open(obsfile, NF90_WRITE, iunit)) + call nclayer_check(nf90_open(obsfile, NF90_NOWRITE, iunit)) call nclayer_check(nf90_inq_dimid(iunit, "nobs", nobsid)) call nclayer_check(nf90_inquire_dimension(iunit, nobsid, len = nobs)) call nclayer_check(nf90_close(iunit)) diff --git a/src/enkf/readsatobs.f90 b/src/enkf/readsatobs.f90 index 44923d2c38..df0534fc90 100644 --- a/src/enkf/readsatobs.f90 +++ b/src/enkf/readsatobs.f90 @@ -1260,7 +1260,7 @@ end subroutine write_satobs_data_bin subroutine write_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & x_fit, x_sprd, x_used, id, gesid) use netcdf, only: nf90_inq_dimid, nf90_open, nf90_close, NF90_NETCDF4, & - nf90_inquire_dimension, NF90_WRITE, nf90_create, nf90_def_dim + nf90_inquire_dimension, NF90_WRITE, NF90_NOWRITE, nf90_create, nf90_def_dim use ncdw_climsg, only: nclayer_check use radinfo, only: iuse_rad,nusis,jpch_rad @@ -1326,7 +1326,7 @@ subroutine write_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & if (.not. fexist) cycle peloop - call nclayer_check(nf90_open(obsfile, NF90_WRITE, iunit)) + call nclayer_check(nf90_open(obsfile, NF90_NOWRITE, iunit)) call nclayer_check(nf90_inq_dimid(iunit, "nobs", nobsid)) call nclayer_check(nf90_inquire_dimension(iunit, nobsid, len = nobs)) call nclayer_check(nf90_close(iunit)) diff --git a/src/gsi/gsi_4dvar.f90 b/src/gsi/gsi_4dvar.f90 index d955bd57d1..e4038dc9ba 100644 --- a/src/gsi/gsi_4dvar.f90 +++ b/src/gsi/gsi_4dvar.f90 @@ -435,8 +435,9 @@ subroutine time_4dvar(idate,step4d) integer(i_kind),intent(in ) :: idate ! Date (yyyymmddhh) real(r_kind) ,intent( out) :: step4d ! Time since start of 4D-Var window (hours) -integer(i_kind) iyr,imo,idy,ihr,nmin_obs,nhrobs,nhrbgn,nhroff +integer(i_kind) iyr,imo,idy,ihr,nmin_obs integer(i_kind),dimension(5) :: idate5 +real(r_kind) nhroff,nhrbgn,nhrobs ihr=idate iyr=ihr/1000000 @@ -457,7 +458,7 @@ subroutine time_4dvar(idate,step4d) end if nhrobs=nmin_obs*r60inv -nhrbgn=NINT(real(iwinbgn,r_kind)*r60inv) +nhrbgn=real(iwinbgn,r_kind)*r60inv nhroff=nhrobs-nhrbgn step4d=real(nhroff,r_kind) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 1c8b73841f..0c39ffb04d 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -145,7 +145,7 @@ module gsimod n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test,oz_univ_static,& regional_ensemble_option,fv3sar_ensemble_opt,merge_two_grid_ensperts, & full_ensemble,pseudo_hybens,pwgtflg,& - beta_s0,s_ens_h,s_ens_v,init_hybrid_ensemble_parameters,& + beta_s0,beta_e0,s_ens_h,s_ens_v,init_hybrid_ensemble_parameters,& readin_localization,write_ens_sprd,eqspace_ensgrid,grid_ratio_ens,& readin_beta,use_localization_grid,use_gfs_ens,q_hyb_ens,i_en_perts_io, & l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB @@ -1268,6 +1268,9 @@ module gsimod ! beta_s(:) = beta_s0 , vertically varying weights given to static B ; ! beta_e(:) = 1 - beta_s0 , vertically varying weights given ensemble derived covariance. ! If (readin_beta) then beta_s and beta_e are read from a file and beta_s0 is not used. +! beta_e0 - default weight given to ensemble background error covariance +! (if .not. readin_beta). if beta_e0<0, then it is set to +! 1.-beta_s0 (this is the default) ! s_ens_h - homogeneous isotropic horizontal ensemble localization scale (km) ! s_ens_v - vertical localization scale (grid units for now) ! s_ens_h, s_ens_v, and beta_s0 are tunable parameters. @@ -1312,7 +1315,7 @@ module gsimod ! namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,nlon_ens,nlat_ens,jcap_ens,& pseudo_hybens,merge_two_grid_ensperts,regional_ensemble_option,fv3sar_bg_opt,fv3sar_ensemble_opt,full_ensemble,pwgtflg,& - jcap_ens_test,beta_s0,s_ens_h,s_ens_v,readin_localization,eqspace_ensgrid,readin_beta,& + jcap_ens_test,beta_s0,beta_e0,s_ens_h,s_ens_v,readin_localization,eqspace_ensgrid,readin_beta,& grid_ratio_ens, & oz_univ_static,write_ens_sprd,use_localization_grid,use_gfs_ens, & i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index e0a6e3fd1f..23a21cd8e5 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -4001,7 +4001,7 @@ subroutine hybens_localization_setup use gfs_stratosphere, only: use_gfs_stratosphere,blend_rm use hybrid_ensemble_parameters, only: grd_ens,jcap_ens,n_ens,grd_loc,sp_loc,& nval_lenz_en,regional_ensemble_option - use hybrid_ensemble_parameters, only: readin_beta,beta_s,beta_e,beta_s0,sqrt_beta_s,sqrt_beta_e + use hybrid_ensemble_parameters, only: readin_beta,beta_s,beta_e,beta_s0,beta_e0,sqrt_beta_s,sqrt_beta_e use hybrid_ensemble_parameters, only: readin_localization,create_hybens_localization_parameters, & vvlocal,s_ens_h,s_ens_hv,s_ens_v,s_ens_vv use gsi_io, only: verbose @@ -4067,7 +4067,11 @@ subroutine hybens_localization_setup if ( .not. readin_beta ) then ! assign all levels to same value, sum = 1.0 beta_s = beta_s0 - beta_e = one - beta_s0 + if (beta_e0 < 0) then + beta_e = one - beta_s0 + else + beta_e = beta_e0 + endif endif if ( regional_ensemble_option == 2 .and. use_gfs_stratosphere .and. .not. readin_beta ) then diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90 index 9ea28c9ee8..61723bbc65 100644 --- a/src/gsi/hybrid_ensemble_parameters.f90 +++ b/src/gsi/hybrid_ensemble_parameters.f90 @@ -87,6 +87,9 @@ module hybrid_ensemble_parameters ! relative weight given to static background B when (readin_beta=.false.) ! when (readin_beta=.true.), the vertical weighting parameters are read from a file, ! instead of being defined based on beta_s0 namelist or default value. +! beta_e0 - default weight given to ensemble background error covariance +! (if .not. readin_beta). if beta_e0<0, then it is set to +! 1.-beta_s0 (this is the default) ! s_ens_h: horizontal localization correlation length (units of km), default = 2828.0 ! s_ens_v: vertical localization correlation length (grid units), default = 30.0 ! generate_ens: if .true., generate ensemble perturbations internally as random samples of background B. @@ -255,7 +258,7 @@ module hybrid_ensemble_parameters ! set passed variables to public public :: generate_ens,n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test,l_hyb_ens,& s_ens_h,oz_univ_static,vvlocal - public :: uv_hyb_ens,q_hyb_ens,s_ens_v,beta_s0,aniso_a_en,s_ens_hv,s_ens_vv + public :: uv_hyb_ens,q_hyb_ens,s_ens_v,beta_s0,beta_e0,aniso_a_en,s_ens_hv,s_ens_vv public :: readin_beta,beta_s,beta_e public :: readin_localization public :: eqspace_ensgrid,grid_ratio_ens @@ -309,7 +312,7 @@ module hybrid_ensemble_parameters logical ens_fast_read integer(i_kind) i_en_perts_io integer(i_kind) n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test - real(r_kind) beta_s0,s_ens_h,s_ens_v,grid_ratio_ens + real(r_kind) beta_s0,beta_e0,s_ens_h,s_ens_v,grid_ratio_ens type(sub2grid_info),save :: grd_ens,grd_loc,grd_sploc,grd_anl,grd_e1,grd_a1 type(spec_vars),save :: sp_ens,sp_loc type(egrid2agrid_parm),save :: p_e2a,p_sploc2ens @@ -401,6 +404,7 @@ subroutine init_hybrid_ensemble_parameters jcap_ens_test=0 nlon_ens=0 beta_s0=one + beta_e0=-one grid_ratio_ens=one s_ens_h = 2828._r_kind ! km (this was optimal value in ! Wang, X.,D. M. Barker, C. Snyder, and T. M. Hamill, 2008: A hybrid diff --git a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/fv3_interface.f90 b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/fv3_interface.f90 index 2eb31b836e..e4e4a466cd 100644 --- a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/fv3_interface.f90 +++ b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/fv3_interface.f90 @@ -133,6 +133,7 @@ subroutine fv3_calc_increment(mype) integer :: j, k ! loop indices within a variable integer :: ivar !! loop index over variables in input_vars & output_vars + real :: taper ! Formats for print statements: 100 format(A,': ',A) @@ -185,6 +186,9 @@ subroutine fv3_calc_increment(mype) an_grid%ilev(k) = real(k) an_grid%hyai(k) = real(k) an_grid%hybi(k) = real(k) + an_grid%ak(k) = meta_ncio%vcoord(k,1) + an_grid%bk(k) = meta_ncio%vcoord(k,2) + an_grid%ck(k) = 0 end do nzp1_init ! Deallocate entire grid. @@ -240,8 +244,26 @@ subroutine fv3_calc_increment(mype) ! Subtract and write an_grid%var3d = an_grid%var3d - fg_grid%var3d + if (mype == 0) print *,trim(input_vars(ivar)),minval(an_grid%var3d),maxval(an_grid%var3d) endif zero_or_read + ! taper humidity, microphysics increments in stratosphere + if (taper_strat .and. (trim(input_vars(ivar)) == 'spfh' .or. & + trim(input_vars(ivar)) == 'icmr' .or. & + trim(input_vars(ivar)) == 'clwmr')) then + if (mype == 0) print *,'k,ak,bk,taper,min/max increment for ',trim(input_vars(ivar)) + do k=1,an_grid%nz + taper = 1.0 + if (k < an_grid%nz/2 .and. (an_grid%ak(k) <= ak_bot .and. an_grid%ak(k) >= ak_top)) then + taper = (an_grid%ak(k) - ak_top)/(ak_bot - ak_top) + else if (an_grid%bk(k) .eq. 0. .and. an_grid%ak(k) < ak_top) then + taper = 0. + endif + an_grid%var3d(:,:,k) = an_grid%var3d(:,:,k)*taper + if (mype == 0) print *,k,an_grid%ak(k),an_grid%bk(k),taper,minval(an_grid%var3d(:,:,k)),maxval(an_grid%var3d(:,:,k)) + enddo + endif + call fv3_netcdf_write_var3d(ncdat,output_vars(ivar),an_grid%var3d) enddo var_loop diff --git a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/gfs_ncio_interface.f90 b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/gfs_ncio_interface.f90 index 4d3d07236b..ddabc66a40 100644 --- a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/gfs_ncio_interface.f90 +++ b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/gfs_ncio_interface.f90 @@ -137,6 +137,8 @@ subroutine gfs_ncio_initialize(meta_ncio,filename) allocate(meta_ncio%lon(meta_ncio%dimx*meta_ncio%dimy)) if (.not. allocated(meta_ncio%lat)) & allocate(meta_ncio%lat(meta_ncio%dimx*meta_ncio%dimy)) + if (.not. allocated(meta_ncio%vcoord)) & + allocate(meta_ncio%vcoord(meta_ncio%dimz+1,2)) call read_vardata(gfile,'lon', tmp2d) meta_ncio%lon = reshape(tmp2d, (/meta_ncio%dimx*meta_ncio%dimy/)) call read_vardata(gfile,'lat', tmp2d) @@ -148,6 +150,10 @@ subroutine gfs_ncio_initialize(meta_ncio,filename) meta_ncio%idvm=1 meta_ncio%ntrac = 8 meta_ncio%ncldt = 5 + call read_attribute(gfile,'ak',tmp1d) + meta_ncio%vcoord(:,1) = tmp1d(:) + call read_attribute(gfile,'bk',tmp1d) + meta_ncio%vcoord(:,2) = tmp1d(:) call read_vardata(gfile,'time',tmp1d) meta_ncio%fhour = nint(tmp1d(1)) diff --git a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/namelist_def.f90 b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/namelist_def.f90 index 4f20d219e9..c0ebc34e7e 100644 --- a/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/namelist_def.f90 +++ b/util/EnKF/gfs/src/calc_increment_ens_ncio.fd/namelist_def.f90 @@ -8,7 +8,7 @@ module namelist_def public :: analysis_filename, firstguess_filename, increment_filename public :: datapath public :: debug - public :: do_icmr + public :: do_icmr, taper_strat, ak_bot, ak_top public :: incvars_to_zero public :: read_namelist public :: write_namelist @@ -26,9 +26,13 @@ module namelist_def character(len=12) :: incvars_to_zero(max_vars) = 'NONE' logical :: do_icmr = .false. + logical :: taper_strat = .false. + ! damp humidity increments between these two levels if taper_strat=T + real :: ak_bot = 10000. ! units Pa + real :: ak_top = 5000. namelist /setup/ datapath, analysis_filename, firstguess_filename, increment_filename, & - nens, debug, imp_physics + nens, debug, imp_physics, ak_top, ak_bot, taper_strat namelist /zeroinc/ incvars_to_zero contains diff --git a/util/EnKF/gfs/src/calc_increment_ncio.fd/calc_increment_ncio.f90 b/util/EnKF/gfs/src/calc_increment_ncio.fd/calc_increment_ncio.f90 index dd14663ef5..4928a62e6c 100755 --- a/util/EnKF/gfs/src/calc_increment_ncio.fd/calc_increment_ncio.f90 +++ b/util/EnKF/gfs/src/calc_increment_ncio.fd/calc_increment_ncio.f90 @@ -21,6 +21,14 @@ PROGRAM calc_increment_ncio ! 4th command line arg is logical for controlling whether microphysics ! increment is computed. +! 5th command line arg is logical for controlling whether delz +! increment should be computed +! 6th command line arg is logical for controlling whether humidity +! and microphysics vars should be tapered to zero in stratosphere. +! The vertical profile of the taper is controlled by ak_top and ak_bot. + +! If delp and/or delz are not in the background history files, then +! their increments are inferred (from ps, T and humidity increments). ! ! attributes: ! language: f95 @@ -48,21 +56,24 @@ PROGRAM calc_increment_ncio real, allocatable, dimension(:,:) :: values_2d_fg,values_2d_anal,values_2d_inc,& ps_fg, ps_anal real, allocatable, dimension(:,:,:) :: values_3d_fg,values_3d_anal,values_3d_inc,& - q_fg, q_anal, tmp_fg, tmp_anal, delzb, delza + taper_vert,q_fg, q_anal, tmp_fg, tmp_anal, delzb, delza type(Dataset) :: dset_anal,dset_fg type(Dimension) :: londim,latdim,levdim integer, dimension(3) :: dimid_3d integer, dimension(1) :: dimid_1d integer varid_lon,varid_lat,varid_lev,varid_ilev,varid_hyai,varid_hybi,& dimid_lon,dimid_lat,dimid_lev,dimid_ilev,ncfileid,ncstatus - logical :: no_mpinc, no_delzinc, has_dpres, has_delz + logical :: no_mpinc, no_delzinc, has_dpres, has_delz, taper_strat character(len=10) :: bufchar - real rd,rv,fv,grav + real rd,rv,fv,grav,ak_bot,ak_top rd = 2.8705e+2 rv = 4.6150e+2 fv = rv/rd-1. ! used in virtual temperature equation grav = 9.80665 + ! damp humidity increments between these two levels if taper_strat=T + ak_bot = 10000. ! units Pa + ak_top = 5000. call getarg(1,filename_fg) ! first guess ncio file call getarg(2,filename_anal) ! analysis ncio file @@ -71,6 +82,8 @@ PROGRAM calc_increment_ncio read(bufchar,'(L)') no_mpinc ! if T, no microphysics increments computed call getarg(5, bufchar) read(bufchar,'(L)') no_delzinc ! if T, no delz increments computed + call getarg(6, bufchar) + read(bufchar,'(L)') taper_strat ! if T, taper sphum,liq_wat,ice_wat in strat write(6,*)'CALC_INCREMENT_NCIO:' write(6,*)'filename_fg=',trim(filename_fg) @@ -78,6 +91,7 @@ PROGRAM calc_increment_ncio write(6,*)'filename_inc=',trim(filename_inc) write(6,*)'no_mpinc',no_mpinc write(6,*)'no_delzinc',no_delzinc + write(6,*)'taper_strat',taper_strat dset_fg = open_dataset(trim(filename_fg),errcode=iret) if (iret .ne. 0) then @@ -259,6 +273,7 @@ PROGRAM calc_increment_ncio ! ps increment. allocate(values_2d_inc(nlons,nlats)) + allocate(taper_vert(nlons,nlats,nlevs)) allocate(values_3d_inc(nlons,nlats,nlevs)) do nvar=1,dset_fg%nvars ndims = dset_fg%variables(nvar)%ndims @@ -269,6 +284,19 @@ PROGRAM calc_increment_ncio values_2d_inc(:,nlats:1:-1) = values_2d_anal - values_2d_fg endif enddo + ! taper function for humidity, ice and liq water increments. + taper_vert=1. + if (taper_strat) print *,'profile to taper strat humid inc (k,ak,bk,taper):' + do k=1,nlevs + if (k < nlevs/2 .and. (ak(k) <= ak_bot .and. ak(k) >= ak_top)) then + taper_vert(:,:,k)= (ak(k) - ak_top)/(ak_bot - ak_top) + else if (bk(k) .eq. 0. .and. ak(k) < ak_top) then + taper_vert(:,:,k) = 0. + endif + if (taper_strat) then + print *,k,ak(k),bk(k),taper_vert(1,1,k) + endif + enddo do nvar=1,dset_fg%nvars ndims = dset_fg%variables(nvar)%ndims @@ -302,7 +330,13 @@ PROGRAM calc_increment_ncio call read_vardata(dset_fg,trim(dset_fg%variables(nvar)%name),values_3d_fg) call read_vardata(dset_anal,trim(dset_fg%variables(nvar)%name),values_3d_anal) ! increment (flip lats) - values_3d_inc(:,nlats:1:-1,:) = values_3d_anal - values_3d_fg + if (taper_strat .and. (trim(ncvarname) .eq. 'sphum_inc' .or. & + trim(ncvarname) .eq. 'liq_wat_inc' .or. & + trim(ncvarname) .eq. 'ice_wat_inc')) then + values_3d_inc(:,nlats:1:-1,:) = taper_vert*(values_3d_anal - values_3d_fg) + else + values_3d_inc(:,nlats:1:-1,:) = values_3d_anal - values_3d_fg + endif call write_ncdata3d(values_3d_inc,ncvarname,nlons,nlats,nlevs,ncfileid,dimid_3d) endif endif ! ndims == 4 diff --git a/util/EnKF/gfs/src/getsigensmeanp_smooth.fd/CMakeLists.txt b/util/EnKF/gfs/src/getsigensmeanp_smooth.fd/CMakeLists.txt index 08ab0defdf..ce102ee833 100644 --- a/util/EnKF/gfs/src/getsigensmeanp_smooth.fd/CMakeLists.txt +++ b/util/EnKF/gfs/src/getsigensmeanp_smooth.fd/CMakeLists.txt @@ -9,4 +9,5 @@ if(BUILD_UTIL) SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_Fortran_FLAGS}" ) include_directories( ${NEMSIOINC} ${SIGIOINC} ${NETCDF_INCLUDE_DIRS} ${FV3GFS_NCIO_INCS}) target_link_libraries( getsigensmeanp_smooth.x ${FV3GFS_NCIO_LIBRARIES} ${BACIO_LIBRARY} ${NEMSIO_LIBRARY} ${BACIO_LIBRARY} ${SIGIO_LIBRARY} ${W3NCO_4_LIBRARY} ${SP_4_LIBRARY} ${MPI_Fortran_LIBRARIES} ${NETCDF_Fortran_LIBRARIES} ${NETCDF_C_LIBRARIES}) + endif() diff --git a/util/EnKF/gfs/src/recenterens_ncio.fd/CMakeLists.txt b/util/EnKF/gfs/src/recenterens_ncio.fd/CMakeLists.txt new file mode 100644 index 0000000000..7e347146f2 --- /dev/null +++ b/util/EnKF/gfs/src/recenterens_ncio.fd/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 2.6) +if(BUILD_UTIL) + file(GLOB LOCAL_SRC ${CMAKE_CURRENT_SOURCE_DIR}/*.f90) + set_source_files_properties( ${LOCAL_SRC} PROPERTIES COMPILE_FLAGS ${UTIL_Fortran_FLAGS} ) + add_executable(recenterens_ncio.x ${LOCAL_SRC} ) + set_target_properties( recenterens_ncio.x PROPERTIES COMPILE_FLAGS ${UTIL_Fortran_FLAGS} ) + SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_Fortran_FLAGS}" ) + include_directories( ${NETCDF_INCLUDES} ${FV3GFS_NCIO_INCS}) + target_link_libraries( recenterens_ncio.x ${FV3GFS_NCIO_LIBRARIES} ${W3NCO_4_LIBRARY} ${MPI_Fortran_LIBRARIES} ${NETCDF_Fortran_LIBRARIES} ${NETCDF_LIBRARIES} ${HDF5_HL_LIBRARIES} ${HDF5_LIBRARIES}) +endif() diff --git a/util/EnKF/gfs/src/recenterens_ncio.fd/recenterens_ncio.f90 b/util/EnKF/gfs/src/recenterens_ncio.fd/recenterens_ncio.f90 new file mode 100644 index 0000000000..6afbb572eb --- /dev/null +++ b/util/EnKF/gfs/src/recenterens_ncio.fd/recenterens_ncio.f90 @@ -0,0 +1,183 @@ +program recenterens_ncio +!$$$ main program documentation block +! +! program: recenterens_ncio recenter +! +! prgmmr: whitaker org: esrl/psd date: 2009-02-23 +! +! abstract: Read ensemble from netcdf history files, +! remove mean specified from another file, add a +! new mean specified from a third file, and write +! out result to a fourth file. 'Partial' recentering +! by weighting old and new means. +! +! program history log: +! 2020-11-04 Initial version. +! +! usage: +! input files: +! +! output files: +! +! attributes: +! language: f95 +! +! +!$$$ + + use module_fv3gfs_ncio, only: open_dataset, create_dataset, read_attribute, & + Dataset, Dimension, close_dataset, has_attr, has_var, & + read_vardata, write_attribute, write_vardata, & + get_dim, quantize_data + + implicit none + + include "mpif.h" + + real,parameter:: zero=0.0_4 + + logical:: quantize + character(len=500) filename_meani,filename_meano,filenamein,filenameout + character(len=3) charnanal, charwgt_ensmean, charwgt_control + character(len=4) charnin + integer iret,mype,mype1,npe,nanals,ierr, iwgt_ensmean, iwgt_control + integer:: latb,lonb,levs,nbits,nvar,ndims + real(4),allocatable, dimension(:,:) :: values_2d, values_2d_i, values_2d_mi,& + values_2d_mo + real(4),allocatable, dimension(:,:,:) :: values_3d, values_3d_i, values_3d_mi,& + values_3d_mo + real(4) compress_err, rwgt_control, rwgt_ensmean + + type(Dataset) :: dseti,dseto,dsetmi,dsetmo + type(Dimension) :: londim,latdim,levdim + +! Initialize mpi + call MPI_Init(ierr) + +! mype is process number, npe is total number of processes. + call MPI_Comm_rank(MPI_COMM_WORLD,mype,ierr) + call MPI_Comm_size(MPI_COMM_WORLD,npe,ierr) + + if (mype==0) call w3tagb('RECENTERENS_NCIO',2011,0319,0055,'NP25') + +! read data from this file + call getarg(1,filenamein) ! increment or analysis + +! subtract this mean + call getarg(2,filename_meani) ! mean increment or analysis + +! then add to this mean + call getarg(3,filename_meano) ! new mean analysis + +! and put in this file. + call getarg(4,filenameout) ! new increment of analysis + +! how many ensemble members to process + call getarg(5,charnin) + read(charnin,'(i4)') nanals + +! weight given to original ens mean + call getarg(6,charwgt_ensmean) + read(charwgt_ensmean,'(i3)') iwgt_ensmean + rwgt_ensmean = iwgt_ensmean/100. + +! weight given to new ens mean + call getarg(7,charwgt_control) + read(charwgt_control,'(i3)') iwgt_control + rwgt_control = iwgt_control/100. + + if (mype==0) then + write(6,*)'RECENTERENS_NCIO: PROCESS ',nanals,' ENSEMBLE MEMBERS' + write(6,*)'filenamein=',trim(filenamein) + write(6,*)'filename_meani=',trim(filename_meani) + write(6,*)'filename_meano=',trim(filename_meano) + write(6,*)'filenameout=',trim(filenameout) + write(6,*)'rwgt_ensmean,rwgt_control=',rwgt_ensmean,rwgt_control + endif + + mype1 = mype+1 + if (mype1 <= nanals) then + + dsetmi = open_dataset(filename_meani,errcode=iret) + + londim = get_dim(dsetmi,'grid_xt'); lonb = londim%len + latdim = get_dim(dsetmi,'grid_yt'); latb = latdim%len + levdim = get_dim(dsetmi,'pfull'); levs = levdim%len + write(charnanal,'(i3.3)') mype1 + dsetmo = open_dataset(filename_meano) + dseti = open_dataset(trim(filenamein)//"_mem"//charnanal) + dseto = create_dataset(trim(filenameout)//"_mem"//charnanal, dseti, copy_vardata=.true.) + do nvar=1,dseti%nvars + ndims = dseti%variables(nvar)%ndims + if (ndims > 2) then + if (ndims == 3 .and. trim(dseti%variables(nvar)%name) /= 'hgtsfc') then + ! pressfc + if (mype == 0) print *,'recentering ',& + trim(dseti%variables(nvar)%name) + call read_vardata(dseti,trim(dseti%variables(nvar)%name),values_2d_i) + call read_vardata(dsetmi,trim(dseti%variables(nvar)%name),values_2d_mi) + call read_vardata(dsetmo,trim(dseti%variables(nvar)%name),values_2d_mo) + values_2d = values_2d_i - values_2d_mi + rwgt_ensmean*values_2d_mi + rwgt_control*values_2d_mo + if (has_attr(dseti, 'nbits', trim(dseti%variables(nvar)%name))) then + call read_attribute(dseti, 'nbits', nbits, & + trim(dseti%variables(nvar)%name)) + quantize = .true. + if (nbits < 1) quantize = .false. + else + quantize = .false. + endif + if (quantize) then + values_2d_mi = values_2d + call quantize_data(values_2d_mi, values_2d, nbits, compress_err) + call write_attribute(dseto,& + 'max_abs_compression_error',compress_err,trim(dseti%variables(nvar)%name)) + endif + call write_vardata(dseto,trim(dseti%variables(nvar)%name),values_2d) + else if (ndims == 4) then + if (mype == 0) print *,'recentering ',& + trim(dseti%variables(nvar)%name) + call read_vardata(dseti,trim(dseti%variables(nvar)%name),values_3d_i) + call read_vardata(dsetmi,trim(dseti%variables(nvar)%name),values_3d_mi) + call read_vardata(dsetmo,trim(dseti%variables(nvar)%name),values_3d_mo) + values_3d = values_3d_i - values_3d_mi + rwgt_ensmean*values_3d_mi + rwgt_control*values_3d_mo + if (has_attr(dseti, 'nbits', trim(dseti%variables(nvar)%name))) then + call read_attribute(dseti, 'nbits', nbits, & + trim(dseti%variables(nvar)%name)) + quantize = .true. + if (nbits < 1) quantize = .false. + else + quantize = .false. + endif + if (quantize) then + values_3d_mi = values_3d + call quantize_data(values_3d_mi, values_3d, nbits, compress_err) + call write_attribute(dseto,& + 'max_abs_compression_error',compress_err,trim(dseti%variables(nvar)%name)) + endif + call write_vardata(dseto,trim(dseti%variables(nvar)%name),values_3d) + endif + endif ! ndims > 2 + enddo ! nvars + + deallocate(values_2d,values_2d_i,values_2d_mi,values_2d_mo) + deallocate(values_3d,values_3d_i,values_3d_mi,values_3d_mo) + call close_dataset(dsetmi) + call close_dataset(dsetmo) + call close_dataset(dseti) + call close_dataset(dseto) + +! Jump here if more mpi processors than files to process + else + write (6,*) 'no files to process for mpi task = ',mype + end if ! end if mype + + call MPI_Barrier(MPI_COMM_WORLD,ierr) + + if (mype==0) call w3tage('RECENTERENS_NCIO') + + call MPI_Finalize(ierr) + if (mype .eq. 0 .and. ierr .ne. 0) then + print *, 'MPI_Finalize error status = ',ierr + end if + +END program recenterens_ncio diff --git a/util/EnKF/gfs/src/recenterncio_hybgain.fd/recenterncio_hybgain.f90 b/util/EnKF/gfs/src/recenterncio_hybgain.fd/recenterncio_hybgain.f90 index aa063b66d4..5ab94f071c 100644 --- a/util/EnKF/gfs/src/recenterncio_hybgain.fd/recenterncio_hybgain.f90 +++ b/util/EnKF/gfs/src/recenterncio_hybgain.fd/recenterncio_hybgain.f90 @@ -6,7 +6,8 @@ program recenterncio_hybgain ! prgmmr: whitaker org: esrl/psd date: 2009-02-23 ! ! abstract: Recenter ensemble analysis files about new -! mean, computed from blended 3DVar and EnKF increments. +! mean, computed from blended 3DVar and EnKF increments +! (optionally applying RTPS inflation). ! ! program history log: ! 2019-02-10 Initial version. @@ -35,20 +36,23 @@ program recenterncio_hybgain external :: MPI_Init, MPI_Comm_rank, MPI_Comm_size, w3tagb, MPI_Abort,& MPI_Barrier, w3tage, MPI_Finalize - character*500 filename_fg,filename_varanal,filename_enkfanal,filenamein,& - filenameout,filename_anal,filename + character*500 filename_varfg,filename_varanal,filename_enkfanal,filenamein,& + filename_enkffg,filenameout,filename_anal,filename,filename_fsprd,filename_asprd character*3 charnanal character(len=4) charnin - integer mype,mype1,npe,nanals,iret,ialpha,ibeta + integer mype,mype1,npe,nanals,iret,ialpha,ibeta,irtps,i,j,k integer:: nlats,nlons,nlevs,nvar,ndims,nbits - real alpha,beta - real(4),allocatable,dimension(:,:) :: values_2d_varanal,values_2d_enkfanal,values_2d_fg,values_2d_anal,& - values_2d_tmp, values_2d - real(4),allocatable,dimension(:,:,:) :: values_3d_varanal,values_3d_enkfanal,values_3d_fg,values_3d_anal,& - values_3d_tmp, values_3d + real alpha,beta,rtps,infmin,infmax,clip + real(4),allocatable,dimension(:,:) :: values_2d_varanal,values_2d_enkfanal,& + values_2d_varfg,values_2d_enkffg,values_2d_anal,& + values_2d,asprd_2d,fsprd_2d,inf_2d,values_2d_tmp + real(4),allocatable,dimension(:,:,:) :: values_3d_varanal,values_3d_enkfanal,& + values_3d_varfg,values_3d_enkffg,values_3d_anal,& + values_3d,asprd_3d,fsprd_3d,inf_3d,values_3d_tmp real(4) compress_err - type(Dataset) :: dseti,dseto,dset_anal,dset_fg,dset_varanal,dset_enkfanal + type(Dataset) :: dseti,dseto,dset_blendanal,dset_varfg,dset_varanal,dset_enkffg,dset_enkfanal,dset_asprd,dset_fsprd type(Dimension) :: londim,latdim,levdim + logical cliptracers,tracer ! Initialize mpi call MPI_Init(iret) @@ -59,66 +63,103 @@ program recenterncio_hybgain if (mype==0) call w3tagb('RECENTERNCIO_HYBGAIN',2011,0319,0055,'NP25') - call getarg(1,filename_fg) ! first guess ensmean background netcdf file + call getarg(1,filename_varfg) ! first guess 3dvar netcdf file call getarg(2,filename_varanal) ! 3dvar analysis - call getarg(3,filename_enkfanal) ! enkf mean analysis - call getarg(4,filename_anal) ! blended analysis (to recenter ensemble around) - call getarg(5,filenamein) ! prefix for input ens member files (append _mem###) - call getarg(6,filenameout) ! prefix for output ens member files (append _mem###) + call getarg(3,filename_enkffg) ! first guess enkf netcdf file + call getarg(4,filename_enkfanal) ! enkf mean analysis + call getarg(5,filename_anal) ! blended analysis (to recenter ensemble around) + call getarg(6,filenamein) ! prefix for input ens member files (append _mem###) + call getarg(7,filenameout) ! prefix for output ens member files (append _mem###) ! blending coefficients - call getarg(7,charnin) + call getarg(8,charnin) read(charnin,'(i4)') ialpha ! wt for varanal (3dvar) alpha = ialpha/1000. - call getarg(8,charnin) + call getarg(9,charnin) read(charnin,'(i4)') ibeta ! wt for enkfanal (enkf) beta = ibeta/1000. -! new_anal = fg + alpha*(varanal-fg) + beta(enkfanal-fg) -! = (1.-alpha-beta)*fg + alpha*varanal + beta*enkfanal + call getarg(10,charnin) + read(charnin,'(i4)') irtps ! rtps relaxation coeff + rtps = irtps/1000. +! new_anal = fg + alpha*(varanal-fg) + beta*(enkfanal-fg) ! how many ensemble members to process - call getarg(9,charnin) + call getarg(11,charnin) read(charnin,'(i4)') nanals + if (rtps > 0) then + call getarg(13,filename_fsprd) ! first guess ensemble spread + call getarg(14,filename_asprd) ! analysis ensemble spread + endif + + infmin=1.0; infmax=10. + clip = tiny(rtps) if (mype==0) then write(6,*)'RECENTERNCIO_HYBGAIN: PROCESS ',nanals,' ENSEMBLE MEMBERS' - write(6,*)'ens mean background in ',trim(filename_fg) + write(6,*)'3dvar background in ',trim(filename_varfg) + write(6,*)'EnKF background in ',trim(filename_enkffg) write(6,*)'3dvar analysis in ',trim(filename_varanal) write(6,*)'EnKF mean analysis in ',trim(filename_enkfanal) write(6,*)'Blended mean analysis to be written to ',trim(filename_anal) write(6,*)'Prefix for member input files ',trim(filenamein) write(6,*)'Prefix for member output files ',trim(filenameout) - write(6,*)'3dvar weight, EnKF weight =',alpha,beta + write(6,*)'3dvar weight, EnKF weight, RTPS relaxation =',alpha,beta,rtps + if (rtps > 0) then + write(6,*)'ens spread background in ',trim(filename_fsprd) + write(6,*)'ens spread posterior in ',trim(filename_asprd) + endif endif mype1 = mype+1 if (mype1 <= nanals) then - dset_fg = open_dataset(filename_fg,errcode=iret) + dset_varfg = open_dataset(filename_varfg,errcode=iret) if (iret == 0 ) then - if (mype == 0) write(6,*)'Read netcdf ',trim(filename_fg) - londim = get_dim(dset_fg,'grid_xt'); nlons = londim%len - latdim = get_dim(dset_fg,'grid_yt'); nlats = latdim%len - levdim = get_dim(dset_fg,'pfull'); nlevs = levdim%len + if (mype == 0) write(6,*)'Read netcdf ',trim(filename_varfg) + londim = get_dim(dset_varfg,'grid_xt'); nlons = londim%len + latdim = get_dim(dset_varfg,'grid_yt'); nlats = latdim%len + levdim = get_dim(dset_varfg,'pfull'); nlevs = levdim%len if (mype == 0) write(6,*)' nlons=',nlons,' nlats=',nlats,' nlevs=',nlevs else - write(6,*) 'error opening ',trim(filename_fg) + write(6,*) 'error opening ',trim(filename_varfg) call MPI_Abort(MPI_COMM_WORLD,98,iret) stop endif - ! readin in 3dvar, enkf analyses, plus ens mean background, blend + ! read in 3dvar, enkf analyses and backgrounds, blend increments + if (alpha > 0) then dset_varanal = open_dataset(filename_varanal,errcode=iret) if (iret /= 0) then print *,'error opening ',trim(filename_varanal) call MPI_Abort(MPI_COMM_WORLD,98,iret) stop endif + endif dset_enkfanal = open_dataset(filename_enkfanal,errcode=iret) if (iret /= 0) then print *,'error opening ',trim(filename_enkfanal) call MPI_Abort(MPI_COMM_WORLD,98,iret) stop endif - if (mype == 0) then - dset_anal = create_dataset(filename_anal, dset_enkfanal, & + dset_enkffg = open_dataset(filename_enkffg,errcode=iret) + if (iret /= 0) then + print *,'error opening ',trim(filename_enkffg) + call MPI_Abort(MPI_COMM_WORLD,98,iret) + stop + endif + if (rtps > 0) then + dset_fsprd = open_dataset(filename_fsprd,errcode=iret) + if (iret /= 0) then + print *,'error opening ',trim(filename_fsprd) + call MPI_Abort(MPI_COMM_WORLD,98,iret) + stop + endif + dset_asprd = open_dataset(filename_asprd,errcode=iret) + if (iret /= 0) then + print *,'error opening ',trim(filename_asprd) + call MPI_Abort(MPI_COMM_WORLD,98,iret) + stop + endif + endif + if (mype == 0 .and. alpha > 0) then + dset_blendanal = create_dataset(filename_anal, dset_enkfanal, & copy_vardata=.true., errcode=iret) if (iret /= 0) then print *,'error opening ',trim(filename_anal) @@ -142,88 +183,151 @@ program recenterncio_hybgain call MPI_Abort(MPI_COMM_WORLD,98,iret) stop endif + allocate(inf_2d(nlons,nlats),fsprd_2d(nlons,nlats),asprd_2d(nlons,nlats)) + allocate(inf_3d(nlons,nlats,nlevs),fsprd_3d(nlons,nlats,nlevs),asprd_3d(nlons,nlats,nlevs)) - do nvar=1,dset_fg%nvars - ndims = dset_fg%variables(nvar)%ndims + do nvar=1,dset_varfg%nvars + ndims = dset_varfg%variables(nvar)%ndims + if (trim(dset_varfg%variables(nvar)%name) == 'spfh' .or. & + trim(dset_varfg%variables(nvar)%name) == 'o3mr' .or. & + trim(dset_varfg%variables(nvar)%name) == 'clwmr' .or. & + trim(dset_varfg%variables(nvar)%name) == 'icmr') then + tracer = .true. + else + tracer = .false. + endif if (ndims > 2) then - if (ndims == 3 .and. trim(dset_fg%variables(nvar)%name) /= 'hgtsfc') then + if (ndims == 3 .and. trim(dset_varfg%variables(nvar)%name) /= 'hgtsfc') then ! pressfc - call read_vardata(dset_fg,trim(dset_fg%variables(nvar)%name),values_2d_fg) - call read_vardata(dset_varanal,trim(dset_fg%variables(nvar)%name),values_2d_varanal) - call read_vardata(dset_enkfanal,trim(dset_fg%variables(nvar)%name),values_2d_enkfanal) - call read_vardata(dseti,trim(dset_fg%variables(nvar)%name),values_2d) + call read_vardata(dset_varfg,trim(dset_varfg%variables(nvar)%name),values_2d_varfg) + call read_vardata(dset_enkffg,trim(dset_enkffg%variables(nvar)%name),values_2d_enkffg) + if (alpha>0) call read_vardata(dset_varanal,trim(dset_varfg%variables(nvar)%name),values_2d_varanal) + call read_vardata(dset_enkfanal,trim(dset_enkffg%variables(nvar)%name),values_2d_enkfanal) + call read_vardata(dseti,trim(dset_enkffg%variables(nvar)%name),values_2d) ! blended analysis - values_2d_anal = (1.-alpha-beta)*values_2d_fg + & - alpha*values_2d_varanal + & - beta*values_2d_enkfanal + values_2d_anal = values_2d_enkffg + beta*(values_2d_enkfanal-values_2d_enkffg) + if (alpha > 0) & + values_2d_anal = values_2d_anal + alpha*(values_2d_varanal-values_2d_varfg) ! recentered ensemble member - values_2d = values_2d - values_2d_enkfanal + values_2d_anal - if (has_attr(dset_fg, 'nbits', trim(dset_fg%variables(nvar)%name))) then - call read_attribute(dset_fg, 'nbits', nbits, & - trim(dset_fg%variables(nvar)%name),errcode=iret) + if (rtps > 0) then ! RTPS inflation + call read_vardata(dset_fsprd,trim(dset_enkffg%variables(nvar)%name),fsprd_2d) + call read_vardata(dset_asprd,trim(dset_enkffg%variables(nvar)%name),asprd_2d) + fsprd_2d = max(fsprd_2d,tiny(fsprd_2d)) + asprd_2d = max(asprd_2d,tiny(asprd_2d)) + inf_2d = rtps*((fsprd_2d-asprd_2d)/asprd_2d) + 1.0 + do j=1,nlats + do i=1,nlons + inf_2d(i,j) = max(infmin,min(inf_2d(i,j),infmax)) + enddo + enddo + values_2d = inf_2d*(values_2d - values_2d_enkfanal) + values_2d_anal + if (mype == 0) & + print *,'min/max ',trim(dset_enkffg%variables(nvar)%name),& + ' inflation = ',minval(inf_2d),maxval(inf_2d) + else + values_2d = values_2d - values_2d_enkfanal + values_2d_anal + endif + if (has_attr(dset_enkffg, 'nbits', trim(dset_enkffg%variables(nvar)%name))) then + call read_attribute(dset_enkffg, 'nbits', nbits, & + trim(dset_enkffg%variables(nvar)%name),errcode=iret) else iret = 1 endif - if (mype == 0) then ! write out blended analysis on root task + if (mype == 0 .and. alpha > 0) then ! write out blended analysis on root task if (iret == 0 .and. nbits > 0) then values_2d_tmp = values_2d_anal call quantize_data(values_2d_tmp, values_2d_anal, nbits, compress_err) - call write_attribute(dset_anal,& - 'max_abs_compression_error',compress_err,trim(dset_fg%variables(nvar)%name)) + call write_attribute(dset_blendanal,& + 'max_abs_compression_error',compress_err,trim(dset_enkffg%variables(nvar)%name)) endif - call write_vardata(dset_anal,trim(dset_fg%variables(nvar)%name),values_2d_anal) + call write_vardata(dset_blendanal,trim(dset_enkffg%variables(nvar)%name),values_2d_anal) endif if (iret == 0 .and. nbits > 0) then values_2d_tmp = values_2d call quantize_data(values_2d_tmp, values_2d, nbits, compress_err) call write_attribute(dseto,& - 'max_abs_compression_error',compress_err,trim(dset_fg%variables(nvar)%name)) + 'max_abs_compression_error',compress_err,trim(dset_enkffg%variables(nvar)%name)) endif - call write_vardata(dseto,trim(dset_fg%variables(nvar)%name),values_2d) + if (tracer) then + if (cliptracers) where (values_2d < clip) values_2d = clip + if (mype == 0) & + print *,'clipping ',trim(dset_enkffg%variables(nvar)%name) + endif + call write_vardata(dseto,trim(dset_enkffg%variables(nvar)%name),values_2d) else if (ndims == 4) then - call read_vardata(dset_fg,trim(dset_fg%variables(nvar)%name),values_3d_fg) - call read_vardata(dset_varanal,trim(dset_fg%variables(nvar)%name),values_3d_varanal) - call read_vardata(dset_enkfanal,trim(dset_fg%variables(nvar)%name),values_3d_enkfanal) - call read_vardata(dseti,trim(dset_fg%variables(nvar)%name),values_3d) + call read_vardata(dset_varfg,trim(dset_varfg%variables(nvar)%name),values_3d_varfg) + call read_vardata(dset_enkffg,trim(dset_enkffg%variables(nvar)%name),values_3d_enkffg) + if (alpha>0) call read_vardata(dset_varanal,trim(dset_varfg%variables(nvar)%name),values_3d_varanal) + call read_vardata(dset_enkfanal,trim(dset_enkffg%variables(nvar)%name),values_3d_enkfanal) + call read_vardata(dseti,trim(dset_enkffg%variables(nvar)%name),values_3d) ! blended analysis - values_3d_anal = (1.-alpha-beta)*values_3d_fg + & - alpha*values_3d_varanal + & - beta*values_3d_enkfanal + values_3d_anal = values_3d_enkffg + beta*(values_3d_enkfanal-values_3d_enkffg) + if (alpha > 0) & + values_3d_anal = values_3d_anal + alpha*(values_3d_varanal-values_3d_varfg) ! recentered ensemble member - values_3d = values_3d - values_3d_enkfanal + values_3d_anal - if (has_attr(dset_fg, 'nbits', trim(dset_fg%variables(nvar)%name))) then - call read_attribute(dset_fg, 'nbits', nbits, & - trim(dset_fg%variables(nvar)%name),errcode=iret) + if (rtps > 0) then ! RTPS inflation + call read_vardata(dset_fsprd,trim(dset_enkffg%variables(nvar)%name),fsprd_3d) + call read_vardata(dset_asprd,trim(dset_enkffg%variables(nvar)%name),asprd_3d) + fsprd_3d = max(fsprd_3d,tiny(fsprd_3d)) + asprd_3d = max(asprd_3d,tiny(asprd_3d)) + inf_3d = rtps*((fsprd_3d-asprd_3d)/asprd_3d) + 1.0 + do k=1,nlevs + do j=1,nlats + do i=1,nlons + inf_3d(i,j,k) = max(infmin,min(inf_3d(i,j,k),infmax)) + enddo + enddo + enddo + values_3d = inf_3d*(values_3d - values_3d_enkfanal) + values_3d_anal + if (mype == 0) & + print *,'min/max ',trim(dset_enkffg%variables(nvar)%name),& + ' inflation = ',minval(inf_3d),maxval(inf_3d) + else + values_3d = values_3d - values_3d_enkfanal + values_3d_anal + endif + if (has_attr(dset_enkffg, 'nbits', trim(dset_enkffg%variables(nvar)%name))) then + call read_attribute(dset_enkffg, 'nbits', nbits, & + trim(dset_enkffg%variables(nvar)%name),errcode=iret) else iret = 1 endif - if (mype == 0) then ! write out blended analysis on root task + if (mype == 0 .and. alpha > 0) then ! write out blended analysis on root task if (iret == 0 .and. nbits > 0) then values_3d_tmp = values_3d_anal call quantize_data(values_3d_tmp, values_3d_anal, nbits, compress_err) - call write_attribute(dset_anal,& - 'max_abs_compression_error',compress_err,trim(dset_fg%variables(nvar)%name)) + call write_attribute(dset_blendanal,& + 'max_abs_compression_error',compress_err,trim(dset_enkffg%variables(nvar)%name)) endif - call write_vardata(dset_anal,trim(dset_fg%variables(nvar)%name),values_3d_anal) + call write_vardata(dset_blendanal,trim(dset_enkffg%variables(nvar)%name),values_3d_anal) endif if (iret == 0 .and. nbits > 0) then values_3d_tmp = values_3d call quantize_data(values_3d_tmp, values_3d, nbits, compress_err) call write_attribute(dseto,& - 'max_abs_compression_error',compress_err,trim(dset_fg%variables(nvar)%name)) + 'max_abs_compression_error',compress_err,trim(dset_enkffg%variables(nvar)%name)) + endif + if (tracer) then + if (cliptracers) where (values_3d < clip) values_3d = clip + if (mype == 0) & + print *,'clipping ',trim(dset_enkffg%variables(nvar)%name) endif - call write_vardata(dseto,trim(dset_fg%variables(nvar)%name),values_3d) + call write_vardata(dseto,trim(dset_enkffg%variables(nvar)%name),values_3d) endif endif ! ndims > 2 enddo ! nvars - if (mype == 0) call close_dataset(dset_anal) + if (mype == 0 .and. alpha > 0) call close_dataset(dset_blendanal) call close_dataset(dseti) call close_dataset(dseto) - call close_dataset(dset_fg) - call close_dataset(dset_varanal) + call close_dataset(dset_enkffg) + call close_dataset(dset_varfg) + if (alpha > 0) call close_dataset(dset_varanal) call close_dataset(dset_enkfanal) + if (rtps > 0) then + call close_dataset(dset_fsprd) + call close_dataset(dset_asprd) + endif write(6,*)'task mype=',mype,' process ',trim(filenameout)//"_mem"//charnanal,' iret=',iret ! Jump here if more mpi processors than files to process @@ -232,9 +336,12 @@ program recenterncio_hybgain end if ! end if mype 100 continue + deallocate(asprd_2d,asprd_3d) + deallocate(fsprd_2d,fsprd_3d) + deallocate(inf_2d,inf_3d) call MPI_Barrier(MPI_COMM_WORLD,iret) - if (mype==0) call w3tage('RECENTERSIGP_HYBGAIN') + if (mype==0) call w3tage('RECENTERSIGP_HYBGAIN2') call MPI_Finalize(iret) if (mype == 0 .and. iret /= 0) then diff --git a/util/EnKF/gfs/src/recentersigp.fd/recentersigp.f90 b/util/EnKF/gfs/src/recentersigp.fd/recentersigp.f90 index bc01dbd984..321845bf35 100644 --- a/util/EnKF/gfs/src/recentersigp.fd/recentersigp.f90 +++ b/util/EnKF/gfs/src/recentersigp.fd/recentersigp.f90 @@ -82,23 +82,23 @@ program recentersigp NSIGO=61 ! read data from this file - call getarg(1,filenamein) + call getarg(1,filenamein) ! increment or analysis ! subtract this mean - call getarg(2,filename_meani) + call getarg(2,filename_meani) ! mean increment or analysis ! then add to this mean - call getarg(3,filename_meano) + call getarg(3,filename_meano) ! new mean analysis ! and put in this file. - call getarg(4,filenameout) + call getarg(4,filenameout) ! new increment of analysis ! how many ensemble members to process call getarg(5,charnin) read(charnin,'(i4)') nanals ! option for increment, read in ens mean guess - call getarg(6,filename_meang) + call getarg(6,filename_meang) ! background ens mean fcst if (mype==0) then @@ -294,6 +294,9 @@ program recentersigp end select values_3d(:,:,:) = zero do j=1,latb +! updated member increment = original member increment - background ens mean - original +! mean increment + gsi control analysis +! = original member increment + gsi control analysis - enkf mean analysis values_3d(:,j,:) = values_3d_i(:,j,:) - values_3d_mb(:,latb-j+1,:) - values_3d_mi(:,j,:) + values_3d_anl(:,latb-j+1,:) end do if (should_zero_increments_for(trim(dseti%variables(nvar)%name))) values_3d = zero diff --git a/util/netcdf_io/calc_analysis.fd/init_calc_analysis.f90 b/util/netcdf_io/calc_analysis.fd/init_calc_analysis.f90 index 2d7c922da7..bd65f08849 100644 --- a/util/netcdf_io/calc_analysis.fd/init_calc_analysis.f90 +++ b/util/netcdf_io/calc_analysis.fd/init_calc_analysis.f90 @@ -41,16 +41,22 @@ subroutine read_nml stop 99 end if + if (fhr > 0) then write(hrstr,'(I0.2)') fhr anal_file = trim(adjustl(datapath)) // '/' // trim(adjustl(analysis_filename)) // '.' // hrstr fcst_file = trim(adjustl(datapath)) // '/' // trim(adjustl(firstguess_filename)) // '.' // hrstr incr_file = trim(adjustl(datapath)) // '/' // trim(adjustl(increment_filename)) // '.' // hrstr + else + anal_file = trim(adjustl(datapath)) // '/' // trim(adjustl(analysis_filename)) + fcst_file = trim(adjustl(datapath)) // '/' // trim(adjustl(firstguess_filename)) + incr_file = trim(adjustl(datapath)) // '/' // trim(adjustl(increment_filename)) + endif if (mype == 0) then write(6,*) 'Analysis File = ', trim(anal_file) write(6,*) 'First Guess File = ', trim(fcst_file) write(6,*) 'Increment File = ', trim(incr_file) - write(6,*) 'Forecast Hour = ', fhr + if (fhr > 0) write(6,*) 'Forecast Hour = ', fhr write(6,*) 'Number of PEs = ', npes write(6,*) 'input guess file and increment file should be in netCDF format' if (use_nemsio_anl) then