From 42dd1fe3b7807c2e015230a1b1884b0b052d9836 Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Fri, 3 Nov 2023 14:52:09 -0500 Subject: [PATCH] Squashed changes to add the dual-resolution EnVar capability for HAFS ensemble. From Xu Lu, OU (luxu@ou.edu) and POC: xuguang.wang@ou.edu --- src/gsi/constants.f90 | 1 + src/gsi/cplr_get_fv3_regional_ensperts.f90 | 17 +- src/gsi/gridmod.F90 | 6 +- src/gsi/gsi_rfv3io_mod.f90 | 359 +++++++++++-- src/gsi/hybrid_ensemble_isotropic.F90 | 15 +- src/gsi/mod_fv3_lola.f90 | 556 +++++++++++++++++++++ 6 files changed, 894 insertions(+), 60 deletions(-) diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index 291a18ec97..c1c5f3f879 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -31,6 +31,7 @@ module constants ! 2016-02-15 Johnson, Y. Wang, X. Wang - define additional constant values for ! radar DA, POC: xuguang.wang@ou.edu ! 2019-09-25 X.Su - put stndrd_atmos_ps constant values +! 2022-03-01 X.Lu & X.Wang - increased max_varname_length for HAFS dual ens. POC: xuguang.wang@ou.edu ! ! Subroutines Included: ! sub init_constants_derived - compute derived constants diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index 2382ff1286..81fb684a73 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -38,6 +38,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! ! 2021-08-10 lei - modify for fv3-lam ensemble spread output ! 2021-11-01 lei - modify for fv3-lam parallel IO + ! 2022-03-01 X.Lu & X.Wang - modify for hafs dual ens. POC: xuguang.wang@ou.edu ! input argument list: ! ! output argument list: @@ -848,7 +849,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g use gridmod, only: eta1_ll,eta2_ll use constants, only: zero,one,fv,zero_single,one_tenth,h300 use hybrid_ensemble_parameters, only: grd_ens,q_hyb_ens - use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt + use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt,dual_res use mpimod, only: mpi_comm_world,mpi_rtype use gsi_rfv3io_mod,only: type_fv3regfilenameg @@ -956,24 +957,24 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g if(fv3sar_ensemble_opt == 0 ) then - call gsi_fv3ncdf_readuv(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput) + call gsi_fv3ncdf_readuv(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput,dual_res) else - call gsi_fv3ncdf_readuv_v1(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput) + call gsi_fv3ncdf_readuv_v1(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput,dual_res) endif if(fv3sar_ensemble_opt == 0) then call gsi_fv3ncdf_read(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& - fv3_filenameginput%dynvars,fv3_filenameginput) + fv3_filenameginput%dynvars,fv3_filenameginput,dual_res) call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& - fv3_filenameginput%tracers,fv3_filenameginput) + fv3_filenameginput%tracers,fv3_filenameginput,dual_res) if( if_model_dbz .or. if_model_fed ) then call gsi_fv3ncdf_read(grd_fv3lam_ens_phyvar_io_nouv,gsibundle_fv3lam_ens_phyvar_nouv,& - fv3_filenameginput%phyvars,fv3_filenameginput) + fv3_filenameginput%phyvars,fv3_filenameginput,dual_res) end if else call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& - fv3_filenameginput%dynvars,fv3_filenameginput) + fv3_filenameginput%dynvars,fv3_filenameginput,dual_res) call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& - fv3_filenameginput%tracers,fv3_filenameginput) + fv3_filenameginput%tracers,fv3_filenameginput,dual_res) endif ier=0 call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus diff --git a/src/gsi/gridmod.F90 b/src/gsi/gridmod.F90 index 559a3f576d..2367899ea5 100644 --- a/src/gsi/gridmod.F90 +++ b/src/gsi/gridmod.F90 @@ -93,6 +93,7 @@ module gridmod ! 2019-09-23 martin - add use_gfs_ncio to read global first guess from netCDF file ! 2020-12-18 Hu - add grid_type_fv3_regional ! 2021-12-30 Hu - add fv3_io_layout_y +! 2022-03-01 X.Lu & X.Wang - add corresponding variables for dual ens for HAFS. POC: xuguang.wang@ou.edu ! ! ! @@ -146,6 +147,7 @@ module gridmod public :: regional_fhr,region_dyi,coeffx,region_dxi,coeffy,nsig_hlf,regional_fmin public :: nsig2,wgtlats,corlats,rbs2,ncepgfs_headv,regional_time,wgtfactlats public :: nlat_regional,nlon_regional,update_regsfc,half_grid,gencode + public :: nlat_regionalens,nlon_regionalens public :: diagnostic_reg,nmmb_reference_grid,filled_grid public :: grid_ratio_nmmb,isd_g,isc_g,dx_gfs,lpl_gfs,nsig5,nmmb_verttype public :: grid_ratio_fv3_regional,fv3_io_layout_y,fv3_regional,fv3_cmaq_regional,grid_type_fv3_regional @@ -329,7 +331,7 @@ module gridmod real(r_kind) rlon_min_dd,rlon_max_dd,rlat_min_dd,rlat_max_dd real(r_kind) dt_ll,pdtop_ll,pt_ll - integer(i_kind) nlon_regional,nlat_regional + integer(i_kind) nlon_regional,nlat_regional,nlon_regionalens,nlat_regionalens real(r_kind) regional_fhr,regional_fmin integer(i_kind) regional_time(6) integer(i_kind) jcap_gfs,nlat_gfs,nlon_gfs @@ -485,6 +487,8 @@ subroutine init_grid update_regsfc = .false. nlon_regional = 0 nlat_regional = 0 + nlon_regionalens = 0 + nlat_regionalens = 0 msig = nsig do k=1,size(nlayers) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index e62cc06f2b..43fb40e8ad 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -18,6 +18,7 @@ module gsi_rfv3io_mod ! This function is needed when fv3 model sets ! io_layout(2)>1 ! 2022-02-15 Lu @ Wang - add time label it for FGAT. POC: xuguang.wang@ou.edu +! 2022-03-01 X.Lu @ X.Wang - add gsi_rfv3io_get_ens_grid_specs for dual ens HAFS. POC: xuguang.wang@ou.edu ! 2022-03-15 Hu - add code to read/write 2m T and Q for they will be ! used as background for surface observation operator ! 2022-04-15 Wang - add IO for regional FV3-CMAQ (RRFS-CMAQ) model @@ -27,6 +28,7 @@ module gsi_rfv3io_mod ! ! subroutines included: ! sub gsi_rfv3io_get_grid_specs +! sub gsi_rfv3io_get_ens_grid_specs ! sub read_fv3_files ! sub read_fv3_netcdf_guess ! sub gsi_fv3ncdf2d_read @@ -49,7 +51,7 @@ module gsi_rfv3io_mod !$$$ end documentation block use kinds, only: r_kind,i_kind - use gridmod, only: nlon_regional,nlat_regional + use gridmod, only: nlon_regional,nlat_regional,nlon_regionalens,nlat_regionalens use constants, only:max_varname_length,max_filename_length use gsi_bundlemod, only : gsi_bundle use general_sub2grid_mod, only: sub2grid_info @@ -82,7 +84,9 @@ module gsi_rfv3io_mod type(type_fv3regfilenameg),allocatable:: bg_fv3regfilenameg(:) integer(i_kind) nx,ny,nz + integer(i_kind) nxens,nyens integer(i_kind),dimension(:),allocatable :: ny_layout_len,ny_layout_b,ny_layout_e + integer(i_kind),dimension(:),allocatable :: ny_layout_lenens,ny_layout_bens,ny_layout_eens real(r_kind),allocatable:: grid_lon(:,:),grid_lont(:,:),grid_lat(:,:),grid_latt(:,:) real(r_kind),allocatable:: ak(:),bk(:) integer(i_kind),allocatable:: ijns2d(:),displss2d(:),ijns(:),displss(:) @@ -122,6 +126,7 @@ module gsi_rfv3io_mod private ! set subroutines to public public :: gsi_rfv3io_get_grid_specs + public :: gsi_rfv3io_get_ens_grid_specs public :: gsi_fv3ncdf_read public :: gsi_fv3ncdf_read_v1 public :: gsi_fv3ncdf_readuv @@ -500,6 +505,179 @@ subroutine gsi_rfv3io_get_grid_specs(ierr) return end subroutine gsi_rfv3io_get_grid_specs +subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr) +!$$$ subprogram documentation block +! . . . . +! subprogram: gsi_rfv3io_get_ens_grid_specs +! modified from gsi_rfv3io_get_grid_specs +! pgrmmr: parrish org: np22 date: 2017-04-03 +! +! abstract: obtain grid dimensions nx,ny and grid definitions +! grid_x,grid_xt,grid_y,grid_yt,grid_lon,grid_lont,grid_lat,grid_latt +! nz,ak(nz),bk(nz) +! +! program history log: +! 2017-04-03 parrish - initial documentation +! 2017-10-10 wu - setup A grid and interpolation coeff with generate_anl_grid +! 2018-02-16 wu - read in time info from file coupler.res +! read in lat, lon at the center and corner of the grid cell +! from file fv3_grid_spec, and vertical grid infor from file +! fv3_akbk +! setup A grid and interpolation/rotation coeff +! input argument list: +! grid_spec +! ak_bk +! lendian_out +! +! output argument list: +! ierr +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr + use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_inquire_variable + use mpimod, only: mype + use mod_fv3_lola, only: definecoef_regular_grids + use gridmod, only:nsig,regional_time,regional_fhr,regional_fmin,aeta1_ll,aeta2_ll + use gridmod, only:nlon_regionalens,nlat_regionalens + use gridmod, only:grid_type_fv3_regional + use kinds, only: i_kind,r_kind + use constants, only: half,zero + use mpimod, only: mpi_comm_world,mpi_itype,mpi_rtype + implicit none + character(:),allocatable,intent(in) :: grid_spec + + integer(i_kind) gfile_grid_spec + integer(i_kind),intent( out) :: ierr + integer(i_kind) i,k,ndimensions,iret,nvariables,nattributes,unlimiteddimid + integer(i_kind) gfile_loc,len + character(len=128) :: name + integer(i_kind) :: nio,nylen + integer(i_kind),allocatable :: gfile_loc_layout(:) + character(len=180) :: filename_layout + integer(i_kind) imiddle,jmiddle + + + iret=nf90_open(trim(grid_spec),nf90_nowrite,gfile_grid_spec) + if(iret/=nf90_noerr) then + write(6,*)' problem opening1 ',trim(grid_spec),', Status = ',iret + ierr=1 + return + endif + iret=nf90_inquire(gfile_grid_spec,ndimensions,nvariables,nattributes,unlimiteddimid) + gfile_loc=gfile_grid_spec + do k=1,ndimensions + iret=nf90_inquire_dimension(gfile_loc,k,name,len) + if(trim(name)=='grid_xt') nxens=len + if(trim(name)=='grid_yt') nyens=len + enddo + allocate(grid_lat(nxens+1,nyens+1)) + allocate(grid_lon(nxens+1,nyens+1)) + allocate(grid_latt(nxens,nyens)) + allocate(grid_lont(nxens,nyens)) + do k=ndimensions+1,nvariables + iret=nf90_inquire_variable(gfile_loc,k,name,len) + if(trim(name)=='grid_lat') then + iret=nf90_get_var(gfile_loc,k,grid_lat) + endif + if(trim(name)=='grid_lon') then + iret=nf90_get_var(gfile_loc,k,grid_lon) + endif + if(trim(name)=='grid_latt') then + iret=nf90_get_var(gfile_loc,k,grid_latt) + endif + if(trim(name)=='grid_lont') then + iret=nf90_get_var(gfile_loc,k,grid_lont) + endif + enddo + iret=nf90_close(gfile_loc) + + nlon_regionalens=nxens + nlat_regionalens=nyens + allocate(ny_layout_lenens(0:fv3_io_layout_y-1)) + allocate(ny_layout_bens(0:fv3_io_layout_y-1)) + allocate(ny_layout_eens(0:fv3_io_layout_y-1)) + ny_layout_lenens=nyens + ny_layout_bens=0 + ny_layout_eens=0 + if(fv3_io_layout_y > 1) then + allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) + do nio=0,fv3_io_layout_y-1 + write(filename_layout,'(a,a,I4.4)') trim(grid_spec),'.',nio + iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio)) + if(iret/=nf90_noerr) then + write(6,*)' problem opening ',trim(filename_layout),', Status =',iret + ierr=1 + return + endif + iret=nf90_inquire(gfile_loc_layout(nio),ndimensions,nvariables,nattributes,unlimiteddimid) + do k=1,ndimensions + iret=nf90_inquire_dimension(gfile_loc_layout(nio),k,name,len) + if(trim(name)=='grid_yt') ny_layout_lenens(nio)=len + enddo + iret=nf90_close(gfile_loc_layout(nio)) + enddo + deallocate(gfile_loc_layout) +! figure out begin and end of each subdomain restart file + nylen=0 + do nio=0,fv3_io_layout_y-1 + ny_layout_bens(nio)=nylen + 1 + nylen=nylen+ny_layout_lenens(nio) + ny_layout_eens(nio)=nylen + enddo + endif + if(mype==0)write(6,*),'nxens,nyens=',nxens,nyens + if(mype==0)write(6,*),'ny_layout_lenens=',ny_layout_lenens + if(mype==0)write(6,*),'ny_layout_bens=',ny_layout_bens + if(mype==0)write(6,*),'ny_layout_eens=',ny_layout_eens + +! +! need to decide the grid orientation of the FV regional model +! +! grid_type_fv3_regional = 0 : decide grid orientation based on +! grid_lat/grid_lon +! 1 : input is E-W N-S grid +! 2 : input is W-E S-N grid +! + if(grid_type_fv3_regional == 0) then + imiddle=nxens/2 + jmiddle=nyens/2 + if( (grid_latt(imiddle,1) < grid_latt(imiddle,nyens)) .and. & + (grid_lont(1,jmiddle) < grid_lont(nxens,jmiddle)) ) then + grid_type_fv3_regional = 2 + else + grid_type_fv3_regional = 1 + endif + endif +! check the grid type + if( grid_type_fv3_regional == 1 ) then + if(mype==0) write(6,*) 'FV3 regional input grid is E-W N-S grid' + grid_reverse_flag=.true. ! grid is revered comparing to usual map view + else if(grid_type_fv3_regional == 2) then + if(mype==0) write(6,*) 'FV3 regional input grid is W-E S-N grid' + grid_reverse_flag=.false. ! grid orientated just like we see on map view + else + write(6,*) 'Error: FV3 regional input grid is unknown grid' + call stop2(678) + endif +! + if(grid_type_fv3_regional == 2) then + call reverse_grid_r(grid_lont,nxens,nyens,1) + call reverse_grid_r(grid_latt,nxens,nyens,1) + call reverse_grid_r(grid_lon,nxens+1,nyens+1,1) + call reverse_grid_r(grid_lat,nxens+1,nyens+1,1) + endif + + call definecoef_regular_grids(nxens,nyens,grid_lon,grid_lont,grid_lat,grid_latt) + deallocate (grid_lon,grid_lat,grid_lont,grid_latt) + return +end subroutine gsi_rfv3io_get_ens_grid_specs + + subroutine read_fv3_files(mype) !$$$ subprogram documentation block ! . . . . @@ -1110,7 +1288,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif end do if (ndynvario2d > 0) then - allocate(fv3lam_io_dynmetvars2d_nouv(ndynvario2d)) + if (.not. allocated(fv3lam_io_dynmetvars2d_nouv)) then + allocate(fv3lam_io_dynmetvars2d_nouv(ndynvario2d)) + end if endif if (ntracerio2d > 0) then allocate(fv3lam_io_tracermetvars2d_nouv(ntracerio2d)) @@ -1421,40 +1601,40 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) end if if( fv3sar_bg_opt == 0) then - call gsi_fv3ncdf_readuv(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it)) + call gsi_fv3ncdf_readuv(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it),.false.) else - call gsi_fv3ncdf_readuv_v1(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it)) + call gsi_fv3ncdf_readuv_v1(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it),.false.) endif if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_read(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv & - & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) + & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it),.false.) call gsi_fv3ncdf_read(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) if( nphyvario3d > 0 )then call gsi_fv3ncdf_read(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv & - & ,fv3filenamegin(it)%phyvars,fv3filenamegin(it)) + & ,fv3filenamegin(it)%phyvars,fv3filenamegin(it),.false.) end if if (laeroana_fv3cmaq) then call gsi_fv3ncdf_read(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) endif if (laeroana_fv3smoke) then call gsi_fv3ncdf_read(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) endif else call gsi_fv3ncdf_read_v1(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv & - & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) + & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it),.false.) call gsi_fv3ncdf_read_v1(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) if (laeroana_fv3cmaq) then call gsi_fv3ncdf_read_v1(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) endif if (laeroana_fv3smoke) then call gsi_fv3ncdf_read_v1(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv & - & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it),.false.) endif endif @@ -2220,7 +2400,7 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io) return end subroutine gsi_fv3ncdf2d_read_v1 -subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) +subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ensgrid) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_read @@ -2253,7 +2433,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid - use mod_fv3_lola, only: fv3_h_to_ll + use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_grid2sub @@ -2262,6 +2442,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) type(gsi_bundle),intent(inout) :: cstate_nouv character(*),intent(in):: filenamein type (type_fv3regfilenameg),intent(in) ::fv3filenamegin + logical, intent(in ) :: ensgrid real(r_kind),allocatable,dimension(:,:):: uu2d real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork character(len=max_varname_length) :: varname,vgsiname @@ -2290,8 +2471,13 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) mm1=mype+1 nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat - nxcase=nx - nycase=ny + if (ensgrid) then + nxcase=nxens + nycase=nyens + else + nxcase=nx + nycase=ny + end if kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate(uu2d(nxcase,nycase)) @@ -2376,11 +2562,20 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - countloc=(/nxcase,ny_layout_len(nio),1/) - allocate(uu2d_layout(nxcase,ny_layout_len(nio))) + if (ensgrid) then + countloc=(/nxcase,ny_layout_lenens(nio)+1,1/) + allocate(uu2d_layout(nxcase,ny_layout_lenens(nio)+1)) + else + countloc=(/nxcase,ny_layout_len(nio),1/) + allocate(uu2d_layout(nxcase,ny_layout_len(nio))) + end if iret=nf90_inq_varid(gfile_loc_layout(nio),trim(adjustl(varname)),var_id) iret=nf90_get_var(gfile_loc_layout(nio),var_id,uu2d_layout,start=startloc,count=countloc) - uu2d(:,ny_layout_b(nio):ny_layout_e(nio))=uu2d_layout + if (ensgrid) then + uu2d(:,ny_layout_bens(nio):ny_layout_eens(nio))=uu2d_layout + else + uu2d(:,ny_layout_b(nio):ny_layout_e(nio))=uu2d_layout + end if deallocate(uu2d_layout) enddo else @@ -2403,7 +2598,11 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) end if endif - call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + if (ensgrid) then + call fv3_h_to_ll_ens(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + else + call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + endif enddo ! ilevtot if(fv3_io_layout_y > 1) then @@ -2423,7 +2622,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) return end subroutine gsi_fv3ncdf_read -subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) +subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ensgrid) !$$$ subprogram documentation block ! . . . . @@ -2458,13 +2657,14 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid - use mod_fv3_lola, only: fv3_h_to_ll + use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none type(sub2grid_info), intent(in):: grd_ionouv character(*),intent(in):: filenamein + logical, intent(in ) :: ensgrid type (type_fv3regfilenameg) :: fv3filenamegin type(gsi_bundle),intent(inout) :: cstate_nouv real(r_kind),allocatable,dimension(:,:):: uu2d @@ -2484,8 +2684,13 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat - nxcase=nx - nycase=ny + if (ensgrid) then + nxcase=nxens + nycase=nyens + else + nxcase=nx + nycase=ny + end if kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate(uu2d(nxcase,nycase)) @@ -2519,7 +2724,11 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) - call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + if (ensgrid) then + call fv3_h_to_ll_ens(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + else + call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + end if enddo ! i call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) @@ -2531,7 +2740,7 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) return end subroutine gsi_fv3ncdf_read_v1 -subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) +subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_readuv @@ -2558,7 +2767,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid - use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth + use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth,fv3_h_to_ll_ens,fv3uv2earthens use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none @@ -2566,6 +2775,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) real(r_kind),dimension(grd_uv%lat2,grd_uv%lon2,grd_uv%nsig),intent(inout)::ges_u real(r_kind),dimension(grd_uv%lat2,grd_uv%lon2,grd_uv%nsig),intent(inout)::ges_v type (type_fv3regfilenameg),intent (in) :: fv3filenamegin + logical, intent(in ) :: ensgrid real(r_kind),dimension(2,grd_uv%nlat,grd_uv%nlon,grd_uv%kbegin_loc:grd_uv%kend_alloc):: hwork character(:), allocatable:: filenamein real(r_kind),allocatable,dimension(:,:):: u2d,v2d @@ -2596,8 +2806,13 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat - nxcase=nx - nycase=ny + if (ensgrid) then + nxcase=nxens + nycase=nyens + else + nxcase=nx + nycase=ny + end if kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate(u2d(nxcase,nycase+1)) @@ -2667,19 +2882,35 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) v_startloc=(/1,1,inative/) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) - allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) + if (ensgrid) then + u_countloc=(/nxcase,ny_layout_lenens(nio)+1,1/) + allocate(u2d_layout(nxcase,ny_layout_lenens(nio)+1)) + else + u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) + allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) + end if call check( nf90_inq_varid(gfile_loc_layout(nio),'u',u_grd_VarId) ) iret=nf90_get_var(gfile_loc_layout(nio),u_grd_VarId,u2d_layout,start=u_startloc,count=u_countloc) - u2d(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) - if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) - deallocate(u2d_layout) - - v_countloc=(/nxcase+1,ny_layout_len(nio),1/) - allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) + if (ensgrid) then + u2d(:,ny_layout_bens(nio):ny_layout_eens(nio))=u2d_layout(:,1:ny_layout_lenens(nio)) + if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_eens(nio)+1)=u2d_layout(:,ny_layout_lenens(nio)+1) + deallocate(u2d_layout) + v_countloc=(/nxcase+1,ny_layout_lenens(nio),1/) + allocate(v2d_layout(nxcase+1,ny_layout_lenens(nio))) + else + u2d(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) + if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) + deallocate(u2d_layout) + v_countloc=(/nxcase+1,ny_layout_len(nio),1/) + allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) + end if call check( nf90_inq_varid(gfile_loc_layout(nio),'v',v_grd_VarId) ) iret=nf90_get_var(gfile_loc_layout(nio),v_grd_VarId,v2d_layout,start=v_startloc,count=v_countloc) - v2d(:,ny_layout_b(nio):ny_layout_e(nio))=v2d_layout + if (ensgrid) then + v2d(:,ny_layout_bens(nio):ny_layout_eens(nio))=v2d_layout + else + v2d(:,ny_layout_b(nio):ny_layout_e(nio))=v2d_layout + end if deallocate(v2d_layout) enddo else @@ -2693,7 +2924,11 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) call reverse_grid_r_uv (u2d,nxcase,nycase+1,1) call reverse_grid_r_uv (v2d,nxcase+1,nycase,1) endif - call fv3uv2earth(u2d(:,:),v2d(:,:),nxcase,nycase,uc2d,vc2d) + if (ensgrid) then + call fv3uv2earthens(u2d(:,:),v2d(:,:),nxcase,nycase,uc2d,vc2d) + else + call fv3uv2earth(u2d(:,:),v2d(:,:),nxcase,nycase,uc2d,vc2d) + end if ! NOTE on transfor to earth u/v: ! The u and v before transferring need to be in E-W/N-S grid, which is @@ -2711,8 +2946,13 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) ! and the last input parameter for fv3_h_to_ll is alway true: ! ! - call fv3_h_to_ll(uc2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) - call fv3_h_to_ll(vc2d,hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + if (ensgrid) then + call fv3_h_to_ll_ens(uc2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + call fv3_h_to_ll_ens(vc2d,hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + else + call fv3_h_to_ll(uc2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + call fv3_h_to_ll(vc2d,hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + end if enddo ! i if(fv3_io_layout_y > 1) then @@ -2734,7 +2974,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) deallocate(worksub) end subroutine gsi_fv3ncdf_readuv -subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) +subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) !$$$ subprogram documentation block ! subprogram: gsi_fv3ncdf_readuv_v1 ! prgmmr: wu w org: np22 date: 2017-11-22 @@ -2762,7 +3002,7 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid - use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth + use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none @@ -2770,6 +3010,7 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) real(r_kind) ,intent(out ) :: ges_u(grd_uv%lat2,grd_uv%lon2,grd_uv%nsig) real(r_kind) ,intent(out ) :: ges_v(grd_uv%lat2,grd_uv%lon2,grd_uv%nsig) type (type_fv3regfilenameg),intent (in) :: fv3filenamegin + logical, intent(in ) :: ensgrid real(r_kind),dimension(2,grd_uv%nlat,grd_uv%nlon,grd_uv%kbegin_loc:grd_uv%kend_alloc):: hwork character(len=:),allocatable :: filenamein real(r_kind),allocatable,dimension(:,:):: us2d,vw2d @@ -2792,8 +3033,13 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat - nxcase=nx - nycase=ny + if (ensgrid) then + nxcase=nxens + nycase=nyens + else + nxcase=nx + nycase=ny + end if kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate (us2d(nxcase,nycase+1),vw2d(nxcase+1,nycase)) @@ -2818,8 +3064,13 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) nz=grd_uv%nsig nzp1=nz+1 inative=nzp1-ilev - us_countloc= (/nlon_regional,nlat_regional+1,1/) - vw_countloc= (/nlon_regional+1,nlat_regional,1/) + if (ensgrid) then + us_countloc= (/nlon_regionalens,nlat_regionalens+1,1/) + vw_countloc= (/nlon_regionalens+1,nlat_regionalens,1/) + else + us_countloc= (/nlon_regional,nlat_regional+1,1/) + vw_countloc= (/nlon_regional+1,nlat_regional,1/) + end if us_startloc=(/1,1,inative+1/) vw_startloc=(/1,1,inative+1/) @@ -2834,11 +3085,19 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) uorv2d(:,j)=half*(us2d(:,j)+us2d(:,j+1)) enddo - call fv3_h_to_ll(uorv2d(:,:),hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + if (ensgrid) then + call fv3_h_to_ll_ens(uorv2d(:,:),hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + else + call fv3_h_to_ll(uorv2d(:,:),hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + end if do j=1,nx uorv2d(j,:)=half*(vw2d(j,:)+vw2d(j+1,:)) enddo - call fv3_h_to_ll(uorv2d(:,:),hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + if (ensgrid) then + call fv3_h_to_ll_ens(uorv2d(:,:),hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + else + call fv3_h_to_ll(uorv2d(:,:),hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + end if enddo ! iilevtoto call general_grid2sub(grd_uv,hwork,worksub) diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index 4bad129a72..63af28cc5e 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -1393,6 +1393,10 @@ subroutine load_ensemble endif ! regional_ensemble_option = 5: ensembles are fv3 regional. + if (l_both_fv3sar_gfs_ens) then ! first read in gfs ensembles for regional + call get_gefs_for_regional +!clthink do we need mpi_bar here? + endif call fv3_regional_enspert%get_fv3_regional_ensperts(en_perts,nelen,ps_bar) @@ -4001,7 +4005,8 @@ subroutine hybens_grid_setup ! 2010-02-20 parrish, adapt for dual resolution ! 2011-01-30 parrish, fix so regional application depends only on parameters regional ! and dual_res. Rename subroutine get_regional_gefs_grid to get_regional_dual_res_grid. -! +! +! 2022-03-01 X.Lu & X.Wang - add vars for hafs dual ens. POC: xuguang.wang@ou.edu ! input argument list: ! ! output argument list: @@ -4027,6 +4032,8 @@ subroutine hybens_grid_setup use control_vectors, only: cvars3d,nc2d,nc3d use gridmod, only: region_lat,region_lon,region_dx,region_dy use hybrid_ensemble_parameters, only:nsclgrp,spc_multwgt,spcwgt_params,global_spectral_filter_sd + use hybrid_ensemble_parameters, only:regional_ensemble_option + use gsi_rfv3io_mod,only:gsi_rfv3io_get_ens_grid_specs implicit none @@ -4035,6 +4042,8 @@ subroutine hybens_grid_setup logical,allocatable::vector(:) real(r_kind) eps,r_e real(r_kind) rlon_a(nlon),rlat_a(nlat),rlon_e(nlon),rlat_e(nlat) + character(:),allocatable:: fv3_spec_grid_filename + integer :: ierr nord_e2a=4 ! soon, move this to hybrid_ensemble_parameters @@ -4121,6 +4130,10 @@ subroutine hybens_grid_setup else if(dual_res) then call get_region_dx_dy_ens(region_dx_ens,region_dy_ens) + if(regional_ensemble_option) then + fv3_spec_grid_filename="fv3_ens_grid_spec" + call gsi_rfv3io_get_ens_grid_specs(fv3_spec_grid_filename,ierr) + endif else region_dx_ens=region_dx region_dy_ens=region_dy diff --git a/src/gsi/mod_fv3_lola.f90 b/src/gsi/mod_fv3_lola.f90 index ebe0816c4a..4ec3c0cb93 100644 --- a/src/gsi/mod_fv3_lola.f90 +++ b/src/gsi/mod_fv3_lola.f90 @@ -18,12 +18,17 @@ module mod_fv3_lola ! fv3_ll_to_h ! 2019-11-01 wu - add checks in generate_anl_grid to present the mean ! longitude correctly to fix problem near lon=0 +! 2022-03-01 X.Lu & X.Wang - add functions for HAFS dual ens capability. POC: +! xuguang.wang@ou.edu ! ! subroutines included: ! sub generate_anl_grid +! sub definecoef_regular_grids ! sub earthuv2fv3 ! sub fv3uv2earth +! sub fv3uv2earthens ! sub fv3_h_to_ll +! sub fv3_h_to_ll_ens ! sub fv3_ll_to_h ! sub rotate2deg ! sub unrotate2deg @@ -65,6 +70,9 @@ module mod_fv3_lola public :: generate_anl_grid,fv3_h_to_ll,fv3_ll_to_h,fv3uv2earth,earthuv2fv3 public :: fv3dx,fv3dx1,fv3dy,fv3dy1,fv3ix,fv3ixp,fv3jy,fv3jyp,a3dx,a3dx1,a3dy,a3dy1,a3ix,a3ixp,a3jy,a3jyp public :: nxa,nya,cangu,sangu,cangv,sangv,nx,ny,bilinear + public :: definecoef_regular_grids,fv3_h_to_ll_ens,fv3uv2earthens + public :: fv3dxens,fv3dx1ens,fv3dyens,fv3dy1ens,fv3ixens,fv3ixpens,fv3jyens,fv3jypens,a3dxens,a3dx1ens,a3dyens,a3dy1ens,a3ixens,a3ixpens,a3jyens,a3jypens + public :: nxe,nye,canguens,sanguens,cangvens,sangvens logical bilinear integer(i_kind) nxa,nya,nx,ny @@ -73,6 +81,12 @@ module mod_fv3_lola real(r_kind) ,allocatable,dimension(:,:):: a3dx,a3dx1,a3dy,a3dy1 real(r_kind) ,allocatable,dimension(:,:):: cangu,sangu,cangv,sangv integer(i_kind),allocatable,dimension(:,:):: a3ix,a3ixp,a3jy,a3jyp + integer(i_kind) nxe,nye + real(r_kind) ,allocatable,dimension(:,:):: fv3dxens,fv3dx1ens,fv3dyens,fv3dy1ens + integer(i_kind),allocatable,dimension(:,:):: fv3ixens,fv3ixpens,fv3jyens,fv3jypens + real(r_kind) ,allocatable,dimension(:,:):: a3dxens,a3dx1ens,a3dyens,a3dy1ens + real(r_kind) ,allocatable,dimension(:,:):: canguens,sanguens,cangvens,sangvens + integer(i_kind),allocatable,dimension(:,:):: a3ixens,a3ixpens,a3jyens,a3jypens contains @@ -574,8 +588,425 @@ subroutine generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt) enddo enddo deallocate( xc,yc,zc,gclat,gclon,gcrlat,gcrlon) + deallocate(rlat_in,rlon_in) end subroutine generate_anl_grid +subroutine definecoef_regular_grids(nxen,nyen,grid_lon,grid_lont,grid_lat,grid_latt) +!$$$ subprogram documentation block +! . . . . +! subprogram: generate_??ens_grid +!clt modified from generate_regular_grid +! prgmmr: parrish +! +! abstract: define rotated lat-lon analysis grid which is centered on fv3 tile +! and oriented to completely cover the tile. +! +! program history log: +! 2017-05-02 parrish +! 2017-10-10 wu - 1. setup analysis A-grid, +! 2. compute/setup FV3 to A grid interpolation parameters +! 3. compute/setup A to FV3 grid interpolation parameters +! 4. setup weightings for wind conversion from FV3 to earth +! 2021-02-01 Lu & Wang - modify variable intent for HAFS dual ens. POC: +! xuguang.wang@ou.edu +! +! input argument list: +! nxen, nyen - number of cells = nxen*nyen +! grid_lon ,grid_lat - longitudes and latitudes of fv3 grid cell corners +! grid_lont,grid_latt - longitudes and latitudes of fv3 grid cell centers +! +! output argument list: +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + use kinds, only: r_kind,i_kind + use constants, only: quarter,one,two,half,zero,deg2rad,rearth,rad2deg + use gridmod, only:grid_ratio_fv3_regional + use mpimod, only: mype + use hybrid_ensemble_parameters, only: nlon_ens,nlat_ens,region_lon_ens,region_lat_ens + implicit none + real(r_kind),allocatable,dimension(:)::xbh_a,xa_a,xa_b + real(r_kind),allocatable,dimension(:)::ybh_a,ya_a,ya_b,yy + real(r_kind),allocatable,dimension(:,:)::xbh_b,ybh_b + real(r_kind) dlat,dlon,dyy,dxx,dyyi,dxxi + real(r_kind) dyyh,dxxh + + real(r_kind),allocatable:: region_lat_tmp(:,:),region_lon_tmp(:,:) + integer(i_kind), intent(in ) :: nxen,nyen ! fv3 tile x- and y-dimensions + real(r_kind) , intent(inout) :: grid_lon(nxen+1,nyen+1) ! fv3 cell corner longitudes + real(r_kind) , intent(inout) :: grid_lont(nxen,nyen) ! fv3 cell center longitudes + real(r_kind) , intent(inout) :: grid_lat(nxen+1,nyen+1) ! fv3 cell corner latitudes + real(r_kind) , intent(inout) :: grid_latt(nxen,nyen) ! fv3 cell center latitudes + integer(i_kind) i,j,ir,jr,n + real(r_kind),allocatable,dimension(:,:) :: xc,yc,zc,gclat,gclon,gcrlat,gcrlon,rlon_in,rlat_in + real(r_kind),allocatable,dimension(:,:) :: glon_an,glat_an + real(r_kind) xcent,ycent,zcent,rnorm,centlat,centlon + integer(i_kind) nlonh,nlath,nxh,nyh + integer(i_kind) ib1,ib2,jb1,jb2,jj + integer (i_kind):: index0 + real(r_kind) region_lat_in(nlat_ens,nlon_ens),region_lon_in(nlat_ens,nlon_ens) + integer(i_kind) nord_e2a + real(r_kind)gxa,gya + + real(r_kind) x(nxen+1,nyen+1),y(nxen+1,nyen+1),z(nxen+1,nyen+1),xr,yr,zr,xu,yu,zu,rlat,rlon + real(r_kind) xv,yv,zv,vval + real(r_kind) cx,cy + real(r_kind) uval,ewval,nsval + + real(r_kind) d(4),ds + integer(i_kind) kk,k + real(r_kind) diff,sq180 + + nord_e2a=4 + bilinear=.false. + +! create xc,yc,zc for the cell centers. + allocate(xc(nxen,nyen)) + allocate(yc(nxen,nyen)) + allocate(zc(nxen,nyen)) + allocate(gclat(nxen,nyen)) + allocate(gclon(nxen,nyen)) + allocate(gcrlat(nxen,nyen)) + allocate(gcrlon(nxen,nyen)) + do j=1,nyen + do i=1,nxen + xc(i,j)=cos(grid_latt(i,j)*deg2rad)*cos(grid_lont(i,j)*deg2rad) + yc(i,j)=cos(grid_latt(i,j)*deg2rad)*sin(grid_lont(i,j)*deg2rad) + zc(i,j)=sin(grid_latt(i,j)*deg2rad) + enddo + enddo + +! compute center as average x,y,z coordinates of corners of domain -- + + xcent=quarter*(xc(1,1)+xc(1,nyen)+xc(nxen,1)+xc(nxen,nyen)) + ycent=quarter*(yc(1,1)+yc(1,nyen)+yc(nxen,1)+yc(nxen,nyen)) + zcent=quarter*(zc(1,1)+zc(1,nyen)+zc(nxen,1)+zc(nxen,nyen)) + + rnorm=one/sqrt(xcent**2+ycent**2+zcent**2) + xcent=rnorm*xcent + ycent=rnorm*ycent + zcent=rnorm*zcent + centlat=asin(zcent)*rad2deg + centlon=atan2(ycent,xcent)*rad2deg + +!! compute new lats, lons + call rotate2deg(grid_lont,grid_latt,gcrlon,gcrlat, & + centlon,centlat,nxen,nyen) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! compute analysis A-grid lats, lons +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!--------------------------obtain analysis grid dimensions nxe,nye + nxe=nlon_ens + nye=nlat_ens + if(mype==0) print *,'nlat,nlon=nye,nxe= ',nlat_ens,nlon_ens + + allocate(rlat_in(nlat_ens,nlon_ens),rlon_in(nlat_ens,nlon_ens)) + allocate(region_lon_tmp(nlat_ens,nlon_ens),region_lat_tmp(nlat_ens,nlon_ens)) + region_lon_tmp=region_lon_ens*rad2deg + region_lat_tmp=region_lat_ens*rad2deg + call rotate2deg(region_lon_tmp,region_lat_tmp,rlon_in,rlat_in, & + centlon,centlat,nlat_ens,nlon_ens) + +!--------------------------obtain analysis grid spacing + dlat=(maxval(gcrlat)-minval(gcrlat))/(nyen-1) + dlon=(maxval(gcrlon)-minval(gcrlon))/(nxen-1) + + +!-----setup analysis A-grid from center of the domain +!--------------------compute all combinations of relative coordinates + + allocate(xbh_a(nxen),xbh_b(nxen,nyen),xa_a(nxe),xa_b(nxe)) + allocate(ybh_a(nyen),ybh_b(nxen,nyen),ya_a(nye),ya_b(nye)) + + nxh=nxen/2 + nyh=nyen/2 + + +!!!!!! fv3 rotated grid; not equal spacing, non_orthogonal !!!!!! + do j=1,nyen + jr=nyen+1-j + do i=1,nxen + ir=nxen+1-i + xbh_b(ir,jr)=gcrlon(i,j)/dlon + end do + end do + do j=1,nyen + jr=nyen+1-j + do i=1,nxen + ir=nxen+1-i + ybh_b(ir,jr)=gcrlat(i,j)/dlat + end do + end do + +!!!! define analysis A grid !!!!!!!!!!!!! + + index0=1 + do j=1,nxe + xa_a(j)= rlon_in(index0,j)/dlon + end do + do i=1,nye + ya_a(i)= rlat_in(i,index0)/dlat + end do + +!!!!!compute fv3 to A grid interpolation parameters !!!!!!!!! + allocate (fv3dxens(nxe,nye),fv3dx1ens(nxe,nye),fv3dyens(nxe,nye),fv3dy1ens(nxe,nye)) + allocate (fv3ixens(nxe,nye),fv3ixpens(nxe,nye),fv3jyens(nxe,nye),fv3jypens(nxe,nye)) + allocate(yy(nyen)) + +! iteration to find the fv3 grid cell + jb1=1 + ib1=1 + do j=1,nye + do i=1,nxe + do n=1,3 + gxa=xa_a(i) + if(gxa < xbh_b(1,jb1))then + gxa= 1 + else if(gxa > xbh_b(nxen,jb1))then + gxa= nxen + else + call grdcrd1(gxa,xbh_b(1,jb1),nxen,1) + endif + ib2=ib1 + ib1=gxa + do jj=1,nyen + yy(jj)=ybh_b(ib1,jj) + enddo + gya=ya_a(j) + if(gya < yy(1))then + gya= 1 + else if(gya > yy(nyen))then + gya= nyen + else + call grdcrd1(gya,yy,nyen,1) + endif + jb2=jb1 + jb1=gya + if(ib1+1 > nxen)then !this block( 6 lines) is copied from GSL gsi repository + ib1=ib1-1 + endif + if(jb1+1 > nyen)then + jb1=jb1-1 + endif + + if((ib1 == ib2) .and. (jb1 == jb2)) exit + if(n==3 ) then +!!!!!!! if not converge, find the nearest corner point + d(1)=(xa_a(i)-xbh_b(ib1,jb1))**2+(ya_a(j)-ybh_b(ib1,jb1))**2 + d(2)=(xa_a(i)-xbh_b(ib1+1,jb1))**2+(ya_a(j)-ybh_b(ib1+1,jb1))**2 + d(3)=(xa_a(i)-xbh_b(ib1,jb1+1))**2+(ya_a(j)-ybh_b(ib1,jb1+1))**2 + d(4)=(xa_a(i)-xbh_b(ib1+1,jb1+1))**2+(ya_a(j)-ybh_b(ib1+1,jb1+1))**2 + kk=1 + do k=2,4 + if(d(k) xa_a(nxe))then + gxa= nxe + else + call grdcrd1(gxa,xa_a,nxe,1) + endif + a3ixens(j,i)=int(gxa) + a3ixens(j,i)=min(max(1,a3ixens(j,i)),nxe) + a3dxens(j,i)=max(zero,min(one,gxa-a3ixens(j,i))) + a3dx1ens(j,i)=one-a3dxens(j,i) + a3ixpens(j,i)=min(nxe,a3ixens(j,i)+1) + end do + end do + + do i=1,nxen + do j=1,nyen + gya=ybh_b(i,j) + if(gya < ya_a(1))then + gya= 1 + else if(gya > ya_a(nye))then + gya= nye + else + call grdcrd1(gya,ya_a,nye,1) + endif + a3jyens(j,i)=int(gya) + a3jyens(j,i)=min(max(1,a3jyens(j,i)),nye) + a3dyens(j,i)=max(zero,min(one,gya-a3jyens(j,i))) + a3dy1ens(j,i)=one-a3dyens(j,i) + a3jypens(j,i)=min(nye,a3jyens(j,i)+1) + end do + end do + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! find coefficients for wind conversion btw FV3 & earth +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + allocate (canguens(nxen,nyen+1),sanguens(nxen,nyen+1),cangvens(nxen+1,nyen),sangvens(nxen+1,nyen)) + +! 1. compute x,y,z at cell cornor from grid_lon, grid_lat + + do j=1,nyen+1 + do i=1,nxen+1 + x(i,j)=cos(grid_lat(i,j)*deg2rad)*cos(grid_lon(i,j)*deg2rad) + y(i,j)=cos(grid_lat(i,j)*deg2rad)*sin(grid_lon(i,j)*deg2rad) + z(i,j)=sin(grid_lat(i,j)*deg2rad) + enddo + enddo + +! 2 find angles to E-W and N-S for U edges + + sq180=180._r_kind**2 + do j=1,nyen+1 + do i=1,nxen +! center lat/lon of the edge + rlat=half*(grid_lat(i,j)+grid_lat(i+1,j)) +! rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) + diff=(grid_lon(i,j)-grid_lon(i+1,j))**2 + if(diff < sq180)then + rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) + else + rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)-360._r_kind) + endif +! vector to center of the edge + xr=cos(rlat*deg2rad)*cos(rlon*deg2rad) + yr=cos(rlat*deg2rad)*sin(rlon*deg2rad) + zr=sin(rlat*deg2rad) +! vector of the edge + xu= x(i+1,j)-x(i,j) + yu= y(i+1,j)-y(i,j) + zu= z(i+1,j)-z(i,j) +! find angle with cross product + uval=sqrt((xu**2+yu**2+zu**2)) + ewval=sqrt((xr**2+yr**2)) + nsval=sqrt((xr*zr)**2+(zr*yr)**2+(xr*xr+yr*yr)**2) + canguens(i,j)=(-yr*xu+xr*yu)/ewval/uval + sanguens(i,j)=(-xr*zr*xu-zr*yr*yu+(xr*xr+yr*yr)*zu) / nsval/uval + enddo + enddo + +! 3 find angles to E-W and N-S for V edges + do j=1,nyen + do i=1,nxen+1 + rlat=half*(grid_lat(i,j)+grid_lat(i,j+1)) +! rlon=half*(grid_lon(i,j)+grid_lon(i,j+1)) + diff=(grid_lon(i,j)-grid_lon(i+1,j))**2 + if(diff < sq180)then + rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)) + else + rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)-360._r_kind) + endif + xr=cos(rlat*deg2rad)*cos(rlon*deg2rad) + yr=cos(rlat*deg2rad)*sin(rlon*deg2rad) + zr=sin(rlat*deg2rad) + xv= x(i,j+1)-x(i,j) + yv= y(i,j+1)-y(i,j) + zv= z(i,j+1)-z(i,j) + vval=sqrt((xv**2+yv**2+zv**2)) + ewval=sqrt((xr**2+yr**2)) + nsval=sqrt((xr*zr)**2+(zr*yr)**2+(xr*xr+yr*yr)**2) + cangvens(i,j)=(-yr*xv+xr*yv)/ewval/vval + sangvens(i,j)=(-xr*zr*xv-zr*yr*yv+(xr*xr+yr*yr)*zv) / nsval/vval + enddo + enddo + deallocate( xc,yc,zc,gclat,gclon,gcrlat,gcrlon) + deallocate(rlat_in,rlon_in) +end subroutine definecoef_regular_grids + subroutine earthuv2fv3(u,v,nx,ny,u_out,v_out) !$$$ subprogram documentation block ! . . . . @@ -679,6 +1110,51 @@ subroutine fv3uv2earth(u,v,nx,ny,u_out,v_out) return end subroutine fv3uv2earth +subroutine fv3uv2earthens(u,v,nxen,nyen,u_out,v_out) +!$$$ subprogram documentation block +! . . . . +! subprogram: fv3uv2earthens +! prgmmr: wu 2017-06-15 +! +! abstract: project fv3 UV to earth UV and interpolate to the center of the +! cells +! +! program history log: +! +! +! input argument list: +! u,v - fv3 winds on the cell boundaries +! nx,ny - dimensions +! +! output argument list: +! u_out,v_out - output earth wind components at center of the cell +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use kinds, only: r_kind,i_kind + use constants, only: half + implicit none + + integer(i_kind), intent(in ) :: nxen,nyen ! fv3 tile x- and y-dimensions + real(r_kind),intent(in ) :: u(nxen,nyen+1),v(nxen+1,nyen) + real(r_kind),intent( out) :: u_out(nxen,nyen),v_out(nxen,nyen) + integer(i_kind) i,j + + do j=1,nyen + do i=1,nxen + u_out(i,j)=half *((u(i,j)*sangvens(i,j)-v(i,j)*sanguens(i,j))/(canguens(i,j)*sangvens(i,j)-sanguens(i,j)*cangvens(i,j)) & + +(u(i,j+1)*sangvens(i+1,j)-v(i+1,j)*sanguens(i,j+1))/(canguens(i,j+1)*sangvens(i+1,j)-sanguens(i,j+1)*cangvens(i+1,j))) + v_out(i,j)=half *((u(i,j)*cangvens(i,j)-v(i,j)*canguens(i,j))/(sanguens(i,j)*cangvens(i,j)-canguens(i,j)*sangvens(i,j)) & + +(u(i,j+1)*cangvens(i+1,j)-v(i+1,j)*canguens(i,j+1))/(sanguens(i,j+1)*cangvens(i+1,j)-canguens(i,j+1)*sangvens(i+1,j))) + end do + end do + return +end subroutine fv3uv2earthens + subroutine fv3_h_to_ll(b_in,a,nb,mb,na,ma,rev_flg) !$$$ subprogram documentation block ! . . . . @@ -753,6 +1229,86 @@ subroutine fv3_h_to_ll(b_in,a,nb,mb,na,ma,rev_flg) return end subroutine fv3_h_to_ll +subroutine fv3_h_to_ll_ens(b_in,a,nb,mb,na,ma,rev_flg) +!$$$ subprogram documentation block +! . . . . +! subprogram: fv3_h_to_ll +! prgmmr: wu 2017-05-30 +! +! abstract: interpolate from rotated fv3 grid to A grid. +! Interpolation choices 1)bilinear both ways +! 2)inverse-distance weighting average +! reverse E-W and N-S directions & reverse i,j for output array a(nlat,nlon) +! +! program history log: +! +! +! input argument list: +! mb,nb - fv3 dimensions +! ma,na - a dimensions +! b - input variable b +! xb,yb - b array x and y coordinates +! xa,ya - a array coordinates (xa in xb units, ya in yb units) +! +! output argument list: +! a - output interpolated array +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + use mpimod, only: mype + use constants, only: zero,one + implicit none + + integer(i_kind),intent(in ) :: mb,nb,ma,na + real(r_kind) ,intent(in ) :: b_in(nb,mb) + logical ,intent(in ) :: rev_flg + real(r_kind) ,intent( out) :: a(ma,na) + + integer(i_kind) i,j,ir,jr,mbp,nbp + real(r_kind) b(nb,mb) + + mbp=mb+1 + nbp=nb+1 + bilinear=.false. + if(rev_flg) then +!!!!!!!!! reverse E-W and N-S + do j=1,mb + jr=mbp-j + do i=1,nb + ir=nbp-i + b(ir,jr)=b_in(i,j) + end do + end do + else + b(:,:)=b_in(:,:) + endif +!!!!!!!!! interpolate to A grid & reverse ij for array a(lat,lon) + if(bilinear)then ! bilinear interpolation + do j=1,ma + do i=1,na + a(j,i)=fv3dx1ens(i,j)*(fv3dy1ens(i,j)*b(fv3ixens(i,j),fv3jyens(i,j)) & + +fv3dyens(i,j)*b(fv3ixens(i,j),fv3jypens(i,j))) & + +fv3dxens(i,j)*(fv3dy1ens(i,j)*b(fv3ixpens(i,j),fv3jyens(i,j)) & + +fv3dyens(i,j)*b(fv3ixpens(i,j),fv3jypens(i,j))) + end do + end do + else ! inverse-distance weighting average + do j=1,ma + do i=1,na + a(j,i)=fv3dxens(i,j)*b(fv3ixens(i,j),fv3jyens(i,j)) & + +fv3dyens(i,j)*b(fv3ixens(i,j),fv3jypens(i,j)) & + +fv3dx1ens(i,j)*b(fv3ixpens(i,j),fv3jyens(i,j)) & + +fv3dy1ens(i,j)*b(fv3ixpens(i,j),fv3jypens(i,j)) + end do + end do + endif + + return +end subroutine fv3_h_to_ll_ens + subroutine fv3_ll_to_h(a,b,nxa,nya,nxb,nyb,rev_flg) !$$$ subprogram documentation block ! . . . .