diff --git a/src/gsibec/gsi/CMakeLists.txt b/src/gsibec/gsi/CMakeLists.txt index f56f77d..0fdca3e 100644 --- a/src/gsibec/gsi/CMakeLists.txt +++ b/src/gsibec/gsi/CMakeLists.txt @@ -101,7 +101,6 @@ m_berror_stats_reg.F90 rapidrefresh_cldsurf_mod.f90 chemmod.f90 radarZ_global_module.f90 -fv3_regional_interface.f90 gsi_rfv3io_mod.f90 mod_fv3_lola.f90 grdcrd.f90 diff --git a/src/gsibec/gsi/balmod.F90 b/src/gsibec/gsi/balmod.F90 index a25a64a..024f36e 100644 --- a/src/gsibec/gsi/balmod.F90 +++ b/src/gsibec/gsi/balmod.F90 @@ -921,6 +921,8 @@ subroutine locatelat_reg(mype) ! 2005-03-28 wu - replace mlath with mlat ! 2005-04-22 treadon - change berror file to 4-byte reals ! 2005-06-06 wu - setup f1 for balance projection (st->t) when fstat=.true. +! 2022-04-20 x.zhang - add switch (usenewgfsberror)for no need to convert +! the unit of clat for using global 127-L BE in regional DA ! ! input argument list: ! mype - mpi task id @@ -937,6 +939,7 @@ subroutine locatelat_reg(mype) use m_kinds, only: r_single use gridmod, only: nlon,nlat,lat2,lon2,istart,jstart,region_lat use constants, only: deg2rad,one + use m_berror_stats, only: usenewgfsberror implicit none ! Declare passed variables @@ -960,9 +963,16 @@ subroutine locatelat_reg(mype) allocate( clat_avn(mlat), clat_avn4(mlat) ) read(lunin)clat_avn4 close(lunin) - do i=1,mlat - clat_avn(i)=clat_avn4(i)*deg2rad - end do + if (usenewgfsberror) then +! 'The unit of clat from global BE is radian, do not need to convert' + do i=1,mlat + clat_avn(i)=clat_avn4(i) + end do + else + do i=1,mlat + clat_avn(i)=clat_avn4(i)*deg2rad + end do + end if deallocate(clat_avn4) diff --git a/src/gsibec/gsi/compute_derived.F90 b/src/gsibec/gsi/compute_derived.F90 index 62ee1d2..0e5aea1 100644 --- a/src/gsibec/gsi/compute_derived.F90 +++ b/src/gsibec/gsi/compute_derived.F90 @@ -137,6 +137,8 @@ subroutine compute_derived(mype,init_pass) use mpeu_util, only: getindex use mpeu_util, only: die, tell use gsi_io, only: verbose + use gen_qsat + implicit none @@ -191,7 +193,7 @@ subroutine compute_derived(mype,init_pass) iderivative=0 ice=.true. do ii=1,nfldsig - call genqsat(ges_qsat(1,1,1,ii),ges_tsen(1,1,1,ii),ges_prsl(1,1,1,ii),lat2,lon2, & + call genqsat(ges_qsat(:,:,:,ii),ges_tsen(:,:,:,ii),ges_prsl(:,:,:,ii),lat2,lon2, & nsig,ice,iderivative) enddo @@ -289,12 +291,13 @@ subroutine compute_derived(mype,init_pass) ! and for getting time derivatives of prognostic variables for ! time extrapolation and non-linear balance constraints. - do nt=1,nfldsig - call get_derivatives(gsi_metguess_bundle(nt),& - gsi_xderivative_bundle(nt), & - gsi_yderivative_bundle(nt)) - enddo - + if(.not. regional) then + do nt=1,nfldsig + call get_derivatives(gsi_metguess_bundle(nt),& + gsi_xderivative_bundle(nt), & + gsi_yderivative_bundle(nt)) + enddo + endif #ifdef USE_ALL_ORIGINAL if(.not. wrf_mass_regional .and. tendsflag)then #else @@ -438,13 +441,13 @@ subroutine compute_derived(mype,init_pass) end if ice=.true. - call genqsat(qsatg,ges_tsen(1,1,1,ntguessig),ges_prsl(1,1,1,ntguessig),lat2,lon2, & + call genqsat(qsatg,ges_tsen(:,:,:,ntguessig),ges_prsl(:,:,:,ntguessig),lat2,lon2, & nsig,ice,iderivative) ! Now load over nfldsig bins for limq (when nobs_bins /= zero) iderivative = 0 do ii=1,nfldsig - call genqsat(ges_qsat(1,1,1,ii),ges_tsen(1,1,1,ii),ges_prsl(1,1,1,ii),lat2,lon2, & + call genqsat(ges_qsat(:,:,:,ii),ges_tsen(:,:,:,ii),ges_prsl(:,:,:,ii),lat2,lon2, & nsig,ice,iderivative) end do endif diff --git a/src/gsibec/gsi/compute_qvar3d.F90 b/src/gsibec/gsi/compute_qvar3d.F90 index 12afcfa..bd7a897 100644 --- a/src/gsibec/gsi/compute_qvar3d.F90 +++ b/src/gsibec/gsi/compute_qvar3d.F90 @@ -55,6 +55,8 @@ subroutine compute_qvar3d use m_berror_stats, only: varq #endif /* USE_ALL_ORIGINAL */ + use gen_qsat + implicit none ! Declare local variables @@ -84,7 +86,7 @@ subroutine compute_qvar3d iderivative = 0 ice=.true. do it=1,nfldsig - call genqsat(ges_qsat(1,1,1,it),ges_tsen(1,1,1,it),ges_prsl(1,1,1,it),lat2,lon2, & + call genqsat(ges_qsat(:,:,:,it),ges_tsen(:,:,:,it),ges_prsl(:,:,:,it),lat2,lon2, & nsig,ice,iderivative) enddo @@ -127,7 +129,7 @@ subroutine compute_qvar3d iderivative = 2 end if ice=.true. - call genqsat(qsatg,ges_tsen(1,1,1,ntguessig),ges_prsl(1,1,1,ntguessig),lat2,lon2,nsig,ice,iderivative) + call genqsat(qsatg,ges_tsen(:,:,:,ntguessig),ges_prsl(:,:,:,ntguessig),lat2,lon2,nsig,ice,iderivative) if (qoption==2) then allocate(rhgues(lat2,lon2,nsig)) diff --git a/src/gsibec/gsi/control2state_ad.F90 b/src/gsibec/gsi/control2state_ad.F90 index 8ce6e71..7b4dbec 100644 --- a/src/gsibec/gsi/control2state_ad.F90 +++ b/src/gsibec/gsi/control2state_ad.F90 @@ -310,6 +310,7 @@ subroutine control2state_ad(rval,bval,grad) ! Adjoint of convert input normalized RH to q to add contribution of moisture ! to t, p , and normalized rh + if(regional .and. ls_prse) rv_prse = 0. if(do_normal_rh_to_q_ad) call normal_rh_to_q_ad(cv_rh,cv_t,rv_prse,rv_q) ! Adjoint to convert ps to 3-d pressure diff --git a/src/gsibec/gsi/fv3_regional_interface.f90 b/src/gsibec/gsi/fv3_regional_interface.f90 deleted file mode 100644 index 18c0a46..0000000 --- a/src/gsibec/gsi/fv3_regional_interface.f90 +++ /dev/null @@ -1,43 +0,0 @@ -subroutine convert_fv3_regional -!$$$ subprogram documentation block -! . . . . -! subprogram: convert_fv3_regional read single fv3 nest -! prgmmr: parrish org: np22 date: 2017-04-09 -! -! abstract: using routines from gsi_rfv3io_mod.f90 module to setup for -! reading tile of forecast fields from an fv3 forecast. -! NOTE: run on single processor, with information stored on unit lendian_out -! -!################################################################################# -!################################################################################# -! Use subroutine convert_nems_nmmb (in wrf_binary_interface.F90) as pattern. -!################################################################################# -!################################################################################# -! -! program history log: -! 2017-04-08 parrish -! 2018-02-16 wu - read in grid and time infor from fv3 files -! read directly from fv3 files and not writeout GSI internal file -! 2020-11-19 Lu & Wang - add time label it for fgat. POC: xuguang.wang@ou.edu -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - - use m_kinds, only: r_single,r_kind,i_kind - use gsi_rfv3io_mod, only: gsi_rfv3io_get_grid_specs - - implicit none - integer(i_kind) ierr - - -!!!!!!!!!!! get grid specs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - call gsi_rfv3io_get_grid_specs(ierr) - if(ierr/=0)then - write(6,*)' problem in convert_fv3_regional - get_grid_specs Status = ',ierr - call stop2 (555) - endif -end subroutine convert_fv3_regional - diff --git a/src/gsibec/gsi/genqsat.F90 b/src/gsibec/gsi/genqsat.F90 index 1f26af5..0d181d4 100644 --- a/src/gsibec/gsi/genqsat.F90 +++ b/src/gsibec/gsi/genqsat.F90 @@ -1,3 +1,5 @@ +module gen_qsat +contains subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) !$$$ subprogram documentation block ! . . . . @@ -62,8 +64,8 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) implicit none logical ,intent(in ) :: ice - real(r_kind),dimension(lat2,lon2,nsig),intent( out) :: qsat - real(r_kind),dimension(lat2,lon2,nsig),intent(in ) :: tsen,prsl + real(r_kind),dimension(:,:,:),intent( out) :: qsat + real(r_kind),dimension(:,:,:),intent(in ) :: tsen,prsl integer(i_kind) ,intent(in ) :: lat2,lon2,nsig,iderivative @@ -234,4 +236,4 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) end do return end subroutine genqsat - +end module gen_qsat diff --git a/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 b/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 index cdd7590..94e23c5 100644 --- a/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 +++ b/src/gsibec/gsi/get_gefs_ensperts_dualres.F90 @@ -75,6 +75,7 @@ subroutine get_gefs_ensperts_dualres (tau) #ifdef USE_ALL_ORIGINAL use m_revBens, only: revBens_ensmean_overwrite #endif /* USE_ALL_ORIGINAL */ + use gen_qsat implicit none integer(i_kind),intent(in) :: tau diff --git a/src/gsibec/gsi/gsi_rfv3io_mod.f90 b/src/gsibec/gsi/gsi_rfv3io_mod.f90 index 29d6dda..9cd3c15 100644 --- a/src/gsibec/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsibec/gsi/gsi_rfv3io_mod.f90 @@ -466,9 +466,13 @@ subroutine gsi_rfv3io_get_grid_specs(ierr) !!! get ak,bk - allocate(aeta1_ll(nsig),aeta2_ll(nsig)) - allocate(eta1_ll(nsig+1),eta2_ll(nsig+1)) - allocate(ak(nz),bk(nz),abk_fv3(nz)) + if(.not.allocated(aeta1_ll))allocate(aeta1_ll(nsig)) + if(.not.allocated(aeta2_ll))allocate(aeta2_ll(nsig)) + if(.not.allocated(eta1_ll))allocate(eta1_ll(nsig+1)) + if(.not.allocated(eta2_ll))allocate(eta2_ll(nsig+1)) + if(.not.allocated(ak))allocate(ak(nz)) + if(.not.allocated(bk))allocate(bk(nz)) + if(.not.allocated(abk_fv3))allocate(abk_fv3(nz)) do k=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,k,name,len) @@ -3635,10 +3639,6 @@ subroutine m_gsi_rfv3io_get_grid_specs(gsi_lats,gsi_lons,ierr) integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout - !coupler_res_filenam='/home/masanori/da/RDASApp_gsib/jedi-assim_test_gsib/rrfs-data_fv3jedi_2022052619/Data/gsibec/coupler.res' - !grid_spec='/home/masanori/da/RDASApp_gsib/jedi-assim_test_gsib/rrfs-data_fv3jedi_2022052619/Data/gsibec/fv3_grid_spec' - !ak_bk='/home/masanori/da/RDASApp_gsib/jedi-assim_test_gsib/rrfs-data_fv3jedi_2022052619/Data/gsibec/fv3_akbk' - coupler_res_filenam='coupler.res' grid_spec='fv3_grid_spec' ak_bk='fv3_akbk' @@ -3837,11 +3837,7 @@ subroutine m_gsi_rfv3io_get_grid_specs(gsi_lats,gsi_lons,ierr) call m_generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt,gsi_lats,gsi_lons) deallocate (grid_lon,grid_lat,grid_lont,grid_latt) - deallocate (ak,bk,abk_fv3) - deallocate(ny_layout_len,ny_layout_b,ny_layout_e) - deallocate(aeta1_ll,aeta2_ll) - deallocate(eta1_ll,eta2_ll) return end subroutine m_gsi_rfv3io_get_grid_specs diff --git a/src/gsibec/gsi/gsimod.F90 b/src/gsibec/gsi/gsimod.F90 index 08b67ea..d1777f7 100644 --- a/src/gsibec/gsi/gsimod.F90 +++ b/src/gsibec/gsi/gsimod.F90 @@ -36,6 +36,7 @@ module gsimod bkgv_flowdep,bkgv_rewgtfct,bkgv_write,fpsproj,nhscrf,adjustozvar,fut2ps,cwcoveqqcov,adjustozhscl,& bkgv_write_cv,bkgv_write_sv use berror, only: simcv !_RT intro for testing + use m_berror_stats, only: usenewgfsberror use compact_diffs, only: noq,init_compact_diffs use gridmod, only: nlat,nlon,nsig,& @@ -417,7 +418,7 @@ module gsimod namelist/bkgerr/vs,nhscrf,hzscl,hswgt,norh,ndeg,noq,bw,norsp,fstat,pert_berr,pert_berr_fct, & bkgv_flowdep,bkgv_rewgtfct,bkgv_write,fpsproj,adjustozvar,fut2ps,cwcoveqqcov,adjustozhscl,& - simcv,bkgv_write_cv,bkgv_write_sv + simcv,bkgv_write_cv,bkgv_write_sv,usenewgfsberror ! STRONGOPTS (strong dynamic constraint) ! reg_tlnmc_type - =1 for 1st version of regional strong constraint @@ -702,10 +703,6 @@ subroutine gsimain_initialize_(nfldsig,nmlfile) if (mype==0) write(6,*)'GSIMOD: tendencies and derivatives are on' endif - if (regional) then - call convert_fv3_regional - endif - ! Initialize variables, create/initialize arrays lendian_in = -1 call create_ges_tendencies(tendsflag,thisrc) diff --git a/src/gsibec/gsi/guess_grids.f90 b/src/gsibec/gsi/guess_grids.f90 index a68b6f2..52c28da 100644 --- a/src/gsibec/gsi/guess_grids.f90 +++ b/src/gsibec/gsi/guess_grids.f90 @@ -17,6 +17,7 @@ module guess_grids use mod_vtrans, only: nvmodes_keep,create_vtrans use mod_strong, only: l_tlnmc use strong_fast_global_mod, only: init_strongvars +use m_mpimod, only: nxpe,nype implicit none private ! @@ -246,6 +247,7 @@ subroutine bkgcov_init_(need) logical, save :: init_pass = .true. call other_set_(need=need) ! a little out of place, but ... call compute_derived(mype,init_pass) ! this belongs in a state set + if (l_tlnmc .and. nvmodes_keep>0) then call create_vtrans(mype,ntguessig) ! if(regional) then @@ -654,6 +656,8 @@ subroutine load_prsges_ ihaveprs(jj)=.true. end do + if(regional) then + if (fv3_regional) then do jj=1,nfldsig do k=1,nsig @@ -669,6 +673,8 @@ subroutine load_prsges_ end do end if ! end if fv3 regional + else + ! load mid-layer pressure by using phillips vertical interpolation if (idsl5/=2) then do jj=1,nfldsig @@ -703,6 +709,8 @@ subroutine load_prsges_ end do endif + end if + ! For regional applications only, load variables containing mean ! surface pressure and pressure profile at the layer midpoints if (regional) then @@ -1101,35 +1109,182 @@ subroutine guess_basics0_ end subroutine guess_basics0_ !-------------------------------------------------------- subroutine guess_basics2_(vname,islot,var) + use gridmod, only: regional character(len=*),intent(in) :: vname integer(i_kind), intent(in) :: islot real(r_kind),dimension(:,:) :: var character(len=*), parameter :: myname_ = myname//'*guess_basics2_' real(r_kind),dimension(:,:),pointer::ptr - integer jj,ier + integer jj,ier,i,j jj=islot call gsi_bundlegetpointer(gsi_metguess_bundle(jj),trim(vname),ptr,ier) if (ier/=0) then call die(myname_,'pointer to '//trim(vname)//" not found",ier) endif + + if(regional) then + if(mype == 0) then + !$omp parallel do default(shared) private(i) + do i = 1, size(var,2) + call put_data_1d(var(:,i)) + enddo + !$omp end parallel do + !$omp parallel do default(shared) private(j) + do j = 1, size(var,1) + call put_data_1d(var(j,:)) + enddo + !$omp end parallel do + else if(mype == nxpe-1) then + !$omp parallel do default(shared) private(i) + do i = 1, size(var,2) + call put_data_1d(var(:,i)) + enddo + !$omp end parallel do + !$omp parallel do default(shared) private(j) + do j = 1, size(var,1) + call put_data_1d_rev(var(j,:)) + enddo + !$omp end parallel do + else if(mype == nxpe*(nype-1)) then + !$omp parallel do default(shared) private(i) + do i = 1, size(var,2) + call put_data_1d_rev(var(:,i)) + enddo + !$omp end parallel do + !$omp parallel do default(shared) private(j) + do j = 1, size(var,1) + call put_data_1d(var(j,:)) + enddo + !$omp end parallel do + else if(mype == nxpe*nype-1) then + !$omp parallel do default(shared) private(i) + do i = 1, size(var,2) + call put_data_1d_rev(var(:,i)) + enddo + !$omp end parallel do + !$omp parallel do default(shared) private(j) + do j = 1, size(var,1) + call put_data_1d_rev(var(j,:)) + enddo + !$omp end parallel do + else if(mype>0 .and. mypenxpe*(nype-1) .and. mype0 .and. mypenxpe-1 .and. mype0 .and. mypenxpe*(nype-1) .and. mype0 .and. mypenxpe-1 .and. mype rmiss_th) then + first_valid = k + exit + end if + end do + + if (first_valid <= 0) return + + if (first_valid > 1) then + x(1:first_valid-1) = x(first_valid) + end if + + end subroutine put_data_1d +!-------------------------------------------------------- + subroutine put_data_1d_rev(x) + implicit none + real(r_kind), intent(inout) :: x(:) + real(r_kind), parameter :: rmiss = -3.334767057904812e38 + real(r_kind), parameter :: rmiss_th = -1.0e30 + integer :: n, k, first_valid + + n = size(x) + + first_valid = 0 + do k = n, 1, -1 + if (x(k) > rmiss_th) then + first_valid = k + exit + end if + end do + + if (first_valid <= 0) return + + if (first_valid > 1) then + x(first_valid+1:n) = x(first_valid) + end if + + end subroutine put_data_1d_rev +!-------------------------------------------------------- + subroutine put_data_2d(x) + implicit none + real(r_kind), intent(inout) :: x(:,:) + real(r_kind), parameter :: rmiss = -3.334767057904812e38 + real(r_kind), parameter :: rmiss_th = -1.0e30 + integer :: n, k, first_valid, i + + n = size(x,1) + + first_valid = 0 + do k = 1, n + if (x(k,1) > rmiss_th) then + first_valid = k + exit + end if + end do + + if (first_valid <= 0) return + + if (first_valid > 1) then + do i = 1, first_valid-1 + x(i,:) = x(first_valid,:) + enddo + end if + + end subroutine put_data_2d +!-------------------------------------------------------- + subroutine put_data_2d_rev(x) + implicit none + real(r_kind), intent(inout) :: x(:,:) + real(r_kind), parameter :: rmiss = -3.334767057904812e38 + real(r_kind), parameter :: rmiss_th = -1.0e30 + integer :: n, k, first_valid, i + + n = size(x,1) + + first_valid = 0 + do k = n, 1, -1 + if (x(k,1) > rmiss_th) then + first_valid = k + exit + end if + end do + + if (first_valid <= 0) return + + if (first_valid > 1) then + do i = first_valid+1, n + x(i,:) = x(first_valid,:) + enddo + end if + + end subroutine put_data_2d_rev !-------------------------------------------------------- end module guess_grids diff --git a/src/gsibec/gsi/mod_fv3_lola.f90 b/src/gsibec/gsi/mod_fv3_lola.f90 index 14913a7..e0f2435 100644 --- a/src/gsibec/gsi/mod_fv3_lola.f90 +++ b/src/gsibec/gsi/mod_fv3_lola.f90 @@ -209,10 +209,8 @@ subroutine generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt) !--------------------------obtain analysis grid dimensions nxa,nya nxa=1+nint((nx-one)/grid_ratio_fv3_regional) nya=1+nint((ny-one)/grid_ratio_fv3_regional) - !nlat=nya - !nlon=nxa - nlat=ny - nlon=nx + nlat=nya + nlon=nxa if(mype==0) print *,'nlat,nlon=nya,nxa= ',nlat,nlon !--------------------------obtain analysis grid spacing @@ -261,6 +259,8 @@ subroutine generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt) if (allocated(region_dx )) deallocate(region_dx ) if (allocated(region_dy )) deallocate(region_dy ) + if (allocated(coeffx )) deallocate(coeffx ) + if (allocated(coeffy )) deallocate(coeffy ) allocate(region_dx(nlat,nlon),region_dy(nlat,nlon)) allocate(region_dxi(nlat,nlon),region_dyi(nlat,nlon)) allocate(coeffx(nlat,nlon),coeffy(nlat,nlon)) @@ -712,10 +712,8 @@ subroutine m_generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt,gsi_l !--------------------------obtain analysis grid dimensions nxa,nya nxa=1+nint((nx-one)/grid_ratio_fv3_regional) nya=1+nint((ny-one)/grid_ratio_fv3_regional) - !nlat=nya - !nlon=nxa - nlat=ny - nlon=nx + nlat=nya + nlon=nxa if(mype==0) print *,'nlat,nlon=nya,nxa= ',nlat,nlon !--------------------------obtain analysis grid spacing @@ -818,15 +816,9 @@ subroutine m_generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt,gsi_l enddo enddo - !call init_general_transform(glat_an,glon_an) - - !deallocate(glat_an,glon_an) - deallocate( xc,yc,zc,gclat,gclon,gcrlat,gcrlon) deallocate(rlat_in,rlon_in) - deallocate(region_dxi,region_dyi) - deallocate(coeffx,coeffy) end subroutine m_generate_anl_grid