From 0b26ca0da26a2763007a22756812e3482bd9289e Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 7 Jan 2022 16:12:33 +0000 Subject: [PATCH 1/7] Github issue #37: fix for broken generate_ens=T option --- src/gsi/bkerror.f90 | 5 +- src/gsi/bkgcov.f90 | 37 +- src/gsi/control2model.f90 | 2 +- src/gsi/control2model_ad.f90 | 2 +- src/gsi/cplr_gfs_ensmod.f90 | 10 +- src/gsi/glbsoi.f90 | 6 +- src/gsi/gsimod.F90 | 9 +- src/gsi/hybrid_ensemble_isotropic.F90 | 100 ++++-- src/gsi/hybrid_ensemble_parameters.f90 | 9 +- src/gsi/netcdfgfs_io.f90 | 446 +++++++++++++++++++++++++ src/gsi/prewgt.f90 | 3 +- 11 files changed, 563 insertions(+), 66 deletions(-) diff --git a/src/gsi/bkerror.f90 b/src/gsi/bkerror.f90 index b3a014069..c7ec05b75 100644 --- a/src/gsi/bkerror.f90 +++ b/src/gsi/bkerror.f90 @@ -62,6 +62,9 @@ subroutine bkerror(grady) use control_vectors, only: mvars,nrf,nrf_var,nrf_3d use timermod, only: timer_ini,timer_fnl use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub + use general_commvars_mod, only: s2g_raf + use general_commvars_mod, only: s2g_cv use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult use hybrid_ensemble_parameters, only: l_hyb_ens implicit none @@ -126,7 +129,7 @@ subroutine bkerror(grady) endif ! Apply variances, as well as vertical & horizontal parts of background error - call bkgcov(mbundle) + call bkgcov(s2g_raf,mbundle) ! The following lines test that indeed proper application of cgkcov ! reproduces results of bkgcov - left as comments (please do not remove diff --git a/src/gsi/bkgcov.f90 b/src/gsi/bkgcov.f90 index 77733368d..10d643ae9 100644 --- a/src/gsi/bkgcov.f90 +++ b/src/gsi/bkgcov.f90 @@ -1,4 +1,4 @@ -subroutine bkgcov(cstate) +subroutine bkgcov(grd,cstate) !$$$ subprogram documentation block ! . . . . ! subprogram: bkgcov perform hor & vert of background error @@ -42,18 +42,19 @@ subroutine bkgcov(cstate) use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use general_sub2grid_mod, only: general_sub2grid,general_grid2sub - use general_commvars_mod, only: s2g_raf + use general_sub2grid_mod, only: sub2grid_info implicit none + type(sub2grid_info), intent(in) :: grd ! Passed Variables type(gsi_bundle),intent(inout) :: cstate ! Local Variables integer(i_kind) n,n3d,istatus,nlevs - real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork + real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork real(r_kind),pointer,dimension(:,:,:):: ptr3d=>NULL() - nlevs=s2g_raf%nlevs_loc + nlevs=grd%nlevs_loc n3d=cstate%n3d ! Multiply by background error variances, and break up skin temp @@ -68,13 +69,13 @@ subroutine bkgcov(cstate) end do ! Convert from subdomain to full horizontal field distributed among processors - call general_sub2grid(s2g_raf,cstate%values,hwork) + call general_sub2grid(grd,cstate%values,hwork) ! Apply horizontal smoother for number of horizontal scales call smoothrf(hwork,nlevs) ! Put back onto subdomains - call general_grid2sub(s2g_raf,hwork,cstate%values) + call general_grid2sub(grd,hwork,cstate%values) ! Apply vertical smoother !$omp parallel do schedule(dynamic,1) private(n,ptr3d,istatus) @@ -90,7 +91,7 @@ subroutine bkgcov(cstate) return end subroutine bkgcov ! ----------------------------------------------------------------------------- -subroutine ckgcov(z,cstate,nval_lenz) +subroutine ckgcov(grd,z,cstate,nval_lenz) !$$$ subprogram documentation block ! . . . . ! subprogram: ckgcov sqrt of bkgcov @@ -132,29 +133,29 @@ subroutine ckgcov(z,cstate,nval_lenz) use gridmod, only: nlat,nlon use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer - use general_sub2grid_mod, only: general_grid2sub - use general_commvars_mod, only: s2g_raf + use general_sub2grid_mod, only: general_grid2sub, sub2grid_info use hybrid_ensemble_parameters, only: l_hyb_ens use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult implicit none ! Passed Variables integer(i_kind) ,intent(in ) :: nval_lenz + type(sub2grid_info), intent(in) :: grd type(gsi_bundle),intent(inout) :: cstate real(r_kind),dimension(nval_lenz),intent(in ) :: z ! Local Variables integer(i_kind) k,n3d,istatus,nlevs - real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork + real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork real(r_kind),dimension(:,:,:),pointer:: ptr3d=>NULL() - nlevs=s2g_raf%nlevs_loc + nlevs=grd%nlevs_loc ! Apply horizontal smoother for number of horizontal scales call sqrt_smoothrf(z,hwork,nlevs) ! Put back onto subdomains - call general_grid2sub(s2g_raf,hwork,cstate%values) + call general_grid2sub(grd,hwork,cstate%values) ! Apply vertical smoother n3d=cstate%n3d @@ -174,7 +175,7 @@ subroutine ckgcov(z,cstate,nval_lenz) return end subroutine ckgcov ! ----------------------------------------------------------------------------- -subroutine ckgcov_ad(z,cstate,nval_lenz) +subroutine ckgcov_ad(grd,z,cstate,nval_lenz) !$$$ subprogram documentation block ! . . . . ! subprogram: ckgcov_ad adjoint of ckgcov @@ -217,11 +218,11 @@ subroutine ckgcov_ad(z,cstate,nval_lenz) use gridmod, only: nlat,nlon use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer - use general_sub2grid_mod, only: general_sub2grid - use general_commvars_mod, only: s2g_raf + use general_sub2grid_mod, only: general_sub2grid,sub2grid_info use hybrid_ensemble_parameters, only: l_hyb_ens use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult implicit none + type(sub2grid_info), intent(in) :: grd ! Passed Variables integer(i_kind) ,intent(in ) :: nval_lenz @@ -230,10 +231,10 @@ subroutine ckgcov_ad(z,cstate,nval_lenz) ! Local Variables integer(i_kind) k,n3d,istatus,nlevs - real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork + real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork real(r_kind),dimension(:,:,:),pointer:: ptr3d=>NULL() - nlevs=s2g_raf%nlevs_loc + nlevs=grd%nlevs_loc ! Apply static betas if(l_hyb_ens) call sqrt_beta_s_mult(cstate) @@ -251,7 +252,7 @@ subroutine ckgcov_ad(z,cstate,nval_lenz) end do ! Convert from subdomain to full horizontal field distributed among processors - call general_sub2grid(s2g_raf,cstate%values,hwork) + call general_sub2grid(grd,cstate%values,hwork) ! Apply horizontal smoother for number of horizontal scales call sqrt_smoothrf_ad(z,hwork,nlevs) diff --git a/src/gsi/control2model.f90 b/src/gsi/control2model.f90 index ec628afe4..b355ade4c 100644 --- a/src/gsi/control2model.f90 +++ b/src/gsi/control2model.f90 @@ -173,7 +173,7 @@ subroutine control2model(xhat,sval,bval) ! Apply sqrt of variance, as well as vertical & horizontal parts of background ! error - call ckgcov(xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) + call ckgcov(grid,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) ! Get pointers to required state variables call gsi_bundlegetpointer (sval(jj),'u' ,sv_u, istatus) diff --git a/src/gsi/control2model_ad.f90 b/src/gsi/control2model_ad.f90 index 639bf87ff..aed25adfa 100644 --- a/src/gsi/control2model_ad.f90 +++ b/src/gsi/control2model_ad.f90 @@ -236,7 +236,7 @@ subroutine control2model_ad(rval,bval,grad) enddo ! Apply adjoint of sqrt-B - call ckgcov_ad(gradz,wbundle,nval_lenz) + call ckgcov_ad(grid,gradz,wbundle,nval_lenz) ! Clean up call gsi_bundledestroy(wbundle,istatus) diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index ac00db4c4..f061f5491 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -1371,6 +1371,7 @@ subroutine put_gfs_ens(this,grd,member,ntindex,pert,iret) use hybrid_ensemble_parameters, only: ensemble_path use hybrid_ensemble_parameters, only: sp_ens use gridmod, only: use_gfs_nemsio, use_gfs_ncio + use netcdfgfs_io, only: write_gfsncatm_pert implicit none @@ -1402,11 +1403,12 @@ subroutine put_gfs_ens(this,grd,member,ntindex,pert,iret) endif !call write_nemsatm(grd,...) else if ( use_gfs_ncio ) then - if ( mype == 0 ) then - write(6,*) 'write_gfsncatm is not adapted to write out perturbations yet' - iret = 999 - endif + !if ( mype == 0 ) then + ! write(6,*) 'write_gfsncatm is not adapted to write out perturbations yet' + ! iret = 999 + !endif !call write_gfsncatm(grd,...) + call write_gfsncatm_pert(grd,sp_ens,filename,mype_atm,pert,ntindex) else call general_write_gfsatm(grd,sp_ens,sp_ens,filename,mype_atm, & pert,ntindex,inithead,iret) diff --git a/src/gsi/glbsoi.f90 b/src/gsi/glbsoi.f90 index a210ba325..faa335628 100644 --- a/src/gsi/glbsoi.f90 +++ b/src/gsi/glbsoi.f90 @@ -150,7 +150,7 @@ subroutine glbsoi use zrnmi_mod, only: zrnmi_initialize use observermod, only: observer_init,observer_set,observer_finalize,ndata use timermod, only: timer_ini, timer_fnl - use hybrid_ensemble_parameters, only: l_hyb_ens,destroy_hybens_localization_parameters + use hybrid_ensemble_parameters, only: l_hyb_ens,destroy_hybens_localization_parameters,create_hybens_localization_parameters,generate_ens,write_generated_ens use hybrid_ensemble_isotropic, only: create_ensemble,load_ensemble,destroy_ensemble, & hybens_localization_setup,hybens_grid_setup use gfs_stratosphere, only: destroy_nmmb_vcoords,use_gfs_stratosphere @@ -278,6 +278,7 @@ subroutine glbsoi ! If l_hyb_ens is true, then read in ensemble perturbations if(l_hyb_ens) then + call create_hybens_localization_parameters call load_ensemble call hybens_localization_setup end if @@ -312,6 +313,9 @@ subroutine glbsoi if (lsensrecompute) jiterlast=jiterend if (l4dvar) jiterlast=jiterstart if (ladtest_obs) jiterlast=jiterstart + if (l_hyb_ens .and. generate_ens .and. write_generated_ens) then + jiterlast=0 ! skip analysis loop + endif ! Main outer analysis loop do jiter=jiterstart,jiterlast diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 1c8b73841..89cb8cedb 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -144,7 +144,7 @@ module gsimod use hybrid_ensemble_parameters,only : l_hyb_ens,uv_hyb_ens,aniso_a_en,generate_ens,& 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,& + full_ensemble,pseudo_hybens,pwgtflg,write_generated_ens,& beta_s0,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, & @@ -715,8 +715,8 @@ module gsimod idmodel,iwrtinc,lwrite4danl,nhr_anal,jiterstart,jiterend,lobserver,lanczosave,llancdone, & lferrscale,print_diag_pcg,tsensible,lread_obs_save,lread_obs_skip, & use_gfs_ozone,check_gfs_ozone_date,regional_ozone,lwrite_predterms,& - lwrite_peakwt,use_gfs_nemsio,use_gfs_ncio,sfcnst_comb,liauon,use_prepb_satwnd,l4densvar,ens_nstarthr,& - use_gfs_stratosphere,pblend0,pblend1,step_start,diag_precon,lrun_subdirs,& + lwrite_peakwt,use_gfs_nemsio,use_gfs_ncio,sfcnst_comb,liauon,use_prepb_satwnd,l4densvar,& + ens_nstarthr,use_gfs_stratosphere,pblend0,pblend1,step_start,diag_precon,lrun_subdirs,& use_sp_eqspace,lnested_loops,lsingleradob,thin4d,use_readin_anl_sfcmask,& luse_obsdiag,id_drifter,id_ship,verbose,print_obs_para,lsingleradar,singleradar,lnobalance, & missing_to_nopcp,minobrangedbz,minobrangedbz,maxobrangedbz,& @@ -1255,6 +1255,7 @@ module gsimod ! oz_univ_static- if true, decouple ozone from other variables and defaults to static B (ozone only) ! aniso_a_en - if true, then use anisotropic localization of hybrid ensemble control variable a_en. ! generate_ens - if true, then generate internal ensemble based on existing background error +! write_generated_ens - if true, then writed generated internal ensemble based on existing background error ! n_ens - number of ensemble members. ! nlon_ens - number of longitudes on ensemble grid (may be different from analysis grid nlon) ! nlat_ens - number of latitudes on ensemble grid (may be different from analysis grid nlat) @@ -1313,7 +1314,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,& - grid_ratio_ens, & + grid_ratio_ens,write_generated_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 e0a6e3fd1..23a55f027 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -1167,7 +1167,7 @@ subroutine load_ensemble use constants, only: zero,one use hybrid_ensemble_parameters, only: n_ens,generate_ens,grd_ens,grd_anl,ntlevs_ens, & pseudo_hybens,regional_ensemble_option,& - i_en_perts_io + i_en_perts_io,write_generated_ens use hybrid_ensemble_parameters, only: nelen,en_perts,ps_bar use gsi_enscouplermod, only: gsi_enscoupler_put_gsi_ens use mpimod, only: mype @@ -1175,27 +1175,38 @@ subroutine load_ensemble use get_wrf_mass_ensperts_mod, only: get_wrf_mass_ensperts_class use get_fv3_regional_ensperts_mod, only: get_fv3_regional_ensperts_class use get_wrf_nmm_ensperts_mod, only: get_wrf_nmm_ensperts_class - use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens + use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens use mpimod, only: mpi_comm_world + use gsi_bundlemod, only: gsi_bundlegetpointer implicit none - type(get_pseudo_ensperts_class) :: pseudo_enspert - type(get_wrf_mass_ensperts_class) :: wrf_mass_enspert - type(get_wrf_nmm_ensperts_class) :: wrf_nmm_enspert - type(get_fv3_regional_ensperts_class) :: fv3_regional_enspert + type(get_pseudo_ensperts_class) :: pseudo_enspert + type(get_wrf_mass_ensperts_class) :: wrf_mass_enspert + type(get_wrf_nmm_ensperts_class) :: wrf_nmm_enspert + type(get_fv3_regional_ensperts_class) :: fv3_regional_enspert type(gsi_bundle),allocatable:: en_bar(:) type(gsi_bundle):: bundle_anl,bundle_ens type(gsi_grid) :: grid_anl,grid_ens - integer(i_kind) i,j,n,ii,m + integer(i_kind) i,j,n,ii,m,m1,m2,k integer(i_kind) istatus real(r_kind),allocatable:: seed(:,:) real(r_kind),pointer,dimension(:,:) :: cv_ps=>NULL() + real(r_kind),pointer,dimension(:,:,:) :: cv_q=>NULL() + real(r_kind),pointer,dimension(:,:,:) :: cv_t=>NULL() real(r_kind) sig_norm,bar_norm character(len=*),parameter::myname_=trim(myname)//'*load_ensemble' + real(r_kind),allocatable,dimension(:,:,:) :: work,work1 ! create simple regular grid + m1=1; m2 = ntlevs_ens + if (write_generated_ens) then + ! only write out one sample + m1=ntlevs_ens/2 + 1 + m2=ntlevs_ens/2 + 1 + endif + if(generate_ens) then @@ -1204,7 +1215,7 @@ subroutine load_ensemble allocate(en_bar(ntlevs_ens)) - do m=1,ntlevs_ens + do m=m1,m2 call gsi_bundlecreate(en_bar(m),grid_ens,'ensemble',istatus, & names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind) enddo @@ -1221,7 +1232,7 @@ subroutine load_ensemble allocate(seed(nval2f,nscl)) seed=-one - do m=1,ntlevs_ens + do m=m1,m2 en_bar(m)%values=zero enddo @@ -1240,7 +1251,7 @@ subroutine load_ensemble call stop2(999) endif - do m=1,ntlevs_ens + do m=m1,m2 do n=1,n_ens call generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed) do ii=1,nelen @@ -1258,6 +1269,50 @@ subroutine load_ensemble enddo enddo +! remove mean, which is locally significantly non-zero, due to sample size. +! with real ensembles, the mean of the actual sample will be removed. + + do m=m1,m2 + do n=1,n_ens + do ii=1,nelen + en_perts(n,m)%valuesr4(ii)=(en_perts(n,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm + enddo + enddo + call gsi_bundledestroy(en_bar(m),istatus) + if(istatus/=0) then + write(6,*)trim(myname_),': trouble destroying en_bar bundle' + call stop2(999) + end if + enddo + + if (write_generated_ens) then + allocate ( work(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) ) + allocate ( work1(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig+1) ) + + do m=m1,m2 + do n=1,n_ens + if (mype == 0) print *,'write out ens perturbation',n + do ii=1,nelen + bundle_ens%values(ii) = en_perts(n,m)%valuesr4(ii)/sig_norm + enddo + call gsi_bundlegetpointer (bundle_ens,'q' ,cv_q ,istatus) + work(:,:,:) = cv_q(:,:,:) ! make a copy + call gsi_bundlegetpointer (bundle_ens,'t' ,cv_t ,istatus) + call gsi_bundlegetpointer (bundle_ens,'ps',cv_ps ,istatus) + ! Get 3d pressure + call getprs_tl(cv_ps,cv_t,work1) + ! Convert input normalized RH to q + call normal_rh_to_q(work,cv_t,work1,cv_q) + call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,bundle_ens,istatus) + if(istatus/=0) then + write(6,*)trim(myname_),': trouble writing perts' + call stop2(999) + endif + enddo + enddo + deallocate(work,work1) + end if + ! do some cleanning call gsi_bundledestroy(bundle_anl,istatus) if(istatus/=0) then @@ -1269,27 +1324,6 @@ subroutine load_ensemble write(6,*)trim(myname_),': trouble destroying bundle_ens bundle' call stop2(999) endif -! remove mean, which is locally significantly non-zero, due to sample size. -! with real ensembles, the mean of the actual sample will be removed. - - do m=1,ntlevs_ens - do n=1,n_ens - do ii=1,nelen - en_perts(n,m)%valuesr4(ii)=(en_perts(n,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm - enddo - call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,en_perts(n,m),istatus) - if(istatus/=0) then - write(6,*)trim(myname_),': trouble writing perts' - call stop2(999) - endif - enddo - - call gsi_bundledestroy(en_bar(m),istatus) - if(istatus/=0) then - write(6,*)trim(myname_),': trouble destroying en_bar bundle' - call stop2(999) - end if - enddo deallocate(en_bar) deallocate(seed) @@ -1499,7 +1533,7 @@ subroutine generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed) ! temporarily redefine nval_lenz nval_lenz_save=nval_lenz nval_lenz=nval2f*nnnn1o*nscl - call ckgcov(z,bundle_anl,nval_lenz) + call ckgcov(grd_anl,z,bundle_anl,nval_lenz) ! restore nval_lenz nval_lenz=nval_lenz_save @@ -4017,7 +4051,7 @@ subroutine hybens_localization_setup if(verbose .and. mype == 0)print_verbose=.true. ! Allocate - call create_hybens_localization_parameters + !call create_hybens_localization_parameters if ( readin_localization .or. readin_beta ) then ! read info from file diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90 index 9ea28c9ee..1748b8dcc 100644 --- a/src/gsi/hybrid_ensemble_parameters.f90 +++ b/src/gsi/hybrid_ensemble_parameters.f90 @@ -92,6 +92,9 @@ module hybrid_ensemble_parameters ! generate_ens: if .true., generate ensemble perturbations internally as random samples of background B. ! (used primarily for testing/debugging) ! if .false., read external ensemble perturbations +! write_generate_ens: if .true., write out ens perturbations generated internally as +! as samples of static B and exit. +! aniso_a_en: if .true., then allow anisotropic localization correlation (not active yet) ! aniso_a_en: if .true., then allow anisotropic localization correlation (not active yet) ! uv_hyb_ens: if .true., then ensemble perturbation wind stored as u,v ! if .false., ensemble perturbation wind stored as psi,chi. @@ -159,6 +162,7 @@ module hybrid_ensemble_parameters ! def aniso_a_en - if true, then use anisotropic rf for localization ! def generate_ens - if true, then create ensemble members internally ! using sqrt of static background error acting on N(0,1) random vectors +! def write_generated_ens - if true, write out generated ensemble members ! def n_ens - number of ensemble members ! def nlon_ens - number of longitudes to use for ensemble members and ensemble control vector ! def nlat_ens - number of latitudes to use for ensemble members and ensemble control vector @@ -254,7 +258,7 @@ module hybrid_ensemble_parameters destroy_hybens_localization_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 + s_ens_h,oz_univ_static,vvlocal,write_generated_ens public :: uv_hyb_ens,q_hyb_ens,s_ens_v,beta_s0,aniso_a_en,s_ens_hv,s_ens_vv public :: readin_beta,beta_s,beta_e public :: readin_localization @@ -294,7 +298,7 @@ module hybrid_ensemble_parameters logical l_hyb_ens,uv_hyb_ens,q_hyb_ens,oz_univ_static,sst_staticB logical aniso_a_en logical full_ensemble,pwgtflg - logical generate_ens + logical generate_ens, write_generated_ens logical dual_res logical pseudo_hybens logical merge_two_grid_ensperts @@ -383,6 +387,7 @@ subroutine init_hybrid_ensemble_parameters sst_staticB=.true. aniso_a_en=.false. generate_ens=.true. + write_generated_ens=.false. pseudo_hybens=.false. merge_two_grid_ensperts=.false. regional_ensemble_option=0 diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index dd3f991da..816ede246 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -52,6 +52,7 @@ module netcdfgfs_io public intrp22 public tran_gfsncsfc public write_gfsncatm + public write_gfsncatm_pert interface read_gfsnc module procedure read_ @@ -90,6 +91,10 @@ module netcdfgfs_io module procedure write_atm_ end interface + interface write_gfsncatm_pert + module procedure write_atm_pert_ + end interface + character(len=*),parameter::myname='netcdfgfs_io' contains @@ -2256,6 +2261,447 @@ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) end subroutine write_atm_ + subroutine write_atm_pert_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) + +!$$$ subprogram documentation block +! . . . +! subprogram: write_ncatm --- Gather, transform, and write out to netcdf +! +! prgmmr: whitaker org: oar date: 2019-10-03 +! +! abstract: This routine gathers fields needed for the GSI analysis +! file from subdomains and then transforms the fields from +! analysis grid to model guess grid, then written to an +! netcdf atmospheric analysis file. +! +! program history log: +! 2019-10-03 whitaker initial version +! +! input argument list: +! filename - file to open and write to +! mype_out - mpi task to write output file +! gfs_bundle - bundle containing fields on subdomains +! ibin - time bin +! +! output argument list: +! +! attributes: +! language: f90 +! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP +! +!$$$ end documentation block + +! !USES: + use kinds, only: r_kind,i_kind + + use constants, only: r1000,fv,one,zero,qcmin,r0_05,t0c + + use mpimod, only: mpi_rtype + use mpimod, only: mpi_comm_world + use mpimod, only: ierror + use mpimod, only: mype + + use guess_grids, only: ifilesig + use guess_grids, only: load_geop_hgt,geop_hgti,ges_geopi + + use gridmod, only: ntracer + use gridmod, only: ncloud + use gridmod, only: strip,jcap_b,bk5 + + use general_commvars_mod, only: load_grid,fill2_ns,filluv2_ns + use general_specmod, only: spec_vars + + use obsmod, only: iadate + + use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl + use general_sub2grid_mod, only: sub2grid_info + use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid + use constants, only: two,pi,half,deg2rad + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use cloud_efr_mod, only: cloud_calc_gfs + + use netcdf, only: nf90_max_name + use module_fv3gfs_ncio, only: open_dataset, close_dataset, Dimension, Dataset,& + read_attribute, write_attribute,get_dim, create_dataset, write_vardata, read_vardata,& + get_idate_from_time_units,quantize_data,get_time_units_from_idate,has_attr,has_var + use ncepnems_io, only: error_msg + + implicit none + +! !INPUT PARAMETERS: + + type(sub2grid_info), intent(in) :: grd + type(spec_vars), intent(in) :: sp_a + character(len=24), intent(in) :: filename ! file to open and write to + integer(i_kind), intent(in) :: mype_out ! mpi task to write output file + type(gsi_bundle), intent(in) :: gfs_bundle + integer(i_kind), intent(in) :: ibin ! time bin + +!------------------------------------------------------------------------- + + real(r_kind),parameter:: r0_001 = 0.001_r_kind + character(6):: fname_ges + character(len=120) :: my_name = 'WRITE_GFSNCATM' + character(len=1) :: null = ' ' + integer(i_kind),dimension(6):: idate,jdate + integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk, kr, nbits + integer(i_kind) :: iret, lonb, latb, levs, istatus + integer(i_kind) :: nfhour + integer(i_kind) :: istop = 104 + integer(i_kind),dimension(5):: mydate + integer(i_kind),dimension(8) :: ida,jda + real(r_kind),dimension(5) :: fha + real(r_kind), allocatable, dimension(:) :: fhour + real(r_kind),allocatable,dimension(:) :: rlats_tmp,rlons_tmp + + real(r_kind),pointer,dimension(:,:) :: sub_ps + real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v,sub_t + real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz,sub_cwmr + + + real(r_kind),dimension(grd%lat1*grd%lon1) :: psm + real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: tsm, usm, vsm + real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: qsm, ozsm + real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: cwsm + real(r_kind),dimension(max(grd%iglobal,grd%itotsub)) :: work1,work2 + real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid + real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons + real(4), allocatable, dimension(:,:) :: values_2d,values_2d_tmp + real(4), allocatable, dimension(:,:,:) :: values_3d,values_3d_tmp,ug3d,vg3d + real(4) compress_err + type(Dataset) :: atmges, atmanl + type(Dimension) :: ncdim + character(len=nf90_max_name) :: time_units + + logical diff_res,eqspace,nocompress + logical,dimension(1) :: vector + type(egrid2agrid_parm) :: p_low,p_high + +!************************************************************************* +! Initialize local variables + mm1=mype+1 + nlatm2=grd%nlat-2 + diff_res=.false. + nocompress = .true. + + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting ps from bundle' + call gsi_bundlegetpointer(gfs_bundle,'sf', sub_u, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting sf from bundle' + call gsi_bundlegetpointer(gfs_bundle,'vp', sub_v, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting vp from bundle' + call gsi_bundlegetpointer(gfs_bundle,'t ', sub_t, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting t from bundle' + call gsi_bundlegetpointer(gfs_bundle,'q', sub_q, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting q from bundle' + call gsi_bundlegetpointer(gfs_bundle,'oz', sub_oz, iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting oz from bundle' + call gsi_bundlegetpointer(gfs_bundle,'cw', sub_cwmr,iret); istatus=istatus+iret + if (iret /= 0 .and. mype == 0) print *,'error getting cw from bundle' + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'write_atm_: ERROR' + write(6,*) 'Missing some of the required fields' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + + if ( sp_a%jcap /= jcap_b ) then + if ( mype == 0 ) write(6, & + '('' dual resolution for nems sp_a%jcap,jcap_b = '',2i6)') & + sp_a%jcap,jcap_b + diff_res = .true. + call stop2(999) + endif + + + ! Single task writes analysis data to analysis file + if ( mype == mype_out ) then + write(fname_ges,'(''sigf'',i2.2)') ifilesig(ibin) + + ! open the netCDF file + atmges = open_dataset(fname_ges,errcode=iret) + if ( iret /= 0 ) call error_msg(trim(my_name),null,null,'open',istop,iret) + ! get dimension sizes + ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len + ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len + ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len + + if ( levs /= grd%nsig ) then + write(6,*) trim(my_name),': problem in data dimension background levs = ',levs,' nsig = ',grd%nsig + call stop2(103) + endif + + ! get time information + idate = get_idate_from_time_units(atmges) + call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later + ! depends on model changes from Jeff Whitaker + nfhour = fhour(1) + + atmanl = create_dataset(filename, atmges, & + copy_vardata=.false., nocompress=nocompress, errcode=iret) + if ( iret /= 0 ) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret) + + ! Update time information (with ibdate) and write it to analysis file + mydate=ibdate + fha(:)=zero ; ida=0; jda=0 + fha(2)=real(nhr_obsbin*(ibin-1)) ! relative time interval in hours + ida(1)=mydate(1) ! year + ida(2)=mydate(2) ! month + ida(3)=mydate(3) ! day + ida(4)=0 ! time zone + ida(5)=mydate(4) ! hour + + ! Move date-time forward by nhr_assimilation hours + call w3movdat(fha,ida,jda) + + jdate(1) = jda(1) ! analysis year + jdate(2) = jda(2) ! analysis month + jdate(3) = jda(3) ! analysis day + jdate(4) = jda(5) ! analysis hour + jdate(5) = iadate(5) ! analysis minute + jdate(6) = 0 ! analysis second + fhour = zero + + call write_vardata(atmanl, 'time', fhour) + time_units = get_time_units_from_idate(jdate) + call write_attribute(atmanl, 'units', time_units, 'time') + + ! Allocate structure arrays to hold data + allocate(ug3d(lonb,latb,levs),vg3d(lonb,latb,levs),values_3d_tmp(lonb,latb,levs),values_2d_tmp(lonb,latb)) + allocate(values_3d(lonb,latb,levs)) + allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb)) + call read_vardata(atmges, 'grid_xt', rlons_tmp, errcode=iret) + call read_vardata(atmges, 'grid_yt', rlats_tmp, errcode=iret) + do j=1,latb + rlats(latb+2-j)=deg2rad*rlats_tmp(j) + enddo + rlats(1)=-half*pi + rlats(latb+2)=half*pi + do j=1,lonb + rlons(j)=deg2rad*rlons_tmp(j) + enddo + deallocate(rlons_tmp, rlats_tmp) + do j=1,lonb + clons(j)=cos(rlons(j)) + slons(j)=sin(rlons(j)) + enddo + + nord_int=4 + eqspace=.false. + call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, & + latb+2,rlats,lonb,rlons,& + nord_int,p_low,.false.,eqspace=eqspace) + call g_create_egrid2agrid(latb+2,rlats,lonb,rlons, & + grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,& + nord_int,p_high,.false.,eqspace=eqspace) + + deallocate(rlats,rlons) + + endif ! if ( mype == mype_out ) + + ! Strip off boundary points from subdomains + call strip(sub_ps ,psm) + call strip(sub_t ,tsm ,grd%nsig) + call strip(sub_q ,qsm ,grd%nsig) + call strip(sub_oz ,ozsm ,grd%nsig) + call strip(sub_cwmr,cwsm ,grd%nsig) + call strip(sub_u ,usm ,grd%nsig) + call strip(sub_v ,vsm ,grd%nsig) + + ! Thermodynamic variable + ! The GSI analysis variable is virtual temperature (Tv). For NEMSIO + ! output we need the sensible temperature. + + ! Generate and write analysis fields + + ! Surface pressure + call mpi_gatherv(psm,grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype==mype_out) then + call load_grid(work1,grid) + values_2d = grid*r1000 + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'pressfc')) then + call read_attribute(atmges, 'nbits', nbits, 'pressfc') + values_2d_tmp = values_2d + call quantize_data(values_2d_tmp, values_2d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'pressfc') + endif + print *,'min/max pressfc pert',minval(values_2d),maxval(values_2d) + call write_vardata(atmanl,'pressfc',values_2d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pressfc','write',istop,iret) + endif + +! u, v + do k=1,grd%nsig + kr = grd%nsig-k+1 + call mpi_gatherv(usm(1,k),grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + call mpi_gatherv(vsm(1,k),grd%ijn(mm1),mpi_rtype,& + work2,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype==mype_out) then + call load_grid(work1,grid) + ug3d(:,:,kr) = grid + call load_grid(work2,grid) + vg3d(:,:,kr) = grid + endif ! mype_out + end do + ! Zonal wind + if (mype==mype_out) then + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'ugrd')) then + call read_attribute(atmges, 'nbits', nbits, 'ugrd') + values_3d_tmp = ug3d + call quantize_data(values_3d_tmp, ug3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'ugrd') + endif + print *,'min/max ugrd pert',minval(ug3d),maxval(ug3d) + call write_vardata(atmanl,'ugrd',ug3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','write',istop,iret) + ! Meridional wind + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'vgrd')) then + call read_attribute(atmges, 'nbits', nbits, 'vgrd') + values_3d_tmp = vg3d + call quantize_data(values_3d_tmp, vg3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'vgrd') + endif + print *,'min/max vgrd pert',minval(vg3d),maxval(vg3d) + call write_vardata(atmanl,'vgrd',vg3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','write',istop,iret) + endif + +! Thermodynamic variable + do k=1,grd%nsig + kr = grd%nsig-k+1 + call mpi_gatherv(tsm(1,k),grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype == mype_out) then + call load_grid(work1,grid) + values_3d(:,:,kr) = grid + endif + end do + if (mype==mype_out) then + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'tmp')) then + call read_attribute(atmges, 'nbits', nbits, 'tmp') + values_3d_tmp = values_3d + call quantize_data(values_3d_tmp, values_3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'tmp') + endif + print *,'min/max tv pert',minval(values_3d),maxval(values_3d) + call write_vardata(atmanl,'tmp',values_3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret) + endif + +! Specific humidity + do k=1,grd%nsig + kr = grd%nsig-k+1 + call mpi_gatherv(qsm(1,k),grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype == mype_out) then + call load_grid(work1,grid) + values_3d(:,:,kr) = grid + endif + end do + if (mype==mype_out) then + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'spfh')) then + call read_attribute(atmges, 'nbits', nbits, 'spfh') + values_3d_tmp = values_3d + call quantize_data(values_3d_tmp, values_3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'spfh') + endif + print *,'min/max spfh pert',minval(values_3d),maxval(values_3d) + call write_vardata(atmanl,'spfh',values_3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','write',istop,iret) + endif + +! Ozone + do k=1,grd%nsig + kr = grd%nsig-k+1 + call mpi_gatherv(ozsm(1,k),grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype == mype_out) then + call load_grid(work1,grid) + values_3d(:,:,kr) = grid + endif + end do + if (mype==mype_out) then + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'o3mr')) then + call read_attribute(atmges, 'nbits', nbits, 'o3mr') + values_3d_tmp = values_3d + call quantize_data(values_3d_tmp, values_3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'o3mr') + endif + print *,'min/max o3mr pert',minval(values_3d),maxval(values_3d) + call write_vardata(atmanl,'o3mr',values_3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','write',istop,iret) + endif + +! Cloud condensate mixing ratio + if (ntracer>2 .or. ncloud>=1) then + + do k=1,grd%nsig + kr = grd%nsig-k+1 + call mpi_gatherv(cwsm(1,k),grd%ijn(mm1),mpi_rtype,& + work1,grd%ijn,grd%displs_g,mpi_rtype,& + mype_out,mpi_comm_world,ierror) + if (mype == mype_out) then + call load_grid(work1,grid) + values_3d(:,:,kr) = grid + endif !mype == mype_out + end do + if (mype==mype_out) then + if (.not. nocompress .and. has_attr(atmges, 'nbits', 'clwmr')) then + call read_attribute(atmges, 'nbits', nbits, 'clwmr') + values_3d_tmp = values_3d + call quantize_data(values_3d_tmp, values_3d, nbits, compress_err) + call write_attribute(atmanl,& + 'max_abs_compression_error',compress_err,'clwmr') + endif + print *,'min/max clwmr pert',minval(values_3d),maxval(values_3d) + call write_vardata(atmanl,'clwmr',values_3d,errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret) + endif + endif !ntracer + +! +! Deallocate local array +! + if (mype==mype_out) then + call close_dataset(atmges, errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret) + + call close_dataset(atmanl, errcode=iret) + if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret) +! +! Deallocate local array +! + if (allocated(values_2d)) deallocate(values_2d) + if (allocated(values_3d)) deallocate(values_3d) + if (allocated(values_2d_tmp)) deallocate(values_2d_tmp) + if (allocated(values_3d_tmp)) deallocate(values_3d_tmp) + if (allocated(ug3d)) deallocate(ug3d) + if (allocated(vg3d)) deallocate(vg3d) +! + write(6,'(a,'': atm pert written for lonb,latb,levs= '',3i6,'',valid hour= '',f4.1,'',idate= '',4i5)') & + trim(my_name),lonb,latb,levs,fhour(1),jdate(1:4) + endif + + end subroutine write_atm_pert_ + subroutine write_sfc_ (filename,mype_sfc,dsfct) !$$$ subprogram documentation block diff --git a/src/gsi/prewgt.f90 b/src/gsi/prewgt.f90 index 38d61050f..a73f974b5 100644 --- a/src/gsi/prewgt.f90 +++ b/src/gsi/prewgt.f90 @@ -487,7 +487,8 @@ subroutine prewgt(mype) end do ! Special case of dssv for qoption=2 and cw - if (qoption==2) call compute_qvar3d + !if (qoption==2) call compute_qvar3d + call compute_qvar3d !!!$omp parallel do schedule(dynamic,1) private(i,n,j,jx,ix,loc) do n=1,nc2d From 03cd899177c078f9d518229683d0eeb172be48e6 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 13 Jan 2022 14:36:49 +0000 Subject: [PATCH 2/7] fix segfault in regression tests --- src/gsi/control2model.f90 | 4 ++-- src/gsi/control2model_ad.f90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/gsi/control2model.f90 b/src/gsi/control2model.f90 index b355ade4c..e5e4272aa 100644 --- a/src/gsi/control2model.f90 +++ b/src/gsi/control2model.f90 @@ -48,7 +48,7 @@ subroutine control2model(xhat,sval,bval) use control_vectors, only: nc2d,mvars use bias_predictors, only: predictors use gsi_4dvar, only: nsubwin, l4dvar, lsqrtb -use gridmod, only: lat2,lon2,nsig +use gridmod, only: grd_a,lat2,lon2,nsig use jfunc, only: nsclen,npclen,ntclen use berror, only: varprd,fpsproj,fut2ps use balmod, only: balance @@ -173,7 +173,7 @@ subroutine control2model(xhat,sval,bval) ! Apply sqrt of variance, as well as vertical & horizontal parts of background ! error - call ckgcov(grid,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) + call ckgcov(grd_a,,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) ! Get pointers to required state variables call gsi_bundlegetpointer (sval(jj),'u' ,sv_u, istatus) diff --git a/src/gsi/control2model_ad.f90 b/src/gsi/control2model_ad.f90 index aed25adfa..3f85c3c37 100644 --- a/src/gsi/control2model_ad.f90 +++ b/src/gsi/control2model_ad.f90 @@ -40,7 +40,7 @@ subroutine control2model_ad(rval,bval,grad) use control_vectors, only: nc2d,mvars use bias_predictors, only: predictors use gsi_4dvar, only: nsubwin, lsqrtb -use gridmod, only: lat2,lon2,nsig +use gridmod, only: grd_a,lat2,lon2,nsig use berror, only: varprd,fpsproj,fut2ps use balmod, only: tbalance use jfunc, only: nsclen,npclen,ntclen,nval_lenz @@ -236,7 +236,7 @@ subroutine control2model_ad(rval,bval,grad) enddo ! Apply adjoint of sqrt-B - call ckgcov_ad(grid,gradz,wbundle,nval_lenz) + call ckgcov_ad(grd_a,gradz,wbundle,nval_lenz) ! Clean up call gsi_bundledestroy(wbundle,istatus) From b62bfb600b32cd9f4853429765cd0fe94e7e3c64 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 16 Nov 2022 19:30:45 +0000 Subject: [PATCH 3/7] fix module name --- src/gsi/netcdfgfs_io.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 02ad9dc38..c23e251a5 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -2393,7 +2393,7 @@ subroutine write_atm_pert_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) use cloud_efr_mod, only: cloud_calc_gfs use netcdf, only: nf90_max_name - use module_fv3gfs_ncio, only: open_dataset, close_dataset, Dimension, Dataset,& + use module_ncio, only: open_dataset, close_dataset, Dimension, Dataset,& read_attribute, write_attribute,get_dim, create_dataset, write_vardata, read_vardata,& get_idate_from_time_units,quantize_data,get_time_units_from_idate,has_attr,has_var use ncepnems_io, only: error_msg From c572799527224db7f2d389d95370c85264869292 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 16 Nov 2022 19:32:55 +0000 Subject: [PATCH 4/7] remove commented out call --- src/gsi/hybrid_ensemble_isotropic.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index d66351d54..30ff7c4bc 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -4055,9 +4055,6 @@ subroutine hybens_localization_setup print_verbose=.false. .and. mype == 0 if(verbose .and. mype == 0)print_verbose=.true. - ! Allocate - !call create_hybens_localization_parameters - if ( readin_localization .or. readin_beta ) then ! read info from file inquire(file=trim(fname),exist=lexist) From ccee50b9c6dac590014898ab9bd900bd6003bba8 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 16 Nov 2022 19:42:18 +0000 Subject: [PATCH 5/7] fix typo --- src/gsi/control2model.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/control2model.f90 b/src/gsi/control2model.f90 index e5e4272aa..60c345bd0 100644 --- a/src/gsi/control2model.f90 +++ b/src/gsi/control2model.f90 @@ -173,7 +173,7 @@ subroutine control2model(xhat,sval,bval) ! Apply sqrt of variance, as well as vertical & horizontal parts of background ! error - call ckgcov(grd_a,,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) + call ckgcov(grd_a,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:))) ! Get pointers to required state variables call gsi_bundlegetpointer (sval(jj),'u' ,sv_u, istatus) From e298e439ac77cfb2444fbb3f966906018d8b6322 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 16 Nov 2022 19:46:23 +0000 Subject: [PATCH 6/7] remove commented out print --- src/gsi/cplr_gfs_ensmod.f90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index 97ef74a52..828726f1d 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -1451,11 +1451,6 @@ subroutine put_gfs_ens(this,grd,member,ntindex,pert,iret) endif !call write_nemsatm(grd,...) else if ( use_gfs_ncio ) then - !if ( mype == 0 ) then - ! write(6,*) 'write_gfsncatm is not adapted to write out perturbations yet' - ! iret = 999 - !endif - !call write_gfsncatm(grd,...) call write_gfsncatm_pert(grd,sp_ens,filename,mype_atm,pert,ntindex) else call general_write_gfsatm(grd,sp_ens,sp_ens,filename,mype_atm, & From 48bbaf2a50452ea9a89a4e6c94c1aed8cdc7308d Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 16 Nov 2022 19:48:48 +0000 Subject: [PATCH 7/7] set my_name to WRITE_GFSNCATM_PERT in write_atm_pert_ --- src/gsi/netcdfgfs_io.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index c23e251a5..55d8ed908 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -2413,7 +2413,7 @@ subroutine write_atm_pert_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) real(r_kind),parameter:: r0_001 = 0.001_r_kind character(6):: fname_ges - character(len=120) :: my_name = 'WRITE_GFSNCATM' + character(len=120) :: my_name = 'WRITE_GFSNCATM_PERT' character(len=1) :: null = ' ' integer(i_kind),dimension(6):: idate,jdate integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk, kr, nbits