module gsi_rfv3io_mod !$$$ module documentation block ! . . . . ! module: gsi_rfv3io_mod ! prgmmr: ! ! abstract: IO routines for regional FV3 ! ! program history log: ! 2017-03-08 parrish - create new module gsi_rfv3io_mod, starting from ! gsi_nemsio_mod as a pattern. ! 2017-10-10 wu - setup A grid and interpolation coeff in generate_anl_grid ! 2018-02-22 wu - add subroutines for read/write fv3_ncdf ! 2019 ting - modifications for use for ensemble IO and cold start files ! 2019-03-13 CAPS(C. Tong) - Port direct radar DA capabilities. ! 2021-11-01 lei - modify for fv3-lam parallel IO ! 2022-01-07 Hu - add code to readi/write subdomain restart files. ! This function is needed when fv3 model sets ! io_layout(2)>1 ! subroutines included: ! sub gsi_rfv3io_get_grid_specs ! sub read_fv3_files ! sub read_fv3_netcdf_guess ! sub gsi_fv3ncdf2d_read ! sub gsi_fv3ncdf_read ! sub gsi_fv3ncdf_readuv ! sub wrfv3_netcdf ! sub gsi_fv3ncdf_writeuv ! sub gsi_fv3ncdf_writeps ! sub gsi_fv3ncdf_write ! sub gsi_fv3ncdf_write_v1 ! sub check ! ! variable definitions: ! ! attributes: ! langauge: f90 ! machine: ! !$$$ end documentation block use kinds, only: r_kind,i_kind use gridmod, only: nlon_regional,nlat_regional use constants, only:max_varname_length use gsi_bundlemod, only : gsi_bundle use general_sub2grid_mod, only: sub2grid_info use gridmod, only: fv3_io_layout_y implicit none public type_fv3regfilenameg public bg_fv3regfilenameg public fv3sar_bg_opt ! directory names (hardwired for now) type type_fv3regfilenameg character(len=:),allocatable :: grid_spec !='fv3_grid_spec' character(len=:),allocatable :: ak_bk !='fv3_akbk' character(len=:),allocatable :: dynvars !='fv3_dynvars' character(len=:),allocatable :: tracers !='fv3_tracer' character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: couplerres!='coupler.res' contains procedure , pass(this):: init=>fv3regfilename_init end type type_fv3regfilenameg integer(i_kind):: fv3sar_bg_opt=0 type(type_fv3regfilenameg):: bg_fv3regfilenameg integer(i_kind) nx,ny,nz integer(i_kind),dimension(:),allocatable :: ny_layout_len,ny_layout_b,ny_layout_e 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(:) integer(i_kind),allocatable:: ijnz(:),displsz_g(:) real(r_kind),dimension(:,: ),allocatable:: ges_ps_bg real(r_kind),dimension(:,: ),allocatable:: ges_ps_inc real(r_kind),dimension(:,:,: ),allocatable:: ges_delp_bg type(sub2grid_info) :: grd_fv3lam_dynvar_ionouv type(sub2grid_info) :: grd_fv3lam_tracer_ionouv type(sub2grid_info) :: grd_fv3lam_uv integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8 character(len=max_varname_length), dimension(ndynvarslist), parameter :: & vardynvars = [character(len=max_varname_length) :: & "u","v","u_w","u_s","v_w","v_s","t","tv","tsen","w","delp","ps","delzinc"] character(len=max_varname_length), dimension(ntracerslist), parameter :: & vartracers = [character(len=max_varname_length) :: & 'q','oz','ql','qi','qr','qs','qg','qnr'] character(len=max_varname_length), dimension(15), parameter :: & varfv3name = [character(len=max_varname_length) :: & 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ps','DZ'], & vgsiname = [character(len=max_varname_length) :: & 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','ps','delzinc'] character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d ! set default to private private ! set subroutines to public public :: gsi_rfv3io_get_grid_specs public :: gsi_fv3ncdf_read public :: gsi_fv3ncdf_read_v1 public :: gsi_fv3ncdf_readuv public :: gsi_fv3ncdf_readuv_v1 public :: read_fv3_files public :: read_fv3_netcdf_guess public :: wrfv3_netcdf public :: gsi_fv3ncdf2d_read_v1 public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv public :: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv integer(i_kind) mype_u,mype_v,mype_t,mype_q,mype_p,mype_delz,mype_oz,mype_ql integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc parameter( & k_f10m =1, & !fact10 k_stype=2, & !soil_type k_vfrac=3, & !veg_frac k_vtype=4, & !veg_type k_zorl =5, & !sfc_rough k_tsea =6, & !sfct ? k_snwdph=7, & !sno ? k_stc =8, & !soil_temp k_smc =9, & !soil_moi k_slmsk=10, & !isli k_orog =11, & !terrain n2d=11 ) logical :: grid_reverse_flag character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars2d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracermetvars2d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_dynvar_nouv !to define names in gsibundle character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_tracer_nouv !to define names in gsibundle type(gsi_bundle):: gsibundle_fv3lam_dynvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracer_nouv contains subroutine fv3regfilename_init(this,grid_spec_input,ak_bk_input,dynvars_input, & tracers_input,sfcdata_input,couplerres_input) implicit None class(type_fv3regfilenameg),intent(inout):: this character(*),optional :: grid_spec_input,ak_bk_input,dynvars_input, & tracers_input,sfcdata_input,couplerres_input if(present(grid_spec_input))then this%grid_spec=grid_spec_input else this%grid_spec='fv3_grid_spec' endif if(present(ak_bk_input))then this%ak_bk=ak_bk_input else this%ak_bk='fv3_ak_bk' endif if(present(dynvars_input))then this%dynvars=dynvars_input else this%dynvars='fv3_dynvars' endif if(present(tracers_input))then this%tracers=tracers_input else this%tracers='fv3_tracer' endif if(present(sfcdata_input))then this%sfcdata=sfcdata_input else this%sfcdata='fv3_sfcdata' endif if(present(couplerres_input))then this%couplerres=couplerres_input else this%couplerres='coupler.res' endif end subroutine fv3regfilename_init subroutine gsi_rfv3io_get_grid_specs(fv3filenamegin,ierr) !$$$ subprogram documentation block ! . . . . ! subprogram: 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: generate_anl_grid use gridmod, only:nsig,regional_time,regional_fhr,regional_fmin,aeta1_ll,aeta2_ll use gridmod, only:nlon_regional,nlat_regional,eta1_ll,eta2_ll 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 integer(i_kind) gfile_grid_spec type (type_fv3regfilenameg) :: fv3filenamegin character(:),allocatable :: grid_spec character(:),allocatable :: ak_bk character(len=:),allocatable :: coupler_res_filenam integer(i_kind),intent( out) :: ierr integer(i_kind) i,k,ndimensions,iret,nvariables,nattributes,unlimiteddimid integer(i_kind) len,gfile_loc character(len=max_varname_length) :: name integer(i_kind) myear,mmonth,mday,mhour,mminute,msecond real(r_kind),allocatable:: abk_fv3(:) integer(i_kind) imiddle,jmiddle ! if fv3_io_layout_y > 1 integer(i_kind) :: nio,nylen integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout coupler_res_filenam=fv3filenamegin%couplerres grid_spec=fv3filenamegin%grid_spec ak_bk=fv3filenamegin%ak_bk !!!!! set regional_time open(24,file=trim(coupler_res_filenam),form='formatted') read(24,*) read(24,*) read(24,*)myear,mmonth,mday,mhour,mminute,msecond close(24) if(mype==0) write(6,*)' myear,mmonth,mday,mhour,mminute,msecond=', myear,mmonth,mday,mhour,mminute,msecond regional_time(1)=myear regional_time(2)=mmonth regional_time(3)=mday regional_time(4)=mhour regional_time(5)=mminute regional_time(6)=msecond regional_fhr=zero ! forecast hour set zero for now regional_fmin=zero ! forecast min set zero for now !!!!!!!!!! grid_spec !!!!!!!!!!!!!!! ierr=0 iret=nf90_open(trim(grid_spec),nf90_nowrite,gfile_grid_spec) if(iret/=nf90_noerr) then write(6,*)' gsi_rfv3io_get_grid_specs: problem opening ',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') nx=len if(trim(name)=='grid_yt') ny=len enddo nlon_regional=nx nlat_regional=ny allocate(ny_layout_len(0:fv3_io_layout_y-1)) allocate(ny_layout_b(0:fv3_io_layout_y-1)) allocate(ny_layout_e(0:fv3_io_layout_y-1)) ny_layout_len=ny ny_layout_b=0 ny_layout_e=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_len(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_b(nio)=nylen + 1 nylen=nylen+ny_layout_len(nio) ny_layout_e(nio)=nylen enddo endif if(mype==0)write(6,*),'nx,ny=',nx,ny if(mype==0)write(6,*),'ny_layout_len=',ny_layout_len if(mype==0)write(6,*),'ny_layout_b=',ny_layout_b if(mype==0)write(6,*),'ny_layout_e=',ny_layout_e !!! get nx,ny,grid_lon,grid_lont,grid_lat,grid_latt,nz,ak,bk allocate(grid_lat(nx+1,ny+1)) allocate(grid_lon(nx+1,ny+1)) allocate(grid_latt(nx,ny)) allocate(grid_lont(nx,ny)) 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 ! ! 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=nx/2 jmiddle=ny/2 if( (grid_latt(imiddle,1) < grid_latt(imiddle,ny)) .and. & (grid_lont(1,jmiddle) < grid_lont(nx,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,nx,ny,1) call reverse_grid_r(grid_latt,nx,ny,1) call reverse_grid_r(grid_lon,nx+1,ny+1,1) call reverse_grid_r(grid_lat,nx+1,ny+1,1) endif iret=nf90_close(gfile_loc) iret=nf90_open(ak_bk,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)'gsi_rfv3io_get_grid_specs: problem opening ',trim(ak_bk),', Status = ',iret ierr=1 return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) if(trim(name)=='xaxis_1') nz=len enddo if(mype==0)write(6,'(" nz=",i5)') nz nsig=nz-1 !!! 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)) do k=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,k,name,len) if(trim(name)=='ak'.or.trim(name)=='AK') then iret=nf90_get_var(gfile_loc,k,abk_fv3) do i=1,nz ak(i)=abk_fv3(nz+1-i) enddo endif if(trim(name)=='bk'.or.trim(name)=='BK') then iret=nf90_get_var(gfile_loc,k,abk_fv3) do i=1,nz bk(i)=abk_fv3(nz+1-i) enddo endif enddo iret=nf90_close(gfile_loc) !!!!! change unit of ak do i=1,nsig+1 eta1_ll(i)=ak(i)*0.001_r_kind eta2_ll(i)=bk(i) enddo do i=1,nsig aeta1_ll(i)=half*(ak(i)+ak(i+1))*0.001_r_kind aeta2_ll(i)=half*(bk(i)+bk(i+1)) enddo if(mype==0)then do i=1,nz write(6,'(" ak,bk(",i3,") = ",2f17.6)') i,ak(i),bk(i) enddo endif !!!!!!! setup A grid and interpolation/rotation coeff. call generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt) deallocate (grid_lon,grid_lat,grid_lont,grid_latt) deallocate (ak,bk,abk_fv3) return end subroutine gsi_rfv3io_get_grid_specs subroutine read_fv3_files(mype) !$$$ subprogram documentation block ! . . . . ! subprogram: read_fv3_files ! prgmmr: wu org: np22 date: 2017-10-10 ! ! abstract: read in from fv3 files and figure out available time levels ! of background fields starting from read_files as a pattern ! temporary setup for one first guess file ! program history log: ! ! input argument list: ! mype - pe number ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,ierror,mpi_rtype,npe use guess_grids, only: nfldsig,nfldsfc,ntguessig,ntguessfc,& ifilesig,ifilesfc,hrdifsig,hrdifsfc,create_gesfinfo use guess_grids, only: hrdifsig_all,hrdifsfc_all use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,nhr_assimilation use constants, only: zero,one,r60inv use obsmod, only: iadate,time_offset use gridmod, only:regional_time implicit none ! Declare passed variables integer(i_kind),intent(in ) :: mype ! Declare local parameters real(r_kind),parameter:: r0_001=0.001_r_kind ! Declare local variables logical(4) fexist character(6) filename integer(i_kind) in_unit integer(i_kind) i,j,iwan,npem1 integer(i_kind) nhr_half integer(i_kind) nminanl,nmings,nming2,ndiff,isecond integer(i_kind),dimension(4):: idateg integer(i_kind),dimension(5):: idate5 real(r_kind) hourg,temp,t4dv real(r_kind),dimension(202,2):: time_ges !----------------------------------------------------------------------------- ! Start read_nems_nmmb_files here. nhr_half=nhr_assimilation/2 if(nhr_half*2 < nhr_assimilation) nhr_half=nhr_half+1 npem1=npe-1 do i=1,202 time_ges(i,1) = 999_r_kind time_ges(i,2) = 999_r_kind end do ! Let a single task query the guess files. if(mype==npem1) then ! Convert analysis time to minutes relative to fixed date call w3fs21(iadate,nminanl) write(6,*)'READ_netcdf_fv3_FILES: analysis date,minutes ',iadate,nminanl ! Check for consistency of times from sigma guess files. in_unit=15 iwan=0 !WWWWWW setup for one first guess file for now ! do i=0,9 !place holder for FGAT i=3 !wwww read in from the external file directly, no internal files sigfxx for FV3 idate5(1)= regional_time(1) idate5(2)= regional_time(2) idate5(3)= regional_time(3) idate5(4)= regional_time(4) idate5(5)= regional_time(5) isecond = regional_time(6) hourg = zero ! forcast hour call w3fs21(idate5,nmings) nming2=nmings+60*hourg write(6,*)'READ_netcdf_fv3_FILES: sigma guess file, nming2 ',hourg,idate5,nming2 t4dv=real((nming2-iwinbgn),r_kind)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) then write(6,*)'ges file not in time range, t4dv=',t4dv ! cycle ! place holder for FGAT endif else ndiff=nming2-nminanl !for test with the 3 hr files with FGAT if(abs(ndiff) > 60*nhr_half ) then write(6,*)'ges file not in time range, ndiff=',ndiff ! cycle ! place holder for FGAT endif endif iwan=iwan+1 time_ges(iwan,1) =real((nming2-iwinbgn),r_kind)*r60inv time_ges(iwan+100,1)=i+r0_001 ! end do ! i !place holder for FGAT time_ges(201,1)=one time_ges(202,1)=one if(iwan > 1)then do i=1,iwan do j=i+1,iwan if(time_ges(j,1) < time_ges(i,1))then temp=time_ges(i+100,1) time_ges(i+100,1)=time_ges(j+100,1) time_ges(j+100,1)=temp temp=time_ges(i,1) time_ges(i,1)=time_ges(j,1) time_ges(j,1)=temp end if end do if(abs(time_ges(i,1)-time_offset) < r0_001)time_ges(202,1) = i end do end if time_ges(201,1) = iwan+r0_001 ! Check for consistency of times from surface guess files. iwan=0 do i=0,99 write(filename,200)i 200 format('sfcf',i2.2) inquire(file=filename,exist=fexist) if(fexist)then idateg(4)=iadate(1); idateg(2)=iadate(2) idateg(3)=iadate(3); idateg(1)=iadate(4) hourg = zero idate5(1)=idateg(4); idate5(2)=idateg(2) idate5(3)=idateg(3); idate5(4)=idateg(1); idate5(5)=0 call w3fs21(idate5,nmings) nming2=nmings+60*hourg write(6,*)'READ_netcdf_fv3_FILES: surface guess file, nming2 ',hourg,idateg,nming2 ndiff=nming2-nminanl if(abs(ndiff) > 60*nhr_half ) then write(6,*)'ges file not in time range, ndiff=',ndiff ! cycle ! place holder for FGAT endif iwan=iwan+1 time_ges(iwan,2) =real((nming2-iwinbgn),r_kind)*r60inv time_ges(iwan+100,2)=i+r0_001 end if if(iwan==1) exit end do time_ges(201,2)=one time_ges(202,2)=one if(iwan > 1)then do i=1,iwan do j=i+1,iwan if(time_ges(j,2) < time_ges(i,2))then temp=time_ges(i+100,2) time_ges(i+100,2)=time_ges(j+100,2) time_ges(j+100,2)=temp temp=time_ges(i,2) time_ges(i,2)=time_ges(j,2) time_ges(j,2)=temp end if end do if(abs(time_ges(i,2)-time_offset) < r0_001)time_ges(202,2) = i end do end if time_ges(201,2) = iwan+r0_001 end if ! Broadcast guess file information to all tasks call mpi_bcast(time_ges,404,mpi_rtype,npem1,mpi_comm_world,ierror) nfldsig = nint(time_ges(201,1)) !!nfldsfc = nint(time_ges(201,2)) nfldsfc = nfldsig ! Allocate space for guess information files call create_gesfinfo do i=1,nfldsig ifilesig(i) = -100 hrdifsig(i) = zero end do do i=1,nfldsfc ifilesfc(i) = -100 hrdifsfc(i) = zero end do ! Load time information for sigma guess field sinfo into output arrays ntguessig = nint(time_ges(202,1)) do i=1,nfldsig hrdifsig(i) = time_ges(i,1) ifilesig(i) = nint(time_ges(i+100,1)) hrdifsig_all(i) = hrdifsig(i) end do if(mype == 0) write(6,*)'READ_netcdf_fv3_FILES: sigma fcst files used in analysis : ',& (ifilesig(i),i=1,nfldsig),(hrdifsig(i),i=1,nfldsig),ntguessig ! Load time information for surface guess field info into output arrays ntguessfc = nint(time_ges(202,2)) do i=1,nfldsfc hrdifsfc(i) = time_ges(i,2) ifilesfc(i) = nint(time_ges(i+100,2)) hrdifsfc_all(i) = hrdifsfc(i) end do ! Below is a temporary fix. The nems_nmmb regional mode does not have a ! surface ! file. Instead the surface fields are passed through the atmospheric guess ! file. Without a separate surface file the code above sets ntguessig and ! nfldsig to zero. This causes problems later in the code when arrays for ! the surface fields are allocated --> one of the array dimensions is nfldsfc ! and it will be zero. This portion of the code should be rewritten, but ! until ! it is, the fix below gets around the above mentioned problem. ntguessfc = ntguessig !!nfldsfc = nfldsig do i=1,nfldsfc hrdifsfc(i) = hrdifsig(i) ifilesfc(i) = ifilesig(i) hrdifsfc_all(i) = hrdifsfc(i) end do if(mype == 0) write(6,*)'READ_nems_nmb_FILES: surface fcst files used in analysis: ',& (ifilesfc(i),i=1,nfldsfc),(hrdifsfc(i),i=1,nfldsfc),ntguessfc ! End of routine return end subroutine read_fv3_files subroutine read_fv3_netcdf_guess(fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: read_fv3_netcdf_guess read fv3 interface file ! prgmmr: wu org: np22 date: 2017-07-06 ! ! abstract: read guess for FV3 regional model ! program history log: ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: npe use mpimod, only: mpi_comm_world use guess_grids, only:ges_prsi use gridmod, only: lat2,lon2,nsig,ijn,eta1_ll,eta2_ll,ijn_s use constants, only: one,fv use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundleinquire, gsi_bundlegetpointer use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundledestroy use general_sub2grid_mod, only: general_sub2grid_create_info use mpeu_util, only: die use guess_grids, only: ntguessig use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval use directDA_radaruse_mod, only: l_use_dbz_directDA use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval use gridmod,only: regional use gridmod,only: l_reg_update_hydro_delz use gridmod, only: grd_a use mpimod, only: mype use gsi_metguess_mod, only: gsi_metguess_get implicit none type (type_fv3regfilenameg),intent (in) :: fv3filenamegin character(len=24),parameter :: myname = 'read_fv3_netcdf_guess' integer(i_kind) k,i,j integer(i_kind) it,ier,istatus real(r_kind),dimension(:,:),pointer::ges_ps=>NULL() real(r_kind),dimension(:,:),pointer::ges_ps_readin=>NULL() real(r_kind),dimension(:,:),pointer::ges_z=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_u=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_v=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_q=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_oz=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_tsen_readin=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_iqr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qs=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qg=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qnr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_w=>NULL() character(len=max_varname_length) :: vartem="" character(len=64),dimension(:,:),allocatable:: names !to be same as in the grid the dummy sub2grid_info character(len=64),dimension(:,:),allocatable:: uvnames integer(i_kind),dimension(:,:),allocatable:: lnames integer(i_kind),dimension(:,:),allocatable:: uvlnames integer(i_kind):: inner_vars,numfields integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer integer(r_kind):: iuv,ndynvario3d,ntracerio3d !clt this block is still maintained for they would be needed for a certain 2d fields IO mype_2d=mod(1,npe) allocate(ijns(npe),ijns2d(npe),ijnz(npe) ) allocate(displss(npe),displss2d(npe),displsz_g(npe) ) do i=1,npe ijns(i)=ijn_s(i)*nsig ijnz(i)=ijn(i)*nsig ijns2d(i)=ijn_s(i)*n2d enddo displss(1)=0 displsz_g(1)=0 displss2d(1)=0 do i=2,npe displss(i)=displss(i-1)+ ijns(i-1) displsz_g(i)=displsz_g(i-1)+ ijnz(i-1) displss2d(i)=displss2d(i-1)+ ijns2d(i-1) enddo it=ntguessig allocate( name_metvars2d(GSI_MetGuess_Bundle(it)%n2d)) allocate( name_metvars3d(GSI_MetGuess_Bundle(it)%n3d)) call gsi_bundleinquire (GSI_MetGuess_Bundle(it),'shortnames::2d', name_metvars2d,istatus) call gsi_bundleinquire (GSI_MetGuess_Bundle(it),'shortnames::3d', name_metvars3d,istatus) if(mype == 0) then do i=1,GSI_MetGuess_Bundle(it)%n2d write(6,*)'metvardeb333-2d name ', trim(name_metvars2d(i)) enddo do i=1,GSI_MetGuess_Bundle(it)%n3d write(6,*)'metvardeb333-3d name ', trim(name_metvars3d(i)) enddo endif !here a strict requirment for the names of "u" and "v" is rquired !maybe more flexibile procedure would relax it iuv=0 ndynvario3d=0 ntracerio3d=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(trim(vartem)=='u'.or.trim(vartem)=='v') then iuv=iuv+1 else if(.not.(trim(vartem)=='iqr')) then if (ifindstrloc(vardynvars,trim(vartem))> 0) then ndynvario3d=ndynvario3d+1 else if (ifindstrloc(vartracers,trim(vartem))> 0) then ntracerio3d=ntracerio3d+1 else write(6,*)'the metvarname1 ',trim(vartem),' has not been considered yet, stop' call stop2(333) endif endif ! iqr is the inital qr, need not to be in IO endif end do if (iuv /= 2.or. ndynvario3d<=0.or.ntracerio3d<=0 ) then write(6,*)"the set up for met variable is not as expected, abort" call stop2(222) endif if (fv3sar_bg_opt == 0.and.ifindstrloc(name_metvars3d,'delp') <= 0)then ndynvario3d=ndynvario3d+1 ! for delp endif if (fv3sar_bg_opt == 1.and.ifindstrloc(name_metvars3d,'delp') > 0)then ndynvario3d=ndynvario3d-1 ! delp should not be used in IO endif if (l_reg_update_hydro_delz.and.fv3sar_bg_opt==0) ndynvario3d=ndynvario3d+1 ! for delzinc allocate(fv3lam_io_dynmetvars3d_nouv(ndynvario3d)) allocate(fv3lam_io_tracermetvars3d_nouv(ntracerio3d)) jdynvar=0 jtracer=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(.not.(trim(vartem)=='u'.or.trim(vartem)=='v'.or.trim(vartem)=='iqr')) then if(trim(vartem)=="tv" ) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)="tsen" else if (ifindstrloc(vardynvars,trim(vartem)) > 0) then if(.not.(fv3sar_bg_opt==1.and.trim(vartem)=="delp")) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)=trim(vartem) endif else if (ifindstrloc(vartracers,trim(vartem)) > 0) then jtracer=jtracer+1 fv3lam_io_tracermetvars3d_nouv(jtracer)=trim(vartem) else write(6,*)'the metvarname ',vartem,' is not expected, stop' call flush(6) call stop2(333) endif endif enddo if(fv3sar_bg_opt == 0.and.ifindstrloc( fv3lam_io_dynmetvars3d_nouv(1:jdynvar),"delp")<= 0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)="delp" endif if (l_reg_update_hydro_delz.and.fv3sar_bg_opt==0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)="delzinc" endif if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d ) then write(6,*)'ndynvario3d is not as expected, stop' call flush(6) call stop2(333) endif if(mype == 0) then write(6,*) ' fv3lam_io_dynmetvars3d_nouv is ',(trim(fv3lam_io_dynmetvars3d_nouv(i)),i=1,ndynvario3d) write(6,*) ' fv3lam_io_tracermevars3d_nouv is ',(trim(fv3lam_io_tracermetvars3d_nouv(i)),i=1,ntracerio3d) endif ndynvario2d=0 ntracerio2d=0 do i=1,size(name_metvars2d) vartem=trim(name_metvars2d(i)) if(.not. (trim(vartem)=='ps'.and.fv3sar_bg_opt==0)) then if (ifindstrloc(vardynvars,trim(vartem))> 0) then ndynvario2d=ndynvario2d+1 else if (ifindstrloc(vartracers,trim(vartem)) > 0) then ntracerio2d=ntracerio2d+1 else if(trim(vartem)=='z') then write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' else write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop' call stop2(333) endif endif end do if (ndynvario2d > 0) then allocate(fv3lam_io_dynmetvars2d_nouv(ndynvario2d)) endif if (ntracerio2d > 0) allocate(fv3lam_io_tracermetvars2d_nouv(ntracerio2d)) jdynvar=0 jtracer=0 do i=1,size(name_metvars2d) vartem=trim(name_metvars2d(i)) if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z"))) then !z is treated separately if (ifindstrloc(vardynvars,trim(vartem)) > 0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem) else if (ifindstrloc(vartracers,trim(vartem)) > 0) then jtracer=jtracer+1 fv3lam_io_tracermetvars2d_nouv(jdynvar)=trim(vartem) else write(6,*)'the metvarname3 ',trim(vartem),' has not been considered yet, stop' call stop2(333) endif endif end do if(mype == 0) then if (allocated(fv3lam_io_dynmetvars2d_nouv)) & write(6,*)' fv3lam_io_dynmetvars2d_nouv is ',(trim(fv3lam_io_dynmetvars2d_nouv(i)), i=1,ndynvario2d) if (allocated(fv3lam_io_tracermetvars2d_nouv))& write(6,*)'fv3lam_io_dynmetvars2d_nouv is ',(trim(fv3lam_io_dynmetvars2d_nouv(i)),i=1,ntracerio3d) endif if (allocated(fv3lam_io_dynmetvars2d_nouv) ) then call gsi_bundlecreate(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_dynvar_nouv',istatus,& names2d=fv3lam_io_dynmetvars2d_nouv,names3d=fv3lam_io_dynmetvars3d_nouv) ndynvario2d=size(fv3lam_io_dynmetvars2d_nouv) else call gsi_bundlecreate(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_dynvar_nouv',istatus, & names3d=fv3lam_io_dynmetvars3d_nouv) ndynvario2d=0 endif if (allocated(fv3lam_io_tracermetvars2d_nouv) ) then call gsi_bundlecreate(gsibundle_fv3lam_tracer_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_tracer_varnouv',istatus,& names2d=fv3lam_io_tracermetvars2d_nouv,names3d=fv3lam_io_tracermetvars2d_nouv) ntracerio2d=size(fv3lam_io_tracermetvars2d_nouv) else call gsi_bundlecreate(gsibundle_fv3lam_tracer_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_tracer_nouv',istatus, & names3d=fv3lam_io_tracermetvars3d_nouv) ntracerio2d=0 endif inner_vars=1 numfields=inner_vars*(ndynvario3d*grd_a%nsig+ndynvario2d) allocate(lnames(1,numfields),names(1,numfields)) ilev=1 do i=1,ndynvario3d do k=1,grd_a%nsig lnames(1,ilev)=k names(1,ilev)=trim(fv3lam_io_dynmetvars3d_nouv(i)) ilev=ilev+1 enddo enddo do i=1,ndynvario2d lnames(1,ilev)=1 names(1,ilev)=trim(fv3lam_io_dynmetvars2d_nouv(i)) ilev=ilev+1 enddo call general_sub2grid_create_info(grd_fv3lam_dynvar_ionouv,inner_vars,grd_a%nlat,& grd_a%nlon,grd_a%nsig,numfields,regional,names=names,lnames=lnames) inner_vars=1 numfields=inner_vars*(ntracerio3d*grd_a%nsig+ntracerio2d) deallocate(lnames,names) allocate(lnames(1,numfields),names(1,numfields)) ilev=1 do i=1,ntracerio3d do k=1,grd_a%nsig lnames(1,ilev)=k names(1,ilev)=trim(fv3lam_io_tracermetvars3d_nouv(i)) ilev=ilev+1 enddo enddo do i=1,ntracerio2d lnames(1,ilev)=1 names(1,ilev)=trim(fv3lam_io_tracermetvars2d_nouv(i)) ilev=ilev+1 enddo call general_sub2grid_create_info(grd_fv3lam_tracer_ionouv,inner_vars,grd_a%nlat,& grd_a%nlon,grd_a%nsig,numfields,regional,names=names,lnames=lnames) inner_vars=2 numfields=grd_a%nsig allocate(uvlnames(inner_vars,numfields),uvnames(inner_vars,numfields)) do k=1,grd_a%nsig uvlnames(1,k)=k uvlnames(2,k)=k uvnames(1,k)='u' uvnames(2,k)='v' enddo call general_sub2grid_create_info(grd_fv3lam_uv,inner_vars,grd_a%nlat,& grd_a%nlon,grd_a%nsig,numfields,regional,names=uvnames,lnames=uvlnames) deallocate(lnames,names,uvlnames,uvnames) ier=0 call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ps' ,ges_ps ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'z' , ges_z ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tv' ,ges_tv ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'oz' ,ges_oz ,istatus );ier=ier+istatus if (l_use_dbz_directDA) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'iqr' ,ges_iqr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus );ier=ier+istatus end if if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 met-fields, ier =',ier) if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_readuv(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin) else call gsi_fv3ncdf_readuv_v1(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin) endif if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_read(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_read(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv,fv3filenamegin%tracers,fv3filenamegin) else call gsi_fv3ncdf_read_v1(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_read_v1(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv,fv3filenamegin%tracers,fv3filenamegin) endif if( fv3sar_bg_opt == 0) then call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'delp' ,ges_delp ,istatus );ier=ier+istatus if(istatus==0) ges_delp=ges_delp*0.001_r_kind endif call gsi_copy_bundle(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)) call gsi_copy_bundle(gsibundle_fv3lam_tracer_nouv,GSI_MetGuess_Bundle(it)) call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'tsen' ,ges_tsen_readin ,istatus );ier=ier+istatus !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig do j=1,lon2 do i=1,lat2 ges_tv(i,j,k)=ges_tsen_readin(i,j,k)*(one+fv*ges_q(i,j,k)) enddo enddo enddo if( fv3sar_bg_opt == 0) then allocate(ges_delp_bg(lat2,lon2,nsig)) allocate(ges_ps_bg(lat2,lon2)) ges_delp_bg=ges_delp ges_prsi(:,:,nsig+1,it)=eta1_ll(nsig+1) do i=nsig,1,-1 ges_prsi(:,:,i,it)=ges_delp(:,:,i)+ges_prsi(:,:,i+1,it) enddo ges_ps(:,:)=ges_prsi(:,:,1,it) ges_ps_bg=ges_ps else call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'ps' ,ges_ps_readin ,istatus );ier=ier+istatus ges_ps_readin=ges_ps_readin*0.001_r_kind !which is from ges_ps=ges_ps_readin ges_ps_bg=ges_ps ges_prsi(:,:,nsig+1,it)=eta1_ll(nsig+1) do k=1,nsig ges_prsi(:,:,k,it)=eta1_ll(k)+eta2_ll(k)*ges_ps enddo endif call gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) if (l_use_dbz_directDA ) then if( fv3sar_bg_opt == 0) then ges_iqr=ges_qr else write(6,*) "FV3 IO READ for 'fv3sar_bg_opt == 0' is only available for now in direct reflectivity DA" stop end if call convert_qx_to_cvpqx(ges_qr, ges_qs, ges_qg, l_use_cvpqx, cvpqx_pval) ! convert Qx call convert_nx_to_cvpnx(ges_qnr, l_cvpnr, cvpnr_pval) ! convert Qnx end if end subroutine read_fv3_netcdf_guess subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf2d_read ! prgmmr: wu w org: np22 date: 2017-10-17 ! ! abstract: read in 2d fields from fv3_sfcdata file in mype_2d ! Scatter the field to each PE ! program history log: ! input argument list: ! it - time index for 2d fields ! ! output argument list: ! ges_z - surface elevation ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype use guess_grids, only: fact10,soil_type,veg_frac,veg_type,sfc_rough, & sfct,sno,soil_temp,soil_moi,isli use gridmod, only: lat2,lon2,itotsub,ijn_s use general_commvars_mod, only: ltosi_s,ltosj_s 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 mod_fv3_lola, only: fv3_h_to_ll,nxa,nya use constants, only: grav implicit none integer(i_kind),intent(in) :: it real(r_kind),intent(in),dimension(:,:),pointer::ges_z type (type_fv3regfilenameg),intent(in) :: fv3filenamegin character(len=max_varname_length) :: name integer(i_kind),allocatable,dimension(:):: dim_id,dim real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: a real(r_kind),allocatable,dimension(:,:,:):: sfcn2d real(r_kind),allocatable,dimension(:,:,:):: sfc real(r_kind),allocatable,dimension(:,:):: sfc1 integer(i_kind) iret,gfile_loc,i,k,len,ndim integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid integer(i_kind) kk,n,ns,j,ii,jj,mm1 character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: dynvars !='fv3_dynvars' ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: sfc_fulldomain integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout sfcdata= fv3filenamegin%sfcdata dynvars= fv3filenamegin%dynvars mm1=mype+1 allocate(a(nya,nxa)) allocate(work(itotsub*n2d)) allocate( sfcn2d(lat2,lon2,n2d)) if(mype==mype_2d ) then allocate(sfc_fulldomain(nx,ny)) 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(sfcdata),'.',nio iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio)) if(iret/=nf90_noerr) then write(6,*)' problem opening3 ',trim(filename_layout),', Status = ',iret return endif enddo gfile_loc=gfile_loc_layout(0) else iret=nf90_open(sfcdata,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening3 ',trim(sfcdata),', Status = ',iret return endif endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(dim(ndimensions)) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!! do i=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,i,name,len) if( trim(name)=='f10m'.or.trim(name)=='F10M' ) then k=k_f10m else if( trim(name)=='stype'.or.trim(name)=='STYPE' ) then k=k_stype else if( trim(name)=='vfrac'.or.trim(name)=='VFRAC' ) then k=k_vfrac else if( trim(name)=='vtype'.or.trim(name)=='VTYPE' ) then k=k_vtype else if( trim(name)=='zorl'.or.trim(name)=='ZORL' ) then k=k_zorl else if( trim(name)=='tsea'.or.trim(name)=='TSEA' ) then k=k_tsea else if( trim(name)=='sheleg'.or.trim(name)=='SHELEG' ) then k=k_snwdph else if( trim(name)=='stc'.or.trim(name)=='STC' ) then k=k_stc else if( trim(name)=='smc'.or.trim(name)=='SMC' ) then k=k_smc else if( trim(name)=='SLMSK'.or.trim(name)=='slmsk' ) then k=k_slmsk else cycle endif iret=nf90_inquire_variable(gfile_loc,i,ndims=ndim) if(ndim < 2) then write(*,*) "wrong dimension number ndim =",ndim call stop2(119) endif if(allocated(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 iret=nf90_inquire_variable(gfile_loc_layout(nio),i,dimids=dim_id) if(allocated(sfc )) deallocate(sfc ) if(dim(dim_id(1)) == nx .and. dim(dim_id(2))==ny_layout_len(nio)) then if(ndim >=3) then allocate(sfc(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) iret=nf90_get_var(gfile_loc_layout(nio),i,sfc) else if (ndim == 2) then allocate(sfc(dim(dim_id(1)),dim(dim_id(2)),1)) iret=nf90_get_var(gfile_loc_layout(nio),i,sfc(:,:,1)) endif else write(*,*) "Mismatch dimension in surfacei reading:",nx,ny_layout_len(nio),dim(dim_id(1)),dim(dim_id(2)) call stop2(119) endif sfc_fulldomain(:,ny_layout_b(nio):ny_layout_e(nio))=sfc(:,:,1) enddo else iret=nf90_inquire_variable(gfile_loc,i,dimids=dim_id) if(allocated(sfc )) deallocate(sfc ) if(dim(dim_id(1)) == nx .and. dim(dim_id(2))==ny) then if(ndim >=3) then !the block of 10 lines is compied from GSL gsi. allocate(sfc(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) iret=nf90_get_var(gfile_loc,i,sfc) else if (ndim == 2) then allocate(sfc(dim(dim_id(1)),dim(dim_id(2)),1)) iret=nf90_get_var(gfile_loc,i,sfc(:,:,1)) endif else write(*,*) "Mismatch dimension in surfacei reading:",nx,ny,dim(dim_id(1)),dim(dim_id(2)) call stop2(119) endif sfc_fulldomain(:,:)=sfc(:,:,1) endif call fv3_h_to_ll(sfc_fulldomain,a,nx,ny,nxa,nya,grid_reverse_flag) kk=0 do n=1,npe ns=displss2d(n)+(k-1)*ijn_s(n) do j=1,ijn_s(n) ns=ns+1 kk=kk+1 ii=ltosi_s(kk) jj=ltosj_s(kk) work(ns)=a(ii,jj) end do end do enddo ! i if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 iret=nf90_close(gfile_loc_layout(nio)) enddo deallocate (gfile_loc_layout) else iret=nf90_close(gfile_loc) endif !!!! read in orog from dynam !!!!!!!!!!!! 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(dynvars),'.',nio iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio)) if(iret/=nf90_noerr) then write(6,*)' problem opening4 ',trim(filename_layout),', Status =',iret return endif enddo gfile_loc=gfile_loc_layout(0) else iret=nf90_open(dynvars,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening4 ',trim(dynvars ),gfile_loc,', Status = ',iret return endif endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) if(allocated(dim )) deallocate(dim ) allocate(dim(ndimensions)) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo do k=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,k,name,len) if(trim(name)=='PHIS' .or. trim(name)=='phis' ) then iret=nf90_inquire_variable(gfile_loc,k,ndims=ndim) if(allocated(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 iret=nf90_inquire_variable(gfile_loc_layout(nio),k,dimids=dim_id) if(allocated(sfc1 )) deallocate(sfc1 ) allocate(sfc1(dim(dim_id(1)),dim(dim_id(2))) ) iret=nf90_get_var(gfile_loc_layout(nio),k,sfc1) sfc_fulldomain(:,ny_layout_b(nio):ny_layout_e(nio))=sfc1 enddo else iret=nf90_inquire_variable(gfile_loc,k,dimids=dim_id) allocate(sfc1(dim(dim_id(1)),dim(dim_id(2))) ) iret=nf90_get_var(gfile_loc,k,sfc1) sfc_fulldomain=sfc1 endif exit endif enddo ! k if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 iret=nf90_close(gfile_loc_layout(nio)) enddo deallocate(gfile_loc_layout) else iret=nf90_close(gfile_loc) endif k=k_orog call fv3_h_to_ll(sfc_fulldomain,a,nx,ny,nxa,nya,grid_reverse_flag) kk=0 do n=1,npe ns=displss2d(n)+(k-1)*ijn_s(n) do j=1,ijn_s(n) ns=ns+1 kk=kk+1 ii=ltosi_s(kk) jj=ltosj_s(kk) work(ns)=a(ii,jj) end do end do deallocate (dim_id,sfc,sfc1,dim) deallocate (sfc_fulldomain) endif ! mype !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,& sfcn2d,ijns2d(mm1),mpi_rtype,mype_2d,mpi_comm_world,ierror) deallocate ( work ) fact10(:,:,it)=sfcn2d(:,:,k_f10m) soil_type(:,:,it)=sfcn2d(:,:,k_stype) veg_frac(:,:,it)=sfcn2d(:,:,k_vfrac) veg_type(:,:,it)=sfcn2d(:,:,k_vtype) sfc_rough(:,:,it)=sfcn2d(:,:,k_zorl) sfct(:,:,it)=sfcn2d(:,:,k_tsea) sno(:,:,it)=sfcn2d(:,:,k_snwdph) soil_temp(:,:,it)=sfcn2d(:,:,k_stc) soil_moi(:,:,it)=sfcn2d(:,:,k_smc) ges_z(:,:)=sfcn2d(:,:,k_orog)/grav isli(:,:,it)=nint(sfcn2d(:,:,k_slmsk)) deallocate (sfcn2d,a) return end subroutine gsi_fv3ncdf2d_read subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv23ncdf2d_readv1 ! prgmmr: T. Lei date: 2019-03-28 ! modified from gsi_fv3ncdf_read and gsi_fv3ncdf2d_read ! ! abstract: read in a 2d field from a netcdf FV3 file in mype_io ! then scatter the field to each PE ! program history log: ! ! input argument list: ! filename - file name to read from ! varname - variable name to read in ! varname2 - variable name to read in ! mype_io - pe to read in the field ! ! output argument list: ! work_sub - output sub domain field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype use gridmod, only: lat2,lon2,nlat,nlon,itotsub,ijn_s,displs_s 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_inq_varid use netcdf, only: nf90_inquire_variable use mod_fv3_lola, only: fv3_h_to_ll use general_commvars_mod, only: ltosi_s,ltosj_s implicit none character(*) ,intent(in ) :: varname,varname2,filenamein real(r_kind) ,intent(out ) :: work_sub(lat2,lon2) integer(i_kind) ,intent(in ) :: mype_io real(r_kind),allocatable,dimension(:,:,:):: uu real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: a integer(i_kind) n,ndim integer(i_kind) gfile_loc,var_id,iret integer(i_kind) kk,j,mm1,ii,jj integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid mm1=mype+1 allocate (work(itotsub)) if(mype==mype_io ) then iret=nf90_open(filenamein,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf2d_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret write(6,*)' gsi_fv3ncdf2d_read_v1: problem opening with varnam ',trim(varname) return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(a(nlat,nlon)) iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) if(iret/=nf90_noerr) then iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname2)),var_id) if(iret/=nf90_noerr) then write(6,*)' wrong to get var_id ',var_id endif endif iret=nf90_inquire_variable(gfile_loc,var_id,ndims=ndim) if(allocated(uu )) deallocate(uu ) allocate(uu(nx,ny,1)) iret=nf90_get_var(gfile_loc,var_id,uu) call fv3_h_to_ll(uu(:,:,1),a,nx,ny,nlon,nlat,grid_reverse_flag) kk=0 do n=1,npe do j=1,ijn_s(n) kk=kk+1 ii=ltosi_s(kk) jj=ltosj_s(kk) work(kk)=a(ii,jj) end do end do iret=nf90_close(gfile_loc) deallocate (uu,a) endif !mype call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,& work_sub,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) deallocate (work) return end subroutine gsi_fv3ncdf2d_read_v1 subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_read ! prgmmr: wu org: np22 date: 2017-10-10 ! lei re-write for parallelization date: 2021-09-29 ! similar for horizontal recurisve filtering ! abstract: read in fields excluding u and v ! program history log: ! ! input argument list: ! filename - file name to read from ! varname - variable name to read in ! varname2 - variable name to read in ! mype_io - pe to read in the field ! ! output argument list: ! work_sub - output sub domain field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,mpi_rtype,mype use mpimod, only: MPI_INFO_NULL 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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none type(sub2grid_info), intent(in):: grd_ionouv type(gsi_bundle),intent(inout) :: cstate_nouv character(*),intent(in):: filenamein type (type_fv3regfilenameg),intent(in) ::fv3filenamegin 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 character(len=max_varname_length) :: name character(len=max_varname_length) :: filenamein2 integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) integer(i_kind) ilev,ilevtot,inative integer(i_kind) kbgn,kend integer(i_kind) gfile_loc,iret,var_id integer(i_kind) nz,nzp1,mm1 ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: uu2d_layout integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout mm1=mype+1 nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat nxcase=nx nycase=ny kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate(uu2d(nxcase,nycase)) 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(filenamein),'.',nio iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret call flush(6) call stop2(333) endif enddo else iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret call flush(6) call stop2(333) endif endif do ilevtot=kbgn,kend vgsiname=grd_ionouv%names(1,ilevtot) if(trim(vgsiname)=='delzinc') cycle !delzinc is not read from DZ ,it's started from hydrostatic height call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) name=trim(varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev startloc=(/1,1,inative/) countloc=(/nxcase,nycase,1/) 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))) 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 deallocate(uu2d_layout) enddo else iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) endif call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) enddo ! ilevtot if(fv3_io_layout_y > 1) then do nio=1,fv3_io_layout_y-1 iret=nf90_close(gfile_loc_layout(nio)) enddo deallocate(gfile_loc_layout) else iret=nf90_close(gfile_loc) endif deallocate (uu2d) call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) return end subroutine gsi_fv3ncdf_read subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_read _v1 ! Lei modified from gsi_fv3ncdf_read ! prgmmr: wu org: np22 date: 2017-10-10 ! ! abstract: read in a field from a netcdf FV3 file in mype_io ! then scatter the field to each PE ! program history log: ! ! input argument list: ! filename - file name to read from ! varname - variable name to read in ! varname2 - variable name to read in ! mype_io - pe to read in the field ! ! output argument list: ! work_sub - output sub domain field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: mpi_rtype,mpi_comm_world,mype,MPI_INFO_NULL use mpimod, only: mpi_comm_world,mpi_rtype,mype 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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll 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 type (type_fv3regfilenameg) :: fv3filenamegin type(gsi_bundle),intent(inout) :: cstate_nouv 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) :: filenamein2 character(len=max_varname_length) :: varname,vgsiname integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) integer(i_kind) kbgn,kend integer(i_kind) var_id integer(i_kind) inative,ilev,ilevtot integer(i_kind) gfile_loc,iret integer(i_kind) nzp1,mm1 mm1=mype+1 nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat nxcase=nx nycase=ny kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate(uu2d(nxcase,nycase)) iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret call flush(6) call stop2(333) endif do ilevtot=kbgn,kend vgsiname=grd_ionouv%names(1,ilevtot) call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev startloc=(/1,1,inative+1/) countloc=(/nxcase,nycase,1/) iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) if(iret/=nf90_noerr) then write(6,*)' wrong to get var_id ',var_id call stop2(333) endif 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) enddo ! i call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) iret=nf90_close(gfile_loc) deallocate (uu2d) return end subroutine gsi_fv3ncdf_read_v1 subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_readuv ! prgmmr: wu w org: np22 date: 2017-11-22 ! ! abstract: read in a field from a netcdf FV3 file in mype_u,mype_v ! then scatter the field to each PE ! program history log: ! ! input argument list: ! ! output argument list: ! ges_u - output sub domain u field ! ges_v - output sub domain v field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,mpi_rtype,mype,mpi_info_null 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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none type(sub2grid_info), intent(in):: grd_uv 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 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 real(r_kind),allocatable,dimension(:,:):: uc2d,vc2d character(len=max_varname_length) :: filenamein2 character(len=max_varname_length) :: varname,vgsiname real(r_kind),allocatable,dimension(:,:,:,:):: worksub integer(i_kind) u_grd_VarId,v_grd_VarId integer(i_kind) nlatcase,nloncase integer(i_kind) nxcase,nycase integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) integer(i_kind) inative,ilev,ilevtot integer(i_kind) kbgn,kend integer(i_kind) gfile_loc,iret integer(i_kind) nz,nzp1,mm1 ! for fv3_io_layout_y > 1 real(r_kind),allocatable,dimension(:,:):: u2d_layout,v2d_layout integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat nxcase=nx nycase=ny kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate(u2d(nxcase,nycase+1)) allocate(v2d(nxcase+1,nycase)) allocate(uc2d(nxcase,nycase)) allocate(vc2d(nxcase,nycase)) allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) filenamein=fv3filenamegin%dynvars 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(filenamein),".",nio iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) if(iret/=nf90_noerr) then write(6,*)'problem opening6 ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret call flush(6) call stop2(333) endif enddo else iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' problem opening6 ',trim(filenamein),', Status = ',iret call flush(6) call stop2(333) endif endif do ilevtot=kbgn,kend vgsiname=grd_uv%names(1,ilevtot) call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_uv%lnames(1,ilevtot) nz=grd_uv%nsig nzp1=nz+1 inative=nzp1-ilev u_countloc=(/nxcase,nycase+1,1/) v_countloc=(/nxcase+1,nycase,1/) u_startloc=(/1,1,inative/) 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)) 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))) 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 deallocate(v2d_layout) enddo else call check( nf90_inq_varid(gfile_loc,'u',u_grd_VarId) ) iret=nf90_get_var(gfile_loc,u_grd_VarId,u2d,start=u_startloc,count=u_countloc) call check( nf90_inq_varid(gfile_loc,'v',v_grd_VarId) ) iret=nf90_get_var(gfile_loc,v_grd_VarId,v2d,start=v_startloc,count=v_countloc) endif if(.not.grid_reverse_flag) then 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) ! NOTE on transfor to earth u/v: ! The u and v before transferring need to be in E-W/N-S grid, which is ! defined as reversed grid here because it is revered from map view. ! ! Have set the following flag for grid orientation ! grid_reverse_flag=true: E-W/N-S grid ! grid_reverse_flag=false: W-E/S-N grid ! ! So for preparing the wind transferring, need to reverse the grid from ! W-E/S-N grid to E-W/N-S grid when grid_reverse_flag=false: ! ! if(.not.grid_reverse_flag) call reverse_grid_r_uv ! ! 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.) enddo ! i if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 iret=nf90_close(gfile_loc_layout(nio)) enddo deallocate(gfile_loc_layout) else iret=nf90_close(gfile_loc) endif deallocate(u2d,v2d,uc2d,vc2d) call general_grid2sub(grd_uv,hwork,worksub) ges_u=worksub(1,:,:,:) ges_v=worksub(2,:,:,:) deallocate(worksub) end subroutine gsi_fv3ncdf_readuv subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) !$$$ subprogram documentation block ! subprogram: gsi_fv3ncdf_readuv_v1 ! prgmmr: wu w org: np22 date: 2017-11-22 ! program history log: ! 2019-04 lei modified from gsi_fv3ncdf_readuv to deal with cold start files . ! abstract: read in a field from a "cold start" netcdf FV3 file in mype_u,mype_v ! then scatter the field to each PE ! program history log: ! ! input argument list: ! ! output argument list: ! ges_u - output sub domain u field ! ges_v - output sub domain v field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use constants, only: half use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,mpi_rtype,mype,mpi_info_null 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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth use general_sub2grid_mod, only: sub2grid_info,general_grid2sub implicit none type(sub2grid_info), intent(in):: grd_uv 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 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 real(r_kind),allocatable,dimension(:,:):: uorv2d real(r_kind),allocatable,dimension(:,:,:,:):: worksub character(len=max_varname_length) :: filenamein2 character(len=max_varname_length) :: varname integer(i_kind) nlatcase,nloncase integer(i_kind) kbgn,kend integer(i_kind) var_id integer(i_kind) gfile_loc,iret integer(i_kind) j,nzp1,mm1 integer(i_kind) ilev,ilevtot,inative integer(i_kind) nxcase,nycase integer(i_kind) us_countloc(3),us_startloc(3) integer(i_kind) vw_countloc(3),vw_startloc(3) allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat nxcase=nx nycase=ny kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate (us2d(nxcase,nycase+1),vw2d(nxcase+1,nycase)) allocate (uorv2d(nxcase,nycase)) filenamein=fv3filenamegin%dynvars iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret call flush(6) call stop2(333) endif do ilevtot=kbgn,kend varname=grd_uv%names(1,ilevtot) filenamein2=fv3filenamegin%dynvars if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_uv%lnames(1,ilevtot) 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/) us_startloc=(/1,1,inative+1/) vw_startloc=(/1,1,inative+1/) ! transfor to earth u/v, interpolate to analysis grid, reverse vertical order iret=nf90_inq_varid(gfile_loc,trim(adjustl("u_s")),var_id) iret=nf90_get_var(gfile_loc,var_id,us2d,start=us_startloc,count=us_countloc) iret=nf90_inq_varid(gfile_loc,trim(adjustl("v_w")),var_id) iret=nf90_get_var(gfile_loc,var_id,vw2d,start=vw_startloc,count=vw_countloc) do j=1,ny 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) 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) enddo ! iilevtoto call general_grid2sub(grd_uv,hwork,worksub) ges_u=worksub(1,:,:,:) ges_v=worksub(2,:,:,:) iret=nf90_close(gfile_loc) deallocate (us2d,vw2d,worksub) end subroutine gsi_fv3ncdf_readuv_v1 subroutine wrfv3_netcdf(fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: wrfv3_netcdf write out FV3 analysis increments ! prgmmr: wu org: np22 date: 2017-10-23 ! ! abstract: write FV3 analysis in netcdf format ! ! program history log: ! 2019-04-18 CAPS(C. Tong) - import direct reflectivity DA capabilities ! 2019-11-22 CAPS(C. Tong) - modify "add_saved" to properly output analyses ! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use guess_grids, only: ntguessig,ges_tsen use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundleputvar use mpeu_util, only: die use gridmod, only: lat2,lon2,nsig use gridmod,only: l_reg_update_hydro_delz use guess_grids, only:geom_hgti,geom_hgti_bg use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt use directDA_radaruse_mod, only: l_use_dbz_directDA use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval use gridmod, only: eta1_ll,eta2_ll implicit none type (type_fv3regfilenameg),intent(in) :: fv3filenamegin ! Declare local constants logical add_saved ! variables for cloud info integer(i_kind) ier,istatus,it real(r_kind),pointer,dimension(:,: ):: ges_ps =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_u =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_v =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL() integer(i_kind) i,k real(r_kind),pointer,dimension(:,:,:):: ges_ql =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qi =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qr =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qs =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qg =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qnr =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_w =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delzinc =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,: ),allocatable:: ges_ps_write real(r_kind), dimension(lat2,lon2,nsig) :: io_arr_qr, io_arr_qs real(r_kind), dimension(lat2,lon2,nsig) :: io_arr_qg, io_arr_qnr real(r_kind), dimension(:,:,:),allocatable ::g_prsi it=ntguessig ier=0 call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ps' ,ges_ps ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus);ier=ier+istatus if (l_use_dbz_directDA) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus);ier=ier+istatus end if if(l_reg_update_hydro_delz) then allocate(ges_delzinc(lat2,lon2,nsig)) do k=1,nsig ges_delzinc(:,:,k)=geom_hgti(:,:,k+1,it)-geom_hgti_bg(:,:,k+1,it)-geom_hgti(:,:,k,it)+geom_hgti_bg(:,:,k,it) enddo call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'delzinc',ges_delzinc,istatus) endif ! additional I/O for direct reflectivity DA capabilities if (l_use_dbz_directDA) then if( fv3sar_bg_opt == 0) then write(6,*) "FV3 IO Write for 'fv3sar_bg_opt == 0 is only available for now in direct relfectivity DA" endif io_arr_qr=ges_qr io_arr_qs=ges_qs io_arr_qg=ges_qg io_arr_qnr=ges_qnr call convert_cvpqx_to_qx(io_arr_qr, io_arr_qs, io_arr_qg, l_use_cvpqx, cvpqx_pval) ! Convert Qx back call convert_cvpnx_to_nx(io_arr_qnr, l_cvpnr, cvpnr_pval, cld_nt_updt, ges_q, io_arr_qr, ges_ps) ! Convert Nx back ges_qr=io_arr_qr ges_qs=io_arr_qs ges_qg=io_arr_qg ges_qnr=io_arr_qnr endif call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_dynvar_nouv) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_tracer_nouv) call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'tsen',ges_tsen(:,:,:,it),istatus) if( fv3sar_bg_opt == 0) then call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'delp' ,ges_delp ,istatus );ier=ier+istatus allocate(g_prsi(lat2,lon2,nsig+1)) allocate(ges_ps_inc(lat2,lon2)) ges_ps_inc=ges_ps-ges_ps_bg g_prsi(:,:,nsig+1)=eta1_ll(nsig+1) do i=nsig,1,-1 g_prsi(:,:,i)=ges_delp_bg(:,:,i)+g_prsi(:,:,i+1) enddo do i=1,nsig+1 g_prsi(:,:,i)=g_prsi(:,:,i)+eta2_ll(i)*ges_ps_inc enddo do i=1,nsig ges_delp(:,:,i)=g_prsi(:,:,i)-g_prsi(:,:,i+1) enddo ges_delp=ges_delp*1000.0_r_kind deallocate(g_prsi,ges_ps_inc) else allocate(ges_ps_write(lat2,lon2)) ges_ps_write=ges_ps*1000.0_r_kind call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'ps',ges_ps_write,istatus) deallocate(ges_ps_write) endif ! write out if (ier/=0) call die('get ges','cannot get pointers for fv3 met-fields, ier =',ier) add_saved=.true. ! write out if(fv3sar_bg_opt == 0) then call gsi_fv3ncdf_write(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv,& add_saved,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_write(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv, & add_saved,fv3filenamegin%tracers,fv3filenamegin) call gsi_fv3ncdf_writeuv(grd_fv3lam_uv,ges_u,ges_v,add_saved,fv3filenamegin) else call gsi_fv3ncdf_write_v1(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv,& add_saved,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_write_v1(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv,& add_saved,fv3filenamegin%tracers,fv3filenamegin) call gsi_fv3ncdf_writeuv_v1(grd_fv3lam_uv,ges_u,ges_v,add_saved,fv3filenamegin) endif if(allocated(g_prsi)) deallocate(g_prsi) deallocate(ny_layout_len) deallocate(ny_layout_b) deallocate(ny_layout_e) ! additional I/O for direct reflectivity DA capabilities end subroutine wrfv3_netcdf subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeuv ! pgrmmr: wu ! ! abstract: gather u/v fields to mype_io, put u/v in FV3 model defined directions & orders ! then write out ! ! program history log: ! ! input argument list: ! varu,varv ! add_saved - true: add analysis increments to readin guess then write out ! - false: write out total analysis fields ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use gridmod, only: nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use netcdf, only: nf90_open,nf90_close,nf90_noerr use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none type(sub2grid_info), intent(in):: grd_uv real(r_kind),dimension(2,grd_uv%nlat,grd_uv%nlon,grd_uv%kbegin_loc:grd_uv%kend_alloc):: hwork logical ,intent(in ) :: add_saved type (type_fv3regfilenameg),intent(in) ::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 integer(i_kind) :: ugrd_VarId,gfile_loc,vgrd_VarId integer(i_kind) i,j,mm1,k,nzp1 integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) nlatcase,nloncase integer(i_kind) nxcase,nycase integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) character(:),allocatable:: filenamein ,varname real(r_kind),allocatable,dimension(:,:,:,:):: worksub real(r_kind),allocatable,dimension(:,:):: work_au,work_av real(r_kind),allocatable,dimension(:,:):: work_bu,work_bv real(r_kind),allocatable,dimension(:,:):: u2d,v2d,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu2,workbv2 ! for fv3_io_layout_y > 1 real(r_kind),allocatable,dimension(:,:):: u2d_layout,v2d_layout integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat nxcase=nx nycase=ny kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate( u2d(nlon_regional,nlat_regional+1)) allocate( v2d(nlon_regional+1,nlat_regional)) allocate( work_bu(nlon_regional,nlat_regional+1)) allocate( work_bv(nlon_regional+1,nlat_regional)) allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) allocate( work_au(nlatcase,nloncase),work_av(nlatcase,nloncase)) do k=1,grd_uv%nsig do j=1,grd_uv%lon2 do i=1,grd_uv%lat2 worksub(1,i,j,k)=ges_u(i,j,k) worksub(2,i,j,k)=ges_v(i,j,k) end do end do end do call general_sub2grid(grd_uv,worksub,hwork) filenamein=fv3filenamegin%dynvars 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(filenamein),".",nio call check( nf90_open(filename_layout,nf90_write,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) ) enddo gfile_loc=gfile_loc_layout(0) else call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) ) endif do ilevtot=kbgn,kend varname=grd_uv%names(1,ilevtot) ilev=grd_uv%lnames(1,ilevtot) nz=grd_uv%nsig nzp1=nz+1 inative=nzp1-ilev u_countloc=(/nxcase,nycase+1,1/) v_countloc=(/nxcase+1,nycase,1/) u_startloc=(/1,1,inative/) v_startloc=(/1,1,inative/) work_au=hwork(1,:,:,ilevtot) work_av=hwork(2,:,:,ilevtot) call check( nf90_inq_varid(gfile_loc,'u',ugrd_VarId) ) call check( nf90_inq_varid(gfile_loc,'v',vgrd_VarId) ) if(add_saved)then allocate( workau2(nlatcase,nloncase),workav2(nlatcase,nloncase)) allocate( workbu2(nlon_regional,nlat_regional+1)) allocate( workbv2(nlon_regional+1,nlat_regional)) !!!!!!!! readin work_b !!!!!!!!!!!!!!!! if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) call check( nf90_get_var(gfile_loc_layout(nio),ugrd_VarId,u2d_layout,start=u_startloc,count=u_countloc) ) work_bu(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) if(nio==fv3_io_layout_y-1) work_bu(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) deallocate(u2d_layout) allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) v_countloc=(/nxcase+1,ny_layout_len(nio),1/) call check( nf90_get_var(gfile_loc_layout(nio),vgrd_VarId,v2d_layout,start=v_startloc,count=v_countloc) ) work_bv(:,ny_layout_b(nio):ny_layout_e(nio))=v2d_layout deallocate(v2d_layout) enddo else call check( nf90_get_var(gfile_loc,ugrd_VarId,work_bu,start=u_startloc,count=u_countloc) ) call check( nf90_get_var(gfile_loc,vgrd_VarId,work_bv,start=v_startloc,count=v_countloc) ) endif if(.not.grid_reverse_flag) then call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,1) call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,1) endif call fv3uv2earth(work_bu,work_bv,nlon_regional,nlat_regional,u2d,v2d) call fv3_h_to_ll(u2d,workau2,nlon_regional,nlat_regional,nloncase,nlatcase,.true.) call fv3_h_to_ll(v2d,workav2,nlon_regional,nlat_regional,nloncase,nlatcase,.true.) !!!!!!!! find analysis_inc: work_a !!!!!!!!!!!!!!!! work_au(:,:)=work_au(:,:)-workau2(:,:) work_av(:,:)=work_av(:,:)-workav2(:,:) call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call earthuv2fv3(u2d,v2d,nlon_regional,nlat_regional,workbu2,workbv2) !!!!!!!! add analysis_inc to readin work_b !!!!!!!!!!!!!!!! work_bu(:,:)=work_bu(:,:)+workbu2(:,:) work_bv(:,:)=work_bv(:,:)+workbv2(:,:) deallocate(workau2,workbu2,workav2,workbv2) else call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call earthuv2fv3(u2d,v2d,nlon_regional,nlat_regional,work_bu(:,:),work_bv(:,:)) endif if(.not.grid_reverse_flag) then call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,1) call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,1) endif if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) u2d_layout=work_bu(:,ny_layout_b(nio):ny_layout_e(nio)+1) call check( nf90_put_var(gfile_loc_layout(nio),ugrd_VarId,u2d_layout,start=u_startloc,count=u_countloc) ) deallocate(u2d_layout) allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) v_countloc=(/nxcase+1,ny_layout_len(nio),1/) v2d_layout=work_bv(:,ny_layout_b(nio):ny_layout_e(nio)) call check( nf90_put_var(gfile_loc_layout(nio),vgrd_VarId,v2d_layout,start=v_startloc,count=v_countloc) ) deallocate(v2d_layout) enddo else call check( nf90_put_var(gfile_loc,ugrd_VarId,work_bu,start=u_startloc,count=u_countloc) ) call check( nf90_put_var(gfile_loc,vgrd_VarId,work_bv,start=v_startloc,count=v_countloc) ) endif enddo !ilevltot if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 call check( nf90_close(gfile_loc_layout(nio)) ) enddo deallocate(gfile_loc_layout) else call check( nf90_close(gfile_loc) ) endif deallocate(work_bu,work_bv,u2d,v2d) deallocate(work_au,work_av) end subroutine gsi_fv3ncdf_writeuv subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeuv ! pgrmmr: wu ! ! abstract: gather u/v fields to mype_io, put u/v in FV3 model defined directions & orders ! then write out ! ! program history log: ! 2019-04-22 lei modified from gsi_nemsio_writeuv_v1 for update ! u_w,v_w,u_s,v_s in the cold start files! ! 2020-03-06 lei added ilev0 fix ! input argument list: ! varu,varv ! add_saved - true: add analysis increments to readin guess then write out ! - false: write out total analysis fields ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use constants, only: half,zero use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null use gridmod, only: nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use netcdf, only: nf90_open,nf90_close,nf90_noerr use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none type(sub2grid_info), intent(in):: grd_uv 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 logical ,intent(in ) :: add_saved type (type_fv3regfilenameg),intent (in) :: fv3filenamegin real(r_kind),dimension(2,grd_uv%nlat,grd_uv%nlon,grd_uv%kbegin_loc:grd_uv%kend_alloc):: hwork character(len=:),allocatable :: filenamein character(len=max_varname_length) :: varname integer(i_kind) :: gfile_loc integer(i_kind) :: u_wgrd_VarId,v_wgrd_VarId integer(i_kind) :: u_sgrd_VarId,v_sgrd_VarId integer(i_kind) i,j,mm1,k,nzp1 integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot real(r_kind),allocatable,dimension(:,:,:,:):: worksub real(r_kind),allocatable,dimension(:,:):: work_au,work_av real(r_kind),allocatable,dimension(:,:):: work_bu_s,work_bv_s real(r_kind),allocatable,dimension(:,:):: work_bu_w,work_bv_w real(r_kind),allocatable,dimension(:,:):: u2d,v2d,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu_s2,workbv_s2 real(r_kind),allocatable,dimension(:,:):: workbu_w2,workbv_w2 integer(i_kind) nlatcase,nloncase,nxcase,nycase integer(i_kind) uw_countloc(3),us_countloc(3),uw_startloc(3),us_startloc(3) integer(i_kind) vw_countloc(3),vs_countloc(3),vw_startloc(3),vs_startloc(3) mm1=mype+1 nloncase=grd_uv%nlon nlatcase=grd_uv%nlat nxcase=nx nycase=ny kbgn=grd_uv%kbegin_loc kend=grd_uv%kend_loc allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) do k=1,grd_uv%nsig do j=1,grd_uv%lon2 do i=1,grd_uv%lat2 worksub(1,i,j,k)=ges_u(i,j,k) worksub(2,i,j,k)=ges_v(i,j,k) end do end do end do call general_sub2grid(grd_uv,worksub,hwork) allocate( u2d(nlon_regional,nlat_regional)) allocate( v2d(nlon_regional,nlat_regional)) allocate( work_bu_s(nlon_regional,nlat_regional+1)) allocate( work_bv_s(nlon_regional,nlat_regional+1)) allocate( work_bu_w(nlon_regional+1,nlat_regional)) allocate( work_bv_w(nlon_regional+1,nlat_regional)) allocate( work_au(nlatcase,nloncase),work_av(nlatcase,nloncase)) if(add_saved) allocate( workau2(nlatcase,nloncase),workav2(nlatcase,nloncase)) allocate( workbu_w2(nlon_regional+1,nlat_regional)) allocate( workbv_w2(nlon_regional+1,nlat_regional)) allocate( workbu_s2(nlon_regional,nlat_regional+1)) allocate( workbv_s2(nlon_regional,nlat_regional+1)) filenamein=fv3filenamegin%dynvars call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) ) do ilevtot=kbgn,kend varname=grd_uv%names(1,ilevtot) ilev=grd_uv%lnames(1,ilevtot) nz=grd_uv%nsig nzp1=nz+1 inative=nzp1-ilev uw_countloc= (/nlon_regional+1,nlat_regional,1/) us_countloc= (/nlon_regional,nlat_regional+1,1/) vw_countloc= (/nlon_regional+1,nlat_regional,1/) vs_countloc= (/nlon_regional,nlat_regional+1,1/) uw_startloc=(/1,1,inative+1/) us_startloc=(/1,1,inative+1/) vw_startloc=(/1,1,inative+1/) vs_startloc=(/1,1,inative+1/) work_au=hwork(1,:,:,ilevtot) work_av=hwork(2,:,:,ilevtot) call check( nf90_inq_varid(gfile_loc,'u_s',u_sgrd_VarId) ) call check( nf90_inq_varid(gfile_loc,'u_w',u_wgrd_VarId) ) call check( nf90_inq_varid(gfile_loc,'v_s',v_sgrd_VarId) ) call check( nf90_inq_varid(gfile_loc,'v_w',v_wgrd_VarId) ) !!!!!!!! readin work_b !!!!!!!!!!!!!!!! call check( nf90_get_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) call check( nf90_get_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) call check( nf90_get_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) call check( nf90_get_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) if(add_saved)then do j=1,nlat_regional u2d(:,j)=half * (work_bu_s(:,j)+ work_bu_s(:,j+1)) enddo do i=1,nlon_regional v2d(i,:)=half*(work_bv_w(i,:)+work_bv_w(i+1,:)) enddo call fv3_h_to_ll(u2d,workau2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) call fv3_h_to_ll(v2d,workav2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! find analysis_inc: work_a !!!!!!!!!!!!!!!! work_au(:,:)=work_au(:,:)-workau2(:,:) work_av(:,:)=work_av(:,:)-workav2(:,:) call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) !!!!!!!! add analysis_inc to readin work_b !!!!!!!!!!!!!!!! do i=2,nlon_regional workbu_w2(i,:)=half*(u2d(i-1,:)+u2d(i,:)) workbv_w2(i,:)=half*(v2d(i-1,:)+v2d(i,:)) enddo workbu_w2(1,:)=u2d(1,:) workbv_w2(1,:)=v2d(1,:) workbu_w2(nlon_regional+1,:)=u2d(nlon_regional,:) workbv_w2(nlon_regional+1,:)=v2d(nlon_regional,:) do j=2,nlat_regional workbu_s2(:,j)=half*(u2d(:,j-1)+u2d(:,j)) workbv_s2(:,j)=half*(v2d(:,j-1)+v2d(:,j)) enddo workbu_s2(:,1)=u2d(:,1) workbv_s2(:,1)=v2d(:,1) workbu_s2(:,nlat_regional+1)=u2d(:,nlat_regional) workbv_s2(:,nlat_regional+1)=v2d(:,nlat_regional) work_bu_w(:,:)=work_bu_w(:,:)+workbu_w2(:,:) work_bu_s(:,:)=work_bu_s(:,:)+workbu_s2(:,:) work_bv_w(:,:)=work_bv_w(:,:)+workbv_w2(:,:) work_bv_s(:,:)=work_bv_s(:,:)+workbv_s2(:,:) else call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) do i=2,nlon_regional work_bu_w(i,:)=half*(u2d(i-1,:)+u2d(i,:)) work_bv_w(i,:)=half*(v2d(i-1,:)+v2d(i,:)) enddo work_bu_w(1,:)=u2d(1,:) work_bv_w(1,:)=v2d(1,:) work_bu_w(nlon_regional+1,:)=u2d(nlon_regional,:) work_bv_w(nlon_regional+1,:)=v2d(nlon_regional,:) do j=2,nlat_regional work_bu_s(:,j)=half*(u2d(:,j-1)+u2d(:,j)) work_bv_s(:,j)=half*(v2d(:,j-1)+v2d(:,j)) enddo work_bu_s(:,1)=u2d(:,1) work_bv_s(:,1)=v2d(:,1) work_bu_s(:,nlat_regional+1)=u2d(:,nlat_regional) work_bv_s(:,nlat_regional+1)=v2d(:,nlat_regional) endif call check( nf90_put_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) call check( nf90_put_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) call check( nf90_put_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) call check( nf90_put_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) enddo ! call check( nf90_close(gfile_loc) ) deallocate(work_bu_w,work_bv_w) deallocate(work_bu_s,work_bv_s) deallocate(work_au,work_av,u2d,v2d) if(add_saved) deallocate(workau2,workav2) if (allocated(workbu_w2)) then deallocate(workbu_w2,workbv_w2) deallocate(workbu_s2,workbv_s2) endif if(allocated(worksub))deallocate(worksub) end subroutine gsi_fv3ncdf_writeuv_v1 subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write ! pgrmmr: wu ! ! abstract: ! ! program history log: ! ! input argument list: ! varu,varv ! add_saved ! mype - mpi task id ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null use mod_fv3_lola, only: fv3_ll_to_h use mod_fv3_lola, only: fv3_h_to_ll use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none type(sub2grid_info), intent(in):: grd_ionouv type(gsi_bundle),intent(inout) :: cstate_nouv logical ,intent(in ) :: add_saved character(len=:), allocatable, intent(in) :: filenamein type (type_fv3regfilenameg),intent (in) :: fv3filenamegin 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) :: filenamein2 character(len=max_varname_length) :: varname,vgsiname integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) :: VarId,gfile_loc integer(i_kind) mm1,nzp1 real(r_kind),allocatable,dimension(:,:):: work_a real(r_kind),allocatable,dimension(:,:):: work_b real(r_kind),allocatable,dimension(:,:):: workb2,worka2 ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: work_b_layout integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout mm1=mype+1 ! Convert from subdomain to full horizontal field distributed among ! processors call general_sub2grid(grd_ionouv,cstate_nouv%values,hwork) nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat nxcase=nx nycase=ny kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate( work_a(nlatcase,nloncase)) allocate( work_b(nlon_regional,nlat_regional)) allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlatcase,nloncase)) 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(filenamein),'.',nio call check( nf90_open(filename_layout,nf90_write,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) ) enddo gfile_loc=gfile_loc_layout(0) else call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) ) endif do ilevtot=kbgn,kend vgsiname=grd_ionouv%names(1,ilevtot) call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev countloc=(/nxcase,nycase,1/) startloc=(/1,1,inative/) work_a=hwork(1,:,:,ilevtot) call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) if(index(vgsiname,"delzinc") > 0) then if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 countloc=(/nxcase,ny_layout_len(nio),1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) call check( nf90_get_var(gfile_loc_layout(nio),VarId,work_b_layout,start = startloc, count = countloc) ) work_b(:,ny_layout_b(nio):ny_layout_e(nio))=work_b_layout deallocate(work_b_layout) enddo else call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) endif call fv3_ll_to_h(work_a(:,:),workb2,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) work_b(:,:)=work_b(:,:)+workb2(:,:) else if(add_saved)then if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 countloc=(/nxcase,ny_layout_len(nio),1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) call check( nf90_get_var(gfile_loc_layout(nio),VarId,work_b_layout,start = startloc, count = countloc) ) work_b(:,ny_layout_b(nio):ny_layout_e(nio))=work_b_layout deallocate(work_b_layout) enddo else call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) endif call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! work_a(:,:)=work_a(:,:)-worka2(:,:) call fv3_ll_to_h(work_a(:,:),workb2,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) work_b(:,:)=work_b(:,:)+workb2(:,:) else call fv3_ll_to_h(work_a(:,:),work_b(:,:),nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) endif endif if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 countloc=(/nxcase,ny_layout_len(nio),1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) work_b_layout=work_b(:,ny_layout_b(nio):ny_layout_e(nio)) call check( nf90_put_var(gfile_loc_layout(nio),VarId,work_b_layout, start = startloc, count = countloc) ) deallocate(work_b_layout) enddo else call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) endif enddo !ilevtotl loop if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 call check(nf90_close(gfile_loc_layout(nio))) enddo deallocate(gfile_loc_layout) else call check(nf90_close(gfile_loc)) endif deallocate(work_b,work_a) deallocate(workb2,worka2) end subroutine gsi_fv3ncdf_write subroutine check(status) use kinds, only: i_kind use netcdf, only: nf90_noerr,nf90_strerror integer(i_kind), intent ( in) :: status if(status /= nf90_noerr) then print *,'ncdf error ', trim(nf90_strerror(status)) stop end if end subroutine check subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write ! pgrmmr: wu ! ! abstract: ! ! program history log: ! 2020-03-05 lei modified from gsi_fv3ncdf_write to gsi_fv3ncdf_write_v1 ! input argument list: ! varu,varv ! add_saved ! mype - mpi task id ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null use mod_fv3_lola, only: fv3_ll_to_h use mod_fv3_lola, only: fv3_h_to_ll use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none type(sub2grid_info), intent(in):: grd_ionouv type(gsi_bundle),intent(inout) :: cstate_nouv logical ,intent(in ) :: add_saved character(*),intent(in):: filenamein type (type_fv3regfilenameg),intent (in) :: fv3filenamegin 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) :: filenamein2 integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) :: VarId,gfile_loc integer(i_kind) mm1,nzp1 real(r_kind),allocatable,dimension(:,:):: work_a real(r_kind),allocatable,dimension(:,:):: work_b real(r_kind),allocatable,dimension(:,:):: workb2,worka2 character(len=max_varname_length) :: varname,vgsiname integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) mm1=mype+1 nloncase=grd_ionouv%nlon nlatcase=grd_ionouv%nlat call general_sub2grid(grd_ionouv,cstate_nouv%values,hwork) nxcase=nx nycase=ny kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc allocate( work_a(nlatcase,nloncase)) allocate( work_b(nlon_regional,nlat_regional)) allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlatcase,nloncase)) call check ( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL)) !clt do ilevtot=kbgn,kend vgsiname=grd_ionouv%names(1,ilevtot) call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev startloc=(/1,1,inative+1/) countloc=(/nxcase,nycase,1/) work_a=hwork(1,:,:,ilevtot) call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) call check( nf90_get_var(gfile_loc,VarId,work_b,start=startloc,count=countloc) ) if(index(vgsiname,"delzinc") > 0) then write(6,*)'delz is not in the cold start fiels with this option, incompatible setup , stop' call stop2(333) endif if(add_saved)then ! for being now only lev between (including ) 2 and nsig+1 of work_b (:,:,lev) ! are updated call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! work_a(:,:)=work_a(:,:)-worka2(:,:) call fv3_ll_to_h(work_a(:,:),workb2,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) work_b(:,:)=work_b(:,:)+workb2(:,:) else call fv3_ll_to_h(work_a(:,:),work_b(:,:),nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) endif call check( nf90_put_var(gfile_loc,VarId,work_b,start=startloc,count=countloc) ) enddo !ilevtot call check(nf90_close(gfile_loc)) deallocate(work_b,work_a) deallocate(worka2,workb2) end subroutine gsi_fv3ncdf_write_v1 subroutine reverse_grid_r(grid,nx,ny,nz) ! ! reverse the first two dimension of the array grid ! use kinds, only: r_kind,i_kind implicit none integer(i_kind), intent(in ) :: nx,ny,nz real(r_kind), intent(inout ) :: grid(nx,ny,nz) real(r_kind) :: tmp_grid(nx,ny) integer(i_kind) :: i,j,k ! do k=1,nz tmp_grid(:,:)=grid(:,:,k) do j=1,ny do i=1,nx grid(i,j,k)=tmp_grid(nx+1-i,ny+1-j) enddo enddo enddo end subroutine reverse_grid_r subroutine reverse_grid_r_uv(grid,nx,ny,nz) ! ! reverse the first two dimension of the array grid ! use kinds, only: r_kind,i_kind implicit none integer(i_kind), intent(in ) :: nx,ny,nz real(r_kind), intent(inout ) :: grid(nx,ny,nz) real(r_kind) :: tmp_grid(nx,ny) integer(i_kind) :: i,j,k ! do k=1,nz tmp_grid(:,:)=grid(:,:,k) do j=1,ny do i=1,nx grid(i,j,k)=-tmp_grid(nx+1-i,ny+1-j) enddo enddo enddo end subroutine reverse_grid_r_uv subroutine convert_qx_to_cvpqx(qr_arr,qs_arr,qg_arr,use_cvpqx,cvpqx_pvalue) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_qx_to_cvpqx ! prgmmr: J. Park(CAPS) date: 2021-05-05 ! ! abstract: convert qx(mixing ratio) to cvpqx using power transform for qr, qs, qg ! ! program history log: ! 2021-05-05 - initial commit ! - this is used when GSI reads qx data from a background file ! (subroutine read_fv3_netcdf_guess) ! - since minimum qr, qs, and qg are set for CVlogq, ! it reads three qx arrays and then processes. ! ! input argument list: ! qr_arr - array of qr ! qs_arr - array of qs ! qg_arr - array of qg ! use_cvpqx - flag to use power transform or not ! cvpqx_pvalue - value to be used for power transform ! ! output argument list: ! qr_arr - updated array of qr after power transform ! qs_arr - updated array of qs after ! qg_arr - updated array of qg after power transfrom ! ! attributes: ! language: f90 ! use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig use guess_grids, only: ges_tsen use mpimod, only: mype use guess_grids, only: ntguessig use constants, only: zero, one_tenth implicit none real(r_kind), intent(inout ) :: qr_arr(lat2,lon2,nsig) real(r_kind), intent(inout ) :: qs_arr(lat2,lon2,nsig) real(r_kind), intent(inout ) :: qg_arr(lat2,lon2,nsig) logical, intent(in ) :: use_cvpqx real(r_kind), intent(in ) :: cvpqx_pvalue integer(i_kind) :: i, j, k, it real(r_kind) :: qr_min, qs_min, qg_min real(r_kind) :: qr_thrshd, qs_thrshd, qg_thrshd ! it=ntguessig ! ! print info message: CVq, CVlogq, and CVpq if(mype==0)then if (use_cvpqx) then if ( cvpqx_pvalue == 0._r_kind ) then ! CVlogq write(6,*)'read_fv3_netcdf_guess: ', & ' reset zero of qr/qs/qg to specified values (~0dbz)', & 'before log transformation. (for dbz assimilation)' write(6,*)'read_fv3_netcdf_guess: convert qr/qs/qg to log transform.' else if ( cvpqx_pvalue > 0._r_kind ) then ! CVpq write(6,*)'read_fv3_netcdf_guess: convert qr/qs/qg with power transform .' end if else ! CVq write(6,*)'read_fv3_netcdf_guess: only reset (qr/qs/qg) to & 0.0 for negative analysis value. (regular qx)' end if end if do k=1,nsig do i=1,lon2 do j=1,lat2 ! Apply power transform if option is ON if (use_cvpqx) then if ( cvpqx_pvalue == 0._r_Kind ) then ! CVlogq if (ges_tsen(j,i,k,it) > 274.15_r_kind) then qr_min=2.9E-6_r_kind qr_thrshd=qr_min * one_tenth qs_min=0.1E-9_r_kind qs_thrshd=qs_min qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth else if (ges_tsen(j,i,k,it) <= 274.15_r_kind .and. & ges_tsen(j,i,k,it) >= 272.15_r_kind) then qr_min=2.0E-6_r_kind qr_thrshd=qr_min * one_tenth qs_min=1.3E-7_r_kind qs_thrshd=qs_min * one_tenth qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth else if (ges_tsen(j,i,k,it) < 272.15_r_kind) then qr_min=0.1E-9_r_kind qr_thrshd=qr_min qs_min=6.3E-6_r_kind qs_thrshd=qs_min * one_tenth qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth end if if ( qr_arr(j,i,k) <= qr_thrshd ) qr_arr(j,i,k) = qr_min if ( qs_arr(j,i,k) <= qs_thrshd ) qs_arr(j,i,k) = qs_min if ( qg_arr(j,i,k) <= qg_thrshd ) qg_arr(j,i,k) = qg_min qr_arr(j,i,k) = log(qr_arr(j,i,k)) qs_arr(j,i,k) = log(qs_arr(j,i,k)) qg_arr(j,i,k) = log(qg_arr(j,i,k)) else if ( cvpqx_pvalue > 0._r_kind ) then ! CVpq qr_arr(j,i,k)=((max(qr_arr(j,i,k),1.0E-6_r_kind))**cvpqx_pvalue-1)/cvpqx_pvalue qs_arr(j,i,k)=((max(qs_arr(j,i,k),1.0E-6_r_kind))**cvpqx_pvalue-1)/cvpqx_pvalue qg_arr(j,i,k)=((max(qg_arr(j,i,k),1.0E-6_r_kind))**cvpqx_pvalue-1)/cvpqx_pvalue end if else ! CVq qr_min=zero qs_min=zero qg_min=zero qr_arr(j,i,k) = max(qr_arr(j,i,k), qr_min) qs_arr(j,i,k) = max(qs_arr(j,i,k), qs_min) qg_arr(j,i,k) = max(qg_arr(j,i,k), qg_min) end if end do end do end do end subroutine convert_qx_to_cvpqx subroutine convert_nx_to_cvpnx(qnx_arr,cvpnr,cvpnr_pvalue) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_nx_to_cvpnx ! prgmmr: J. Park(CAPS) date: 2021-05-05 ! ! abstract: convert nx (number concentration) to cvpnx using power transform ! ! program history log: ! 2021-05-05 - initial commit ! - this is used when GSI reads nx data from a background file ! (subroutine read_fv3_netcdf_guess) ! - this can be used for other nx variables ! ! input argument list: ! qnx_arr - array of qnx ! cvpnr - flag to use power transform or not ! cvpnr_pvalue - value to be used for power transform ! ! output argument list: ! qnx_arr - updated array of qnx after power transform ! ! attributes: ! language: f90 ! use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig use mpimod, only: mype use constants, only: zero, one_tenth implicit none real(r_kind), intent(inout ) :: qnx_arr(lat2,lon2,nsig) logical, intent(in ) :: cvpnr real(r_kind), intent(in ) :: cvpnr_pvalue integer(i_kind) :: i, j, k ! ! print info message: CVpnr if (mype==0 .and. cvpnr)then write(6,*)'read_fv3_netcdf_guess: convert qnx with power transform .' end if do k=1,nsig do i=1,lon2 do j=1,lat2 ! Treatment on qnx ; power transform if (cvpnr) then qnx_arr(j,i,k)=((max(qnx_arr(j,i,k),1.0E-2_r_kind)**cvpnr_pvalue)-1)/cvpnr_pvalue endif end do end do end do end subroutine convert_nx_to_cvpnx subroutine convert_cvpqx_to_qx(qr_arr,qs_arr,qg_arr,use_cvpqx,cvpqx_pvalue) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_cvpqx_to_qx ! prgmmr: J. Park(CAPS) date: 2021-05-05 ! ! abstract: convert cvpqx to qx for qr, qs, qg ! ! program history log: ! 2021-05-05 - initial commit ! - this is used when GSI writes qx data to a background file ! (subroutine wrfv3_netcdf) ! - since minimum qr, qs, and qg are set for CVlogq, ! it reads three qx arrays and then processes. ! ! input argument list: ! qr_arr - array of qr ! qs_arr - array of qs ! qg_arr - array of qg ! use_cvpqx - flag to use power transform or not ! cvpqx_pvalue - value to be used for power transform ! ! output argument list: ! qr_arr - updated array of qr after power transform ! qs_arr - updated array of qs after ! qg_arr - updated array of qg after power transfrom ! ! attributes: ! language: f90 ! use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig use guess_grids, only: ges_tsen use mpimod, only: mype use guess_grids, only: ntguessig use constants, only: zero, one_tenth,r0_01 implicit none real(r_kind), intent(inout ) :: qr_arr(lat2,lon2,nsig) real(r_kind), intent(inout ) :: qs_arr(lat2,lon2,nsig) real(r_kind), intent(inout ) :: qg_arr(lat2,lon2,nsig) logical, intent(in ) :: use_cvpqx real(r_kind), intent(in ) :: cvpqx_pvalue integer(i_kind) :: i, j, k, it real(r_kind), dimension(lat2,lon2,nsig) :: tmparr_qr, tmparr_qs real(r_kind), dimension(lat2,lon2,nsig) :: tmparr_qg real(r_kind) :: qr_min, qs_min, qg_min real(r_kind) :: qr_tmp, qs_tmp, qg_tmp real(r_kind) :: qr_thrshd, qs_thrshd, qg_thrshd ! it=ntguessig ! ! print info message: CVq, CVlogq, and CVpq if(mype==0)then if (use_cvpqx) then if ( cvpqx_pvalue == 0._r_kind ) then ! CVlogq write(6,*)'wrfv3_netcdf: convert log(qr/qs/qg) back to qr/qs/qg.' write(6,*)'wrfv3_netcdf: then reset (qr/qs/qg) to 0.0 for some cases.' else if ( cvpqx_pvalue > 0._r_kind ) then ! CVpq write(6,*)'wrfv3_netcdf: convert power transformed (qr/qs/qg) back to qr/qs/qg.' write(6,*)'wrfv3_netcdf: then reset (qr/qs/qg) to 0.0 for some cases.' end if else ! CVq write(6,*)'wrfv3_netcdf: only reset (qr/qs/qg) to 0.0 for negative analysis value. (regular qx)' end if end if ! Initialized temporary arrays with ges. Will be recalculated later if cvlogq or cvpq is used tmparr_qr =qr_arr tmparr_qs =qs_arr tmparr_qg =qg_arr do k=1,nsig do i=1,lon2 do j=1,lat2 ! initialize hydrometeors as zero qr_tmp=zero qs_tmp=zero qg_tmp=zero if ( use_cvpqx ) then if ( cvpqx_pvalue == 0._r_kind ) then ! CVlogq if (ges_tsen(j,i,k,it) > 274.15_r_kind) then qr_min=2.9E-6_r_kind qr_thrshd=qr_min * one_tenth qs_min=0.1E-9_r_kind qs_thrshd=qs_min qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth else if (ges_tsen(j,i,k,it) <= 274.15_r_kind .and. & ges_tsen(j,i,k,it) >= 272.15_r_kind ) then qr_min=2.0E-6_r_kind qr_thrshd=qr_min * one_tenth qs_min=1.3E-7_r_kind qs_thrshd=qs_min * one_tenth qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth else if (ges_tsen(j,i,k,it) < 272.15_r_kind) then qr_min=0.1E-9_r_kind qr_thrshd=qr_min qs_min=6.3E-6_r_kind qs_thrshd=qs_min * one_tenth qg_min=3.1E-7_r_kind qg_thrshd=qg_min * one_tenth end if qr_tmp=exp(qr_arr(j,i,k)) qs_tmp=exp(qs_arr(j,i,k)) qg_tmp=exp(qg_arr(j,i,k)) ! if no update or very tiny value of qr/qs/qg, re-set/clear it ! off to zero if ( abs(qr_tmp - qr_min) < (qr_min*r0_01) ) then qr_tmp=zero else if (qr_tmp < qr_thrshd) then qr_tmp=zero end if if ( abs(qs_tmp - qs_min) < (qs_min*r0_01) ) then qs_tmp=zero else if (qs_tmp < qs_thrshd) then qs_tmp=zero end if if ( abs(qg_tmp - qg_min) < (qg_min*r0_01) ) then qg_tmp=zero else if (qg_tmp < qg_thrshd) then qg_tmp=zero end if else if ( cvpqx_pvalue > 0._r_kind ) then ! CVpq qr_tmp=max((cvpqx_pvalue*qr_arr(j,i,k)+1)**(1/cvpqx_pvalue)-1.0E-6_r_kind,0.0_r_kind) qs_tmp=max((cvpqx_pvalue*qs_arr(j,i,k)+1)**(1/cvpqx_pvalue)-1.0E-6_r_kind,0.0_r_kind) qg_tmp=max((cvpqx_pvalue*qg_arr(j,i,k)+1)**(1/cvpqx_pvalue)-1.0E-6_r_kind,0.0_r_kind) !Set a upper limit to hydrometeors to disable overshooting qr_tmp=min(qr_tmp,1E-2_r_kind) qs_tmp=min(qs_tmp,1.0E-2_r_kind) qg_tmp=min(qg_tmp,2E-2_r_kind) end if else ! For CVq qr_min=zero qs_min=zero qg_min=zero qr_tmp=qr_arr(j,i,k)-1.0E-8_r_kind qs_tmp=qs_arr(j,i,k)-1.0E-8_r_kind qg_tmp=qg_arr(j,i,k)-1.0E-8_r_kind qr_tmp=max(qr_tmp,qr_min) qs_tmp=max(qs_tmp,qs_min) qg_tmp=max(qg_tmp,qg_min) end if ! cvpqx tmparr_qr(j,i,k)=qr_tmp tmparr_qs(j,i,k)=qs_tmp tmparr_qg(j,i,k)=qg_tmp end do end do end do qr_arr=tmparr_qr qs_arr=tmparr_qs qg_arr=tmparr_qg end subroutine convert_cvpqx_to_qx subroutine convert_cvpnx_to_nx(qnx_arr,cvpnr,cvpnr_pvalue,cloud_nt_updt,q_arr,qr_arr,ps_arr) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_cvpnx_to_nx ! prgmmr: J. Park(CAPS) date: 2021-05-05 ! ! abstract: convert cvpnx to nx (number concentration) ! ! program history log: ! 2021-05-05 - initial commit ! - this is used when GSI writes nx data from a background file ! (subroutine wrfv3_netcdf) ! - this can be used for other nx variables ! ! input argument list: ! qnx_arr - array of qnx ! cvpnr - flag to use power transform or not ! cvpnr_pvalue - value to be used for power transform ! cloud_nt_updt - integer flag to use re-initialisation of QNRAIN with analyzed qr and n0r ! q_arr - array of qv, used only if cloud_nt_up_dt is 2 ! qr_arr - array of qr, used only if cloud_nt_up_dt is 2 ! ps_arr - array of ps, used only if cloud_nt_up_dt is 2 ! ! output argument list: ! qnx_arr - updated array of qnx after power transform ! ! attributes: ! language: f90 ! use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig use mpimod, only: mype use constants, only: zero, one, one_tenth use directDA_radaruse_mod, only: init_mm_qnr use guess_grids, only: ges_tsen use guess_grids, only: ntguessig use gridmod, only: pt_ll, aeta1_ll use constants, only: r10, r100, rd implicit none real(r_kind), intent(inout ) :: qnx_arr(lat2,lon2,nsig) logical, intent(in ) :: cvpnr real(r_kind), intent(in ) :: cvpnr_pvalue integer(i_kind), intent(in ) :: cloud_nt_updt real(r_kind), intent(in ) :: q_arr(lat2,lon2,nsig) real(r_kind), intent(in ) :: qr_arr(lat2,lon2,nsig) real(r_kind), intent(in ) :: ps_arr(lat2,lon2) real(r_kind), dimension(lat2,lon2,nsig) :: tmparr_qnr integer(i_kind) :: i, j, k, it real(r_kind) :: qnr_tmp real(r_kind) :: P1D,T1D,Q1D,RHO,QR1D real(r_kind),parameter:: D608=0.608_r_kind ! it=ntguessig ! ! print info message: CVpnr if (mype==0 .and. cvpnr)then write(6,*)'wrfv3_netcdf: convert power transformed (qnx) back to qnx.' end if ! Initialized temp arrays with ges. tmparr_qnr=qnx_arr do k=1,nsig do i=1,lon2 do j=1,lat2 ! re-initialisation of QNRAIN with analyzed qr and N0r(which is single-moment parameter) ! equation is used in subroutine init_MM of initlib3d.f90 in arps package qnr_tmp = zero if ( cloud_nt_updt == 2 ) then T1D=ges_tsen(j,i,k,it) ! sensible temperature (K) P1D=r100*(aeta1_ll(k)*(r10*ps_arr(j,i)-pt_ll)+pt_ll) ! pressure hPa --> Pa Q1D=q_arr(j,i,k)/(one-q_arr(j,i,k)) ! mixing ratio RHO=P1D/(rd*T1D*(one+D608*Q1D)) ! air density in kg m^-3 QR1D=qr_arr(j,i,k) CALL init_mm_qnr(RHO,QR1D,qnr_tmp) qnr_tmp = max(qnr_tmp, zero) else if (cvpnr) then ! power transform qnr_tmp=max((cvpnr_pvalue*qnx_arr(j,i,k)+1)**(1/cvpnr_pvalue)-1.0E-2_r_kind,0.0_r_kind) else qnr_tmp=qnx_arr(j,i,k) end if end if tmparr_qnr(j,i,k)=qnr_tmp end do end do end do qnx_arr=tmparr_qnr end subroutine convert_cvpnx_to_nx subroutine gsi_copy_bundle(bundi,bundo) use gsi_bundlemod, only:gsi_bundleinquire, gsi_bundlegetpointer,gsi_bundleputvar implicit none ! copy the variables in the gsi_metguess_bundle_inout to gsi_bundle_inout or ! vice versa, according to icopy_flag ! !INPUT PARAMETERS: type(gsi_bundle), intent(in ) :: bundi ! !INPUT/OUTPUT PARAMETERS: type(gsi_bundle), intent(inout) :: bundo character(len=max_varname_length),dimension(:),allocatable:: src_name_vars2d character(len=max_varname_length),dimension(:),allocatable:: src_name_vars3d character(len=max_varname_length),dimension(:),allocatable:: target_name_vars2d character(len=max_varname_length),dimension(:),allocatable:: target_name_vars3d character(len=max_varname_length) ::varname real(r_kind),dimension(:,:,:),pointer:: pvar3d=>NULL() real(r_kind),dimension(:,:,:),pointer:: pvar2d =>NULL() integer(i_kind):: src_nc3d,src_nc2d,target_nc3d,target_nc2d integer(i_kind):: ivar,jvar,istatus src_nc3d=bundi%n3d src_nc2d=bundi%n2d target_nc3d=bundo%n3d target_nc2d=bundo%n2d allocate(src_name_vars3d(src_nc3d),src_name_vars2d(src_nc2d)) allocate(target_name_vars3d(target_nc3d),target_name_vars2d(target_nc2d)) call gsi_bundleinquire(bundi,'shortnames::3d',src_name_vars3d,istatus) call gsi_bundleinquire(bundi,'shortnames::2d',src_name_vars2d,istatus) call gsi_bundleinquire(bundo,'shortnames::3d',target_name_vars3d,istatus) call gsi_bundleinquire(bundo,'shortnames::2d',target_name_vars2d,istatus) do ivar=1,src_nc3d varname=trim(src_name_vars3d(ivar)) do jvar=1,target_nc3d if(index(target_name_vars3d(jvar),varname) > 0) then call GSI_BundleGetPointer (bundi,varname,pvar3d,istatus) call gsi_bundleputvar (bundo,varname,pvar3d,istatus) exit endif enddo enddo do ivar=1,src_nc2d varname=trim(src_name_vars2d(ivar)) do jvar=1,target_nc2d if(index(target_name_vars2d(jvar),varname) > 0) then call GSI_BundleGetPointer (bundi,varname,pvar2d,istatus) call gsi_bundleputvar (bundo,varname,pvar2d,istatus) exit endif enddo enddo deallocate(src_name_vars3d,src_name_vars2d) deallocate(target_name_vars3d,target_name_vars2d) return end subroutine gsi_copy_bundle subroutine getfv3lamfilevname(vgsinamein,fv3filenamegref,filenameout,vname) type (type_fv3regfilenameg),intent (in) :: fv3filenamegref character(len=*):: vgsinamein character(len=*),intent(out):: vname character(len=*),intent(out):: filenameout if (ifindstrloc(vgsiname,vgsinamein)<= 0) then write(6,*)'the name ',vgsinamein ,'cannot be treated correctly in getfv3lamfilevname,stop' call stop2(333) endif if(ifindstrloc(vardynvars,vgsinamein)> 0) then filenameout=fv3filenamegref%dynvars else if(ifindstrloc(vartracers,vgsinamein)> 0 ) then filenameout=fv3filenamegref%tracers else write(6,*)'the filename corresponding to var ',trim(vgsinamein),' is not found, stop ' call stop2(333) endif vname=varfv3name(ifindstrloc(vgsiname,vgsinamein)) if(trim(vname)=="T".and. fv3sar_bg_opt==1) then vname="t" endif return end subroutine getfv3lamfilevname function ifindstrloc(str_array,strin) integer(i_kind) ifindstrloc character(len=max_varname_length),dimension(:) :: str_array character(len=*) :: strin integer(i_kind) i ifindstrloc=0 do i=1,size(str_array) if(trim(str_array(i)) == trim(strin)) then ifindstrloc=i exit endif enddo end function ifindstrloc end module gsi_rfv3io_mod