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 ! 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 check ! ! variable definitions: ! ! attributes: ! langauge: f90 ! machine: ! !$$$ end documentation block use kinds, only: r_kind,i_kind use gridmod, only: nlon_regional,nlat_regional 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 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(:) ! 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 :: 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 integer(i_kind) mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql 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 ) 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,aeta1_ll,aeta2_ll use gridmod, only:nlon_regional,nlat_regional,eta1_ll,eta2_ll 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=128) :: name integer(i_kind) myear,mmonth,mday,mhour,mminute,msecond real(r_kind),allocatable:: abk_fv3(:) 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 !!!!!!!!!! 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 if(mype==0)write(6,*),'nx,ny=',nx,ny !!! 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 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 guess_grids, only: ges_tsen,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_bundlegetpointer use mpeu_util, only: die use guess_grids, only: ntguessig 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_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_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_oz=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_tv=>NULL() character(len=:),allocatable :: dynvars !='fv3_dynvars' character(len=:),allocatable :: tracers !='fv3_tracer' dynvars= fv3filenamegin%dynvars tracers= fv3filenamegin%tracers if(npe< 8) then call die('read_fv3_netcdf_guess','not enough PEs to read in fv3 fields' ) endif mype_u=0 mype_v=1 mype_t=2 mype_p=3 mype_q=4 mype_ql=5 mype_oz=6 mype_2d=7 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 ! do it=1,nfldsig 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), '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), 'ql' ,ges_ql ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'oz' ,ges_oz ,istatus );ier=ier+istatus 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(dynvars,ges_u,ges_v) else call gsi_fv3ncdf_readuv_v1(dynvars,ges_u,ges_v) endif if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_read(dynvars,'T','t',ges_tsen(1,1,1,it),mype_t) else call gsi_fv3ncdf_read_v1(dynvars,'t','T',ges_tsen(1,1,1,it),mype_t) endif if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_read(dynvars,'DELP','delp',ges_prsi,mype_p) ges_prsi(:,:,nsig+1,it)=eta1_ll(nsig+1) do i=nsig,1,-1 ges_prsi(:,:,i,it)=ges_prsi(:,:,i,it)*0.001_r_kind+ges_prsi(:,:,i+1,it) enddo ges_ps(:,:)=ges_prsi(:,:,1,it) else call gsi_fv3ncdf2d_read_v1(dynvars,'ps','PS',ges_ps,mype_p) 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 if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_read(tracers,'SPHUM','sphum',ges_q,mype_q) ! call gsi_fv3ncdf_read(tracers,'LIQ_WAT','liq_wat',ges_ql,mype_ql) call gsi_fv3ncdf_read(tracers,'O3MR','o3mr',ges_oz,mype_oz) else call gsi_fv3ncdf_read_v1(tracers,'sphum','SPHUM',ges_q,mype_q) call gsi_fv3ncdf_read_v1(tracers,'o3mr','O3MR',ges_oz,mype_oz) endif !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig do j=1,lon2 do i=1,lat2 ges_tv(i,j,k)=ges_tsen(i,j,k,it)*(one+fv*ges_q(i,j,k)) enddo enddo enddo call gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) 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=128) :: 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' 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 iret=nf90_open(sfcdata,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening3 ',trim(sfcdata),', Status = ',iret return 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(allocated(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) iret=nf90_inquire_variable(gfile_loc,i,dimids=dim_id) if(allocated(sfc )) deallocate(sfc ) allocate(sfc(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) iret=nf90_get_var(gfile_loc,i,sfc) call fv3_h_to_ll(sfc(:,:,1),a,nx,ny,nxa,nya) 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 iret=nf90_close(gfile_loc) !!!! read in orog from dynam !!!!!!!!!!!! iret=nf90_open(trim(dynvars ),nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening4 ',trim(dynvars ),gfile_loc,', Status = ',iret return 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)) 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) exit endif enddo ! k iret=nf90_close(gfile_loc) k=k_orog call fv3_h_to_ll(sfc1,a,nx,ny,nxa,nya) 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) 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 integer(i_kind),allocatable,dimension(:):: dim_id,dim 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(trim(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(dim(ndimensions)) 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(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) iret=nf90_inquire_variable(gfile_loc,var_id,dimids=dim_id) 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) 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,dim,dim_id) 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(filenamein,varname,varname2,work_sub,mype_io) !$$$ subprogram documentation block ! . . . . ! subprogram: 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: ierror,mpi_comm_world,npe,mpi_rtype,mype use gridmod, only: lat2,lon2,nsig,nlat,nlon,itotsub,ijn_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 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,nsig) integer(i_kind) ,intent(in ) :: mype_io character(len=128) :: name real(r_kind),allocatable,dimension(:,:,:):: uu integer(i_kind),allocatable,dimension(:):: dim_id,dim real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: a integer(i_kind) n,ns,k,len,ndim integer(i_kind) gfile_loc,iret integer(i_kind) nz,nzp1,kk,j,mm1,i,ir,ii,jj integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid mm1=mype+1 allocate (work(itotsub*nsig)) if(mype==mype_io ) then iret=nf90_open(trim(filenamein),nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret write(6,*)' gsi_fv3ncdf_read:problem opening5 with varnam ',trim(varname) return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(dim(ndimensions)) allocate(a(nlat,nlon)) 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)==varname .or. trim(name)==varname2) then iret=nf90_inquire_variable(gfile_loc,k,ndims=ndim) if(allocated(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) iret=nf90_inquire_variable(gfile_loc,k,dimids=dim_id) if(allocated(uu )) deallocate(uu ) allocate(uu(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) iret=nf90_get_var(gfile_loc,k,uu) exit endif enddo ! k nz=nsig nzp1=nz+1 do i=1,nz ir=nzp1-i call fv3_h_to_ll(uu(:,:,i),a,dim(dim_id(1)),dim(dim_id(2)),nlon,nlat) kk=0 do n=1,npe ns=displss(n)+(ir-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 iret=nf90_close(gfile_loc) deallocate (uu,a,dim,dim_id) endif !mype call mpi_scatterv(work,ijns,displss,mpi_rtype,& work_sub,ijns(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) deallocate (work) return end subroutine gsi_fv3ncdf_read subroutine gsi_fv3ncdf_read_v1(filenamein,varname,varname2,work_sub,mype_io) !$$$ 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: ierror,mpi_comm_world,npe,mpi_rtype,mype use gridmod, only: lat2,lon2,nsig,nlat,nlon,itotsub,ijn_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 netcdf, only: nf90_inq_varid 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,nsig) integer(i_kind) ,intent(in ) :: mype_io character(len=128) :: name real(r_kind),allocatable,dimension(:,:,:):: uu real(r_kind),allocatable,dimension(:,:,:):: temp0 integer(i_kind),allocatable,dimension(:):: dim real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: a integer(i_kind) n,ns,k,len,var_id integer(i_kind) gfile_loc,iret integer(i_kind) nztmp,nzp1,kk,j,mm1,i,ir,ii,jj integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid mm1=mype+1 allocate (work(itotsub*nsig)) if(mype==mype_io ) then iret=nf90_open(trim(filenamein),nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret write(6,*)' gsi_fv3ncdf_read_v1: problem opening5 with varnam ',trim(varname) return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(dim(ndimensions)) allocate(a(nlat,nlon)) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo allocate(uu(nx,ny,nsig)) allocate(temp0(nx,ny,nsig+1)) 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_get_var(gfile_loc,var_id,temp0) uu(:,:,:)=temp0(:,:,2:(nsig+1)) nztmp=nsig nzp1=nztmp+1 do i=1,nztmp ir=nzp1-i call fv3_h_to_ll(uu(:,:,i),a,nx,ny,nlon,nlat) kk=0 do n=1,npe ns=displss(n)+(ir-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 iret=nf90_close(gfile_loc) deallocate (uu,a,dim) deallocate (temp0) endif !mype call mpi_scatterv(work,ijns,displss,mpi_rtype,& work_sub,ijns(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) deallocate (work) return end subroutine gsi_fv3ncdf_read_v1 subroutine gsi_fv3ncdf_readuv(dynvarsfile,ges_u,ges_v) !$$$ 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: ierror,mpi_comm_world,npe,mpi_rtype,mype use gridmod, only: lat2,lon2,nsig,itotsub,ijn_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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,nya,nxa,fv3uv2earth use general_commvars_mod, only: ltosi_s,ltosj_s implicit none character(*) ,intent(in ):: dynvarsfile real(r_kind) ,intent(out ) :: ges_u(lat2,lon2,nsig) real(r_kind) ,intent(out ) :: ges_v(lat2,lon2,nsig) character(len=128) :: name real(r_kind),allocatable,dimension(:,:,:):: uu,temp1 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(:,:):: u,v integer(i_kind) n,ns,k,len,ndim integer(i_kind) gfile_loc,iret integer(i_kind) nz,nzp1,kk,j,mm1,i,ir,ii,jj integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid allocate (work(itotsub*nsig)) mm1=mype+1 if(mype==mype_u .or. mype==mype_v) then iret=nf90_open(dynvarsfile,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening6 ',trim(dynvarsfile),', Status = ',iret return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(dim(ndimensions)) allocate(a(nya,nxa)) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo allocate(u(dim(1),dim(4))) allocate(v(dim(1),dim(4))) iret=nf90_inq_varid(gfile_loc,trim(adjustl("xaxis_1")),k) !thinkdeb iret=nf90_get_var(gfile_loc,k,u(:,1)) do k=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,k,name,len) if(trim(name)=='u'.or.trim(name)=='U' .or. & trim(name)=='v'.or.trim(name)=='V' ) then iret=nf90_inquire_variable(gfile_loc,k,ndims=ndim) if(allocated(dim_id )) deallocate(dim_id ) allocate(dim_id(ndim)) iret=nf90_inquire_variable(gfile_loc,k,dimids=dim_id) ! NOTE: dimension of variables on native fv3 grid. ! u and v have an extra row in one of the dimensions if(allocated(uu)) deallocate(uu) allocate(uu(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) iret=nf90_get_var(gfile_loc,k,uu) if(trim(name)=='u'.or.trim(name)=='U') then allocate(temp1(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) temp1=uu else if(trim(name)=='v'.or.trim(name)=='V') then exit endif endif enddo ! k ! transfor to earth u/v, interpolate to analysis grid, reverse vertical order nz=nsig nzp1=nz+1 do i=1,nz ir=nzp1-i call fv3uv2earth(temp1(:,:,i),uu(:,:,i),nx,ny,u,v) if(mype==mype_u)then call fv3_h_to_ll(u,a,nx,ny,nxa,nya) else call fv3_h_to_ll(v,a,nx,ny,nxa,nya) endif kk=0 do n=1,npe ns=displss(n)+(ir-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 deallocate(temp1,a) deallocate (dim,dim_id,uu,v,u) iret=nf90_close(gfile_loc) endif ! mype !! scatter to ges_u,ges_v !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns,displss,mpi_rtype,& ges_u,ijns(mm1),mpi_rtype,mype_u,mpi_comm_world,ierror) call mpi_scatterv(work,ijns,displss,mpi_rtype,& ges_v,ijns(mm1),mpi_rtype,mype_v,mpi_comm_world,ierror) deallocate(work) end subroutine gsi_fv3ncdf_readuv subroutine gsi_fv3ncdf_readuv_v1(dynvarsfile,ges_u,ges_v) !$$$ 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: ierror,mpi_comm_world,npe,mpi_rtype,mype use gridmod, only: lat2,lon2,nsig,itotsub,ijn_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 netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,nya,nxa,fv3uv2earth use general_commvars_mod, only: ltosi_s,ltosj_s implicit none character(*) ,intent(in ):: dynvarsfile real(r_kind) ,intent(out ) :: ges_u(lat2,lon2,nsig) real(r_kind) ,intent(out ) :: ges_v(lat2,lon2,nsig) character(len=128) :: name real(r_kind),allocatable,dimension(:,:,:):: uu,temp0 integer(i_kind),allocatable,dimension(:):: dim real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: a real(r_kind),allocatable,dimension(:,:):: uorv integer(i_kind) n,ns,k,len,ndim,var_id integer(i_kind) gfile_loc,iret integer(i_kind) nztmp,nzp1,kk,j,mm1,i,ir,ii,jj integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid allocate (work(itotsub*nsig)) mm1=mype+1 if(mype==mype_u .or. mype==mype_v) then iret=nf90_open(dynvarsfile,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_readuv_v1: problem opening ',trim(dynvarsfile),', Status = ',iret return endif iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid) allocate(dim(ndimensions)) allocate(a(nya,nxa)) do k=1,ndimensions iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo allocate(uorv(nx,ny)) if(mype == mype_u) then allocate(uu(nx,ny+1,nsig)) else ! for mype_v allocate(uu(nx+1,ny,nsig)) endif ! transfor to earth u/v, interpolate to analysis grid, reverse vertical order if(mype == mype_u) then iret=nf90_inq_varid(gfile_loc,trim(adjustl("u_s")),var_id) iret=nf90_inquire_variable(gfile_loc,var_id,ndims=ndim) allocate(temp0(nx,ny+1,nsig+1)) iret=nf90_get_var(gfile_loc,var_id,temp0) uu(:,:,:)=temp0(:,:,2:nsig+1) deallocate(temp0) endif if(mype == mype_v) then allocate(temp0(nx+1,ny,nsig+1)) iret=nf90_inq_varid(gfile_loc,trim(adjustl("v_w")),var_id) iret=nf90_inquire_variable(gfile_loc,var_id,ndims=ndim) iret=nf90_get_var(gfile_loc,var_id,temp0) uu(:,:,:)=(temp0(:,:,2:nsig+1)) deallocate (temp0) endif nztmp=nsig nzp1=nztmp+1 do i=1,nztmp ir=nzp1-i if(mype == mype_u)then do j=1,ny uorv(:,j)=half*(uu(:,j,i)+uu(:,j+1,i)) enddo call fv3_h_to_ll(uorv(:,:),a,nx,ny,nxa,nya) else do j=1,nx uorv(j,:)=half*(uu(j,:,i)+uu(j+1,:,i)) enddo call fv3_h_to_ll(uorv(:,:),a,nx,ny,nxa,nya) endif kk=0 do n=1,npe ns=displss(n)+(ir-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 deallocate(a) deallocate (uu,uorv) iret=nf90_close(gfile_loc) endif ! mype !! scatter to ges_u,ges_v !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns,displss,mpi_rtype,& ges_u,ijns(mm1),mpi_rtype,mype_u,mpi_comm_world,ierror) call mpi_scatterv(work,ijns,displss,mpi_rtype,& ges_v,ijns(mm1),mpi_rtype,mype_v,mpi_comm_world,ierror) deallocate(work) 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: ! ! 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 use mpeu_util, only: die implicit none type (type_fv3regfilenameg),intent(in) :: fv3filenamegin ! Declare local constants logical add_saved character(len=:),allocatable :: dynvars !='fv3_dynvars' character(len=:),allocatable :: tracers !='fv3_tracer' ! 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() dynvars=fv3filenamegin%dynvars tracers=fv3filenamegin%tracers 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 (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(dynvars,'T',ges_tsen(1,1,1,it),mype_t,add_saved) call gsi_fv3ncdf_write(tracers,'sphum',ges_q ,mype_q,add_saved) call gsi_fv3ncdf_writeuv(dynvars,ges_u,ges_v,mype_v,add_saved) call gsi_fv3ncdf_writeps(dynvars,'delp',ges_ps,mype_p,add_saved) else call gsi_fv3ncdf_write(dynvars,'t',ges_tsen(1,1,1,it),mype_t,add_saved) call gsi_fv3ncdf_write(tracers,'sphum',ges_q ,mype_q,add_saved) call gsi_fv3ncdf_writeuv_v1(dynvars,ges_u,ges_v,mype_v,add_saved) call gsi_fv3ncdf_writeps_v1(dynvars,'ps',ges_ps,mype_p,add_saved) endif end subroutine wrfv3_netcdf subroutine gsi_fv3ncdf_writeuv(dynvars,varu,varv,mype_io,add_saved) !$$$ 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,ierror,npe,mype use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1,nsig, & ijn,displs_g,itotsub,iglobal, & nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use general_commvars_mod, only: ltosi,ltosj 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 implicit none character(len=*),intent(in) :: dynvars !='fv3_dynvars' real(r_kind) ,intent(in ) :: varu(lat2,lon2,nsig) real(r_kind) ,intent(in ) :: varv(lat2,lon2,nsig) integer(i_kind),intent(in ) :: mype_io logical ,intent(in ) :: add_saved integer(i_kind) :: ugrd_VarId,gfile_loc,vgrd_VarId integer(i_kind) i,j,mm1,n,k,ns,kr,m real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:,:):: work_sub,work_au,work_av real(r_kind),allocatable,dimension(:,:,:):: work_bu,work_bv real(r_kind),allocatable,dimension(:,:):: u,v,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu2,workbv2 mm1=mype+1 allocate( work(max(iglobal,itotsub)*nsig),work_sub(lat1,lon1,nsig)) !!!!!! gather analysis u !! revers k !!!!!!!!!!! do k=1,nsig kr=nsig+1-k do i=1,lon1 do j=1,lat1 work_sub(j,i,kr)=varu(j+1,i+1,k) end do end do enddo call mpi_gatherv(work_sub,ijnz(mm1),mpi_rtype, & work,ijnz,displsz_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then allocate( work_au(nlat,nlon,nsig),work_av(nlat,nlon,nsig)) ns=0 do m=1,npe do k=1,nsig do n=displs_g(m)+1,displs_g(m)+ijn(m) ns=ns+1 work_au(ltosi(n),ltosj(n),k)=work(ns) end do enddo enddo endif ! mype !!!!!! gather analysis v !! reverse k !!!!!!!!!!!!!!!!!! do k=1,nsig kr=nsig+1-k do i=1,lon1 do j=1,lat1 work_sub(j,i,kr)=varv(j+1,i+1,k) end do end do enddo call mpi_gatherv(work_sub,ijnz(mm1),mpi_rtype, & work,ijnz,displsz_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then ns=0 do m=1,npe do k=1,nsig do n=displs_g(m)+1,displs_g(m)+ijn(m) ns=ns+1 work_av(ltosi(n),ltosj(n),k)=work(ns) end do enddo enddo deallocate(work,work_sub) allocate( u(nlon_regional,nlat_regional+1)) allocate( v(nlon_regional+1,nlat_regional)) allocate( work_bu(nlon_regional,nlat_regional+1,nsig)) allocate( work_bv(nlon_regional+1,nlat_regional,nsig)) call check( nf90_open(trim(dynvars ),nf90_write,gfile_loc) ) 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(nlat,nlon),workav2(nlat,nlon)) allocate( workbu2(nlon_regional,nlat_regional+1)) allocate( workbv2(nlon_regional+1,nlat_regional)) !!!!!!!! readin work_b !!!!!!!!!!!!!!!! call check( nf90_get_var(gfile_loc,ugrd_VarId,work_bu) ) call check( nf90_get_var(gfile_loc,vgrd_VarId,work_bv) ) do k=1,nsig call fv3uv2earth(work_bu(1,1,k),work_bv(1,1,k),nlon_regional,nlat_regional,u,v) call fv3_h_to_ll(u,workau2,nlon_regional,nlat_regional,nlon,nlat) call fv3_h_to_ll(v,workav2,nlon_regional,nlat_regional,nlon,nlat) !!!!!!!! find analysis_inc: work_a !!!!!!!!!!!!!!!! work_au(:,:,k)=work_au(:,:,k)-workau2(:,:) work_av(:,:,k)=work_av(:,:,k)-workav2(:,:) call fv3_ll_to_h(work_au(:,:,k),u,nlon,nlat,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,.true.) call earthuv2fv3(u,v,nlon_regional,nlat_regional,workbu2,workbv2) !!!!!!!! add analysis_inc to readin work_b !!!!!!!!!!!!!!!! work_bu(:,:,k)=work_bu(:,:,k)+workbu2(:,:) work_bv(:,:,k)=work_bv(:,:,k)+workbv2(:,:) enddo deallocate(workau2,workbu2,workav2,workbv2) else do k=1,nsig call fv3_ll_to_h(work_au(:,:,k),u,nlon,nlat,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,.true.) call earthuv2fv3(u,v,nlon_regional,nlat_regional,work_bu(:,:,k),work_bv(:,:,k)) enddo endif deallocate(work_au,work_av,u,v) print *,'write out u/v to ',trim(dynvars ) call check( nf90_put_var(gfile_loc,ugrd_VarId,work_bu) ) call check( nf90_put_var(gfile_loc,vgrd_VarId,work_bv) ) call check( nf90_close(gfile_loc) ) deallocate(work_bu,work_bv) end if !mype_io if(allocated(work))deallocate(work) if(allocated(work_sub))deallocate(work_sub) end subroutine gsi_fv3ncdf_writeuv subroutine gsi_fv3ncdf_writeps(filename,varname,var,mype_io,add_saved) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeps ! pgrmmr: wu ! ! abstract: write out analyzed "delp" to fv_core.res.nest02.tile7.nc ! ! 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,ierror,mype use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1,nsig use gridmod, only: ijn,displs_g,itotsub,iglobal use gridmod, only: nlon_regional,nlat_regional,eta1_ll,eta2_ll use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll use general_commvars_mod, only: ltosi,ltosj use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var implicit none real(r_kind) ,intent(in ) :: var(lat2,lon2) integer(i_kind),intent(in ) :: mype_io logical ,intent(in ) :: add_saved character(*) ,intent(in ) :: varname,filename integer(i_kind) :: VarId,gfile_loc integer(i_kind) i,j,mm1,k,kr,kp real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: work_sub,work_a real(r_kind),allocatable,dimension(:,:,:):: work_b,work_bi real(r_kind),allocatable,dimension(:,:):: workb2,worka2 mm1=mype+1 allocate( work(max(iglobal,itotsub)),work_sub(lat1,lon1) ) do i=1,lon1 do j=1,lat1 work_sub(j,i)=var(j+1,i+1) end do end do call mpi_gatherv(work_sub,ijn(mm1),mpi_rtype, & work,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then allocate( work_a(nlat,nlon)) do i=1,iglobal work_a(ltosi(i),ltosj(i))=work(i) end do allocate( work_bi(nlon_regional,nlat_regional,nsig+1)) allocate( work_b(nlon_regional,nlat_regional,nsig)) call check( nf90_open(trim(filename),nf90_write,gfile_loc) ) call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) if(add_saved)then allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlat,nlon)) !!!!!!!! read in guess delp !!!!!!!!!!!!!! call check( nf90_get_var(gfile_loc,VarId,work_b) ) work_bi(:,:,1)=eta1_ll(nsig+1) do i=2,nsig+1 work_bi(:,:,i)=work_b(:,:,i-1)*0.001_r_kind+work_bi(:,:,i-1) enddo call fv3_h_to_ll(work_bi(:,:,nsig+1),worka2,nlon_regional,nlat_regional,nlon,nlat) !!!!!!! analysis_inc Psfc: work_a work_a(:,:)=work_a(:,:)-worka2(:,:) call fv3_ll_to_h(work_a,workb2,nlon,nlat,nlon_regional,nlat_regional,.true.) do k=1,nsig+1 kr=nsig+2-k !!!!!!! ges_prsi+hydrostatic analysis_inc !!!!!!!!!!!!!!!! work_bi(:,:,k)=work_bi(:,:,k)+eta2_ll(kr)*workb2(:,:) enddo else call fv3_ll_to_h(work_a,workb2,nlon,nlat,nlon_regional,nlat_regional,.true.) do k=1,nsig+1 kr=nsig+2-k !!!!!!! Psfc_ges+hydrostatic analysis_inc !!!!!!!!!!!!!!!! work_bi(:,:,k)=eta1_ll(kr)+eta2_ll(kr)*workb2(:,:) enddo endif ! delp do k=nsig,1,-1 kp=k+1 work_b(:,:,k)=(work_bi(:,:,kp)-work_bi(:,:,k))*1000._r_kind enddo call check( nf90_put_var(gfile_loc,VarId,work_b) ) call check( nf90_close(gfile_loc) ) deallocate(worka2,workb2) deallocate(work_b,work_a,work_bi) end if !mype_io deallocate(work,work_sub) end subroutine gsi_fv3ncdf_writeps subroutine gsi_fv3ncdf_writeuv_v1(dynvars,varu,varv,mype_io,add_saved) !$$$ 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! ! 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,ierror,npe,mype use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1,nsig, & ijn,displs_g,itotsub,iglobal, & nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use general_commvars_mod, only: ltosi,ltosj 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 implicit none character(len=*),intent(in) :: dynvars !='fv3_dynvars' real(r_kind) ,intent(in ) :: varu(lat2,lon2,nsig) real(r_kind) ,intent(in ) :: varv(lat2,lon2,nsig) integer(i_kind),intent(in ) :: mype_io logical ,intent(in ) :: add_saved 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,n,k,ns,kr,m real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:,:):: work_sub,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(:,:):: u,v,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu_s2,workbv_s2 real(r_kind),allocatable,dimension(:,:):: workbu_w2,workbv_w2 mm1=mype+1 allocate( work(max(iglobal,itotsub)*nsig),work_sub(lat1,lon1,nsig)) !!!!!! gather analysis u !! revers k !!!!!!!!!!! do k=1,nsig kr=nsig+1-k do i=1,lon1 do j=1,lat1 work_sub(j,i,kr)=varu(j+1,i+1,k) end do end do enddo call mpi_gatherv(work_sub,ijnz(mm1),mpi_rtype, & work,ijnz,displsz_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then allocate( work_au(nlat,nlon,nsig),work_av(nlat,nlon,nsig)) ns=0 do m=1,npe do k=1,nsig do n=displs_g(m)+1,displs_g(m)+ijn(m) ns=ns+1 work_au(ltosi(n),ltosj(n),k)=work(ns) end do enddo enddo endif ! mype !!!!!! gather analysis v !! reverse k !!!!!!!!!!!!!!!!!! do k=1,nsig kr=nsig+1-k do i=1,lon1 do j=1,lat1 work_sub(j,i,kr)=varv(j+1,i+1,k) end do end do enddo call mpi_gatherv(work_sub,ijnz(mm1),mpi_rtype, & work,ijnz,displsz_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then ns=0 do m=1,npe do k=1,nsig do n=displs_g(m)+1,displs_g(m)+ijn(m) ns=ns+1 work_av(ltosi(n),ltosj(n),k)=work(ns) end do enddo enddo deallocate(work,work_sub) !clt u and v would contain winds at either D-grid or A-grid !clt do not diretly use them in between fv3uv2eath and fv3_h_to_ll unless paying !attention to the actual storage layout call check( nf90_open(trim(dynvars ),nf90_write,gfile_loc) ) allocate( u(nlon_regional,nlat_regional)) allocate( v(nlon_regional,nlat_regional)) allocate( work_bu_s(nlon_regional,nlat_regional+1,nsig)) allocate( work_bv_s(nlon_regional,nlat_regional+1,nsig)) allocate( work_bu_w(nlon_regional+1,nlat_regional,nsig)) allocate( work_bv_w(nlon_regional+1,nlat_regional,nsig)) 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) ) if(add_saved)then allocate( workau2(nlat,nlon),workav2(nlat,nlon)) 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)) !!!!!!!! readin work_b !!!!!!!!!!!!!!!! call check( nf90_get_var(gfile_loc,u_sgrd_VarId,work_bu_s) ) call check( nf90_get_var(gfile_loc,u_wgrd_VarId,work_bu_w) ) call check( nf90_get_var(gfile_loc,v_sgrd_VarId,work_bv_s) ) call check( nf90_get_var(gfile_loc,v_wgrd_VarId,work_bv_w) ) do k=1,nsig do j=1,nlat_regional u(:,j)=half * (work_bu_s(:,j,k)+ work_bu_s(:,j+1,k)) enddo do i=1,nlon_regional v(i,:)=half*(work_bv_w(i,:,k)+work_bv_w(i+1,:,k)) enddo call fv3_h_to_ll(u,workau2,nlon_regional,nlat_regional,nlon,nlat) call fv3_h_to_ll(v,workav2,nlon_regional,nlat_regional,nlon,nlat) !!!!!!!! find analysis_inc: work_a !!!!!!!!!!!!!!!! work_au(:,:,k)=work_au(:,:,k)-workau2(:,:) work_av(:,:,k)=work_av(:,:,k)-workav2(:,:) call fv3_ll_to_h(work_au(:,:,k),u,nlon,nlat,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,.true.) !!!!!!!! add analysis_inc to readin work_b !!!!!!!!!!!!!!!! do i=2,nlon_regional-1 workbu_w2(i,:)=half*(u(i,:)+u(i+1,:)) workbv_w2(i,:)=half*(v(i,:)+v(i+1,:)) enddo workbu_w2(1,:)=u(1,:) workbv_w2(1,:)=v(1,:) workbu_w2(nlon_regional+1,:)=u(nlon_regional,:) workbv_w2(nlon_regional+1,:)=v(nlon_regional,:) do j=2,nlat_regional-1 workbu_s2(:,j)=half*(u(:,j)+u(:,j+1)) workbv_s2(:,j)=half*(v(:,j)+v(:,j+1)) enddo workbu_s2(:,1)=u(:,1) workbv_s2(:,1)=v(:,1) workbu_s2(:,nlat_regional+1)=u(:,nlat_regional) workbv_s2(:,nlat_regional+1)=v(:,nlat_regional) work_bu_w(:,:,k)=work_bu_w(:,:,k)+workbu_w2(:,:) work_bu_s(:,:,k)=work_bu_s(:,:,k)+workbu_s2(:,:) work_bv_w(:,:,k)=work_bv_w(:,:,k)+workbv_w2(:,:) work_bv_s(:,:,k)=work_bv_s(:,:,k)+workbv_s2(:,:) enddo deallocate(workau2,workav2) deallocate(workbu_w2,workbv_w2) deallocate(workbu_s2,workbv_s2) else do k=1,nsig call fv3_ll_to_h(work_au(:,:,k),u,nlon,nlat,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,.true.) do i=2,nlon_regional-1 work_bu_w(i,:,k)=half*(u(i,:)+u(i+1,:)) work_bv_w(i,:,k)=half*(v(i,:)+v(i+1,:)) enddo work_bu_w(1,:,k)=u(1,:) work_bv_w(1,:,k)=v(1,:) work_bu_w(nlon_regional+1,:,k)=u(nlon_regional,:) work_bv_w(nlon_regional+1,:,k)=v(nlon_regional,:) do j=2,nlat_regional-1 work_bu_s(:,j,k)=half*(u(:,j)+u(:,j+1)) work_bv_s(:,j,k)=half*(v(:,j)+v(:,j+1)) enddo work_bu_s(:,1,k)=u(:,1) work_bv_s(:,1,k)=v(:,1) work_bu_s(:,nlat_regional+1,k)=u(:,nlat_regional) work_bv_s(:,nlat_regional+1,k)=v(:,nlat_regional) enddo endif deallocate(work_au,work_av,u,v) print *,'write out u/v to ',trim(dynvars ) call check( nf90_put_var(gfile_loc,u_wgrd_VarId,work_bu_w) ) call check( nf90_put_var(gfile_loc,u_sgrd_VarId,work_bu_s) ) call check( nf90_put_var(gfile_loc,v_wgrd_VarId,work_bv_w) ) call check( nf90_put_var(gfile_loc,v_sgrd_VarId,work_bv_s) ) call check( nf90_close(gfile_loc) ) deallocate(work_bu_w,work_bv_w) deallocate(work_bu_s,work_bv_s) end if !mype_io if(allocated(work))deallocate(work) if(allocated(work_sub))deallocate(work_sub) end subroutine gsi_fv3ncdf_writeuv_v1 subroutine gsi_fv3ncdf_writeps_v1(filename,varname,var,mype_io,add_saved) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeps ! pgrmmr: wu ! ! abstract: write out analyzed "delp" to fv_core.res.nest02.tile7.nc ! ! program history log: ! 2019-04 lei, modified from gsi_nemsio_writeps to deal with cold start files ! ! 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,ierror,mype use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1 use gridmod, only: ijn,displs_g,itotsub,iglobal use gridmod, only: nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll use general_commvars_mod, only: ltosi,ltosj use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var implicit none real(r_kind) ,intent(in ) :: var(lat2,lon2) integer(i_kind),intent(in ) :: mype_io logical ,intent(in ) :: add_saved character(*) ,intent(in ) :: varname,filename integer(i_kind) :: VarId,gfile_loc integer(i_kind) i,j,mm1 real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:):: work_sub,work_a real(r_kind),allocatable,dimension(:,:):: work_b,work_bi real(r_kind),allocatable,dimension(:,:):: workb2,worka2 mm1=mype+1 allocate( work(max(iglobal,itotsub)),work_sub(lat1,lon1) ) do i=1,lon1 do j=1,lat1 work_sub(j,i)=var(j+1,i+1) end do end do call mpi_gatherv(work_sub,ijn(mm1),mpi_rtype, & work,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then allocate( work_a(nlat,nlon)) do i=1,iglobal work_a(ltosi(i),ltosj(i))=work(i) end do allocate( work_bi(nlon_regional,nlat_regional)) allocate( work_b(nlon_regional,nlat_regional)) call check( nf90_open(trim(filename),nf90_write,gfile_loc) ) call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) work_a=work_a*1000.0_r_kind if(add_saved)then allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlat,nlon)) !!!!!!!! read in guess delp !!!!!!!!!!!!!! call check( nf90_get_var(gfile_loc,VarId,work_b) ) call fv3_h_to_ll(work_b,worka2,nlon_regional,nlat_regional,nlon,nlat) !!!!!!! analysis_inc Psfc: work_a work_a(:,:)=work_a(:,:)-worka2(:,:) call fv3_ll_to_h(work_a,workb2,nlon,nlat,nlon_regional,nlat_regional,.true.) work_b(:,:)=work_b(:,:)+workb2(:,:) else call fv3_ll_to_h(work_a,work_b,nlon,nlat,nlon_regional,nlat_regional,.true.) endif call check( nf90_put_var(gfile_loc,VarId,work_b) ) call check( nf90_close(gfile_loc) ) deallocate(worka2,workb2) deallocate(work_b,work_a,work_bi) end if !mype_io deallocate(work,work_sub) end subroutine gsi_fv3ncdf_writeps_v1 subroutine gsi_fv3ncdf_write(filename,varname,var,mype_io,add_saved) !$$$ 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,ierror,npe,mype use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1,nsig use gridmod, only: ijn,displs_g,itotsub,iglobal use mod_fv3_lola, only: fv3_ll_to_h use mod_fv3_lola, only: fv3_h_to_ll use general_commvars_mod, only: ltosi,ltosj use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var implicit none real(r_kind) ,intent(in ) :: var(lat2,lon2,nsig) integer(i_kind),intent(in ) :: mype_io logical ,intent(in ) :: add_saved character(*) ,intent(in ) :: varname,filename integer(i_kind) :: VarId,gfile_loc integer(i_kind) i,j,mm1,k,kr,ns,n,m real(r_kind),allocatable,dimension(:):: work real(r_kind),allocatable,dimension(:,:,:):: work_sub,work_a real(r_kind),allocatable,dimension(:,:,:):: work_b real(r_kind),allocatable,dimension(:,:):: workb2,worka2 mm1=mype+1 allocate( work(max(iglobal,itotsub)*nsig),work_sub(lat1,lon1,nsig)) !!!!!!!! reverse z !!!!!!!!!!!!!! do k=1,nsig kr=nsig+1-k do i=1,lon1 do j=1,lat1 work_sub(j,i,kr)=var(j+1,i+1,k) end do end do enddo call mpi_gatherv(work_sub,ijnz(mm1),mpi_rtype, & work,ijnz,displsz_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then allocate( work_a(nlat,nlon,nsig)) ns=0 do m=1,npe do k=1,nsig do n=displs_g(m)+1,displs_g(m)+ijn(m) ns=ns+1 work_a(ltosi(n),ltosj(n),k)=work(ns) end do enddo enddo allocate( work_b(nlon_regional,nlat_regional,nsig)) call check( nf90_open(trim(filename),nf90_write,gfile_loc) ) call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) if(add_saved)then allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlat,nlon)) call check( nf90_get_var(gfile_loc,VarId,work_b) ) do k=1,nsig call fv3_h_to_ll(work_b(:,:,k),worka2,nlon_regional,nlat_regional,nlon,nlat) !!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! work_a(:,:,k)=work_a(:,:,k)-worka2(:,:) call fv3_ll_to_h(work_a(1,1,k),workb2,nlon,nlat,nlon_regional,nlat_regional,.true.) work_b(:,:,k)=work_b(:,:,k)+workb2(:,:) enddo deallocate(worka2,workb2) else do k=1,nsig call fv3_ll_to_h(work_a(1,1,k),work_b(1,1,k),nlon,nlat,nlon_regional,nlat_regional,.true.) enddo endif print *,'write out ',trim(varname),' to ',trim(filename) call check( nf90_put_var(gfile_loc,VarId,work_b) ) call check( nf90_close(gfile_loc) ) deallocate(work_b,work_a) end if !mype_io deallocate(work,work_sub) 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 end module gsi_rfv3io_mod