From a21fb6610d1bf5370f318995141bbd5a4b7753e4 Mon Sep 17 00:00:00 2001 From: TingLei-daprediction Date: Fri, 5 Nov 2021 13:22:27 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#218 GSI parallel IO for FV3-LAM --- src/enkf/gridio_fv3reg.f90 | 1 + src/gsi/cplr_get_fv3_regional_ensperts.f90 | 1162 +++++---- src/gsi/gsi_rfv3io_mod.f90 | 2694 ++++++++++---------- src/gsi/gsimod.F90 | 4 + src/gsi/mod_fv3_lola.f90 | 7 + 5 files changed, 1953 insertions(+), 1915 deletions(-) diff --git a/src/enkf/gridio_fv3reg.f90 b/src/enkf/gridio_fv3reg.f90 index cfd8b35b26..0dd03dfb82 100644 --- a/src/enkf/gridio_fv3reg.f90 +++ b/src/enkf/gridio_fv3reg.f90 @@ -591,6 +591,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid tsen_ind = getindex(vars3d, 'tsen') ! Tv (3D) q_ind = getindex(vars3d, 'q') ! Q (3D) cw_ind = getindex(vars3d, 'cw') ! CWM for WRF-NMM + oz_ind = getindex(vars3d, 'oz') ! Oz (3D) w_ind = getindex(vars3d, 'w') ! W for WRF-ARW ph_ind = getindex(vars3d, 'ph') ! PH for WRF-ARW diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index 1d79aef6e0..b54f727b71 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -1,14 +1,23 @@ module get_fv3_regional_ensperts_mod use abstract_get_fv3_regional_ensperts_mod,only: abstract_get_fv3_regional_ensperts_class use kinds, only : i_kind + use general_sub2grid_mod, only: sub2grid_info + use constants, only:max_varname_length type, extends(abstract_get_fv3_regional_ensperts_class) :: get_fv3_regional_ensperts_class contains procedure, pass(this) :: get_fv3_regional_ensperts => get_fv3_regional_ensperts_run procedure, pass(this) :: ens_spread_dualres_regional => ens_spread_dualres_regional_fv3_regional procedure, pass(this) :: general_read_fv3_regional - procedure, pass(this) :: general_read_fv3_regional_dirZDA end type get_fv3_regional_ensperts_class + type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_uv + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars3d_nouv + ! copy of cvars3d excluding uv 3-d fields + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars3d_nouv + ! copy of cvars3d excluding uv 3-d fields + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars2d_nouv + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars2d_nouv contains + subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) !$$$ subprogram documentation block ! . . . . @@ -24,6 +33,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! 2019-04-22 CAPS(C. Tong) - add direct reflectivity DA option ! ! 2021-08-10 lei - modify for fv3-lam ensemble spread output + ! 2021-11-01 lei - modify for fv3-lam parallel IO ! input argument list: ! ! output argument list: @@ -34,372 +44,477 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! !$$$ end documentation block - use kinds, only: r_kind,i_kind,r_single - use constants, only: zero,one,half,zero_single,rd_over_cp,one_tenth - use mpimod, only: mpi_comm_world,ierror,mype - use hybrid_ensemble_parameters, only: n_ens,grd_ens - use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path - use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d - use gsi_bundlemod, only: gsi_bundlecreate - use gsi_bundlemod, only: gsi_grid - use gsi_bundlemod, only: gsi_bundle - use gsi_bundlemod, only: gsi_bundlegetpointer - use gsi_bundlemod, only: gsi_bundledestroy - use gsi_bundlemod, only: gsi_gridcreate - use gsi_4dvar, only: ens_fhrlevs - use gsi_rfv3io_mod, only: type_fv3regfilenameg - use hybrid_ensemble_parameters, only: write_ens_sprd - 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 - - implicit none - class(get_fv3_regional_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) - integer(i_kind), intent(in ):: nelen - real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar - - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr - - real(r_single),pointer,dimension(:,:,:):: w3 - real(r_single),pointer,dimension(:,:):: w2 - real(r_kind),pointer,dimension(:,:,:):: x3 - real(r_kind),pointer,dimension(:,:):: x2 - type(gsi_bundle),allocatable,dimension(:):: en_bar - type(gsi_grid):: grid_ens - real(r_kind):: bar_norm,sig_norm,kapr,kap1 - - integer(i_kind):: i,j,k,n,mm1,istatus - integer(i_kind):: ic2,ic3 - integer(i_kind):: m - - - character(255) ensfilenam_str - type(type_fv3regfilenameg)::fv3_filename - - call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) - ! Allocate bundle to hold mean of ensemble members - allocate(en_bar(ntlevs_ens)) - do m=1,ntlevs_ens - call gsi_bundlecreate(en_bar(m),grid_ens,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind) - if(istatus/=0) then - write(6,*)' get_fv3_regional_ensperts_netcdf: trouble creating en_bar bundle' - call stop2(9991) - endif - enddo ! for m + use kinds, only: r_kind,i_kind,r_single + use constants, only: zero,one,half,zero_single,rd_over_cp,one_tenth + use mpimod, only: mpi_comm_world,ierror,mype + use hybrid_ensemble_parameters, only: n_ens,grd_ens + use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path + use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d + use gsi_bundlemod, only: gsi_bundlecreate + use gsi_bundlemod, only: gsi_grid + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use gsi_bundlemod, only: gsi_bundledestroy + use gsi_bundlemod, only: gsi_gridcreate + use gsi_4dvar, only: ens_fhrlevs + use gsi_rfv3io_mod, only: type_fv3regfilenameg + use hybrid_ensemble_parameters, only: write_ens_sprd + 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 + use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info + use gridmod,only: regional + use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv + use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv + + + implicit none + class(get_fv3_regional_ensperts_class), intent(inout) :: this + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + integer(i_kind), intent(in ):: nelen + real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar - ! print info message for dirZDA - if(mype==0)then - if (l_use_cvpqx) then - if ( cvpqx_pval == 0._r_kind ) then ! CVlogq - write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg to log transform.' - else if ( cvpqx_pval > 0._r_kind ) then ! CVpq - write(6,*) 'general_read_fv3_regional_dirZDA: ', & - 'reset minimum of qr/qs/qg to specified values before power transform.' - write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg with power transform.' - end if - end if - if ( cld_nt_updt > 0 .and. l_cvpnr) then - write(6,*) 'general_read_fv3_regional_dirZDA: ', & - 'reset minimum of qnr to a specified value before power transform' - write(6,*) 'general_read_fv3_regional_dirZDA: convert qnr with power transform .' - end if - end if + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr + real(r_single),pointer,dimension(:,:,:):: w3 =>NULL() + real(r_single),pointer,dimension(:,:):: w2 =>NULL() + real(r_kind),pointer,dimension(:,:,:):: x3 =>NULL() + real(r_kind),pointer,dimension(:,:):: x2 =>NULL() + type(gsi_bundle),allocatable,dimension(:):: en_bar + type(gsi_grid):: grid_ens + real(r_kind):: bar_norm,sig_norm,kapr,kap1 + + character(len=64),dimension(:,:),allocatable:: names + character(len=64),dimension(:,:),allocatable:: uvnames + integer(i_kind),dimension(:,:),allocatable:: lnames + integer(i_kind),dimension(:,:),allocatable:: uvlnames + + integer(i_kind):: i,j,k,n,mm1,istatus + integer(i_kind):: ndynvario2d,ntracerio2d + integer(r_kind):: ndynvario3d,ntracerio3d + integer(i_kind):: inner_vars,numfields + integer(i_kind):: ilev,ic2,ic3 + integer(i_kind):: m - do m=1,ntlevs_ens + + character(255) ensfilenam_str + type(type_fv3regfilenameg)::fv3_filename + + call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) + ! Allocate bundle to hold mean of ensemble members + allocate(en_bar(ntlevs_ens)) + +!clt setup varnames for IO + ndynvario2d=0 + ntracerio2d=0 + ndynvario3d=size(fv3lam_io_dynmetvars3d_nouv) + ntracerio3d=size(fv3lam_io_tracermetvars3d_nouv) + if (allocated(fv3lam_io_dynmetvars2d_nouv))then + ndynvario2d=size(fv3lam_io_dynmetvars2d_nouv) + endif + if (allocated(fv3lam_io_tracermetvars2d_nouv)) then + ntracerio2d=size(fv3lam_io_tracermetvars2d_nouv) + endif + allocate(fv3lam_ens_io_dynmetvars3d_nouv(ndynvario3d)) + allocate(fv3lam_ens_io_tracermetvars3d_nouv(ndynvario3d)) + fv3lam_ens_io_dynmetvars3d_nouv=fv3lam_io_dynmetvars3d_nouv + fv3lam_ens_io_tracermetvars3d_nouv=fv3lam_io_tracermetvars3d_nouv + if (ndynvario2d > 0 ) then + allocate(fv3lam_ens_io_dynmetvars2d_nouv(ndynvario2d)) + fv3lam_ens_io_dynmetvars2d_nouv=fv3lam_io_dynmetvars2d_nouv + endif + if (ntracerio2d > 0 ) then + allocate(fv3lam_ens_io_tracermetvars2d_nouv(ntracerio2d)) + fv3lam_ens_io_tracermetvars2d_nouv =fv3lam_io_tracermetvars3d_nouv + endif + inner_vars=1 + numfields=inner_vars*(ndynvario3d*grd_ens%nsig+ndynvario2d) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,ndynvario3d + do k=1,grd_ens%nsig + lnames(1,ilev)=k + names(1,ilev)=fv3lam_ens_io_dynmetvars3d_nouv(i) + ilev=ilev+1 + enddo + enddo + do i=1,ndynvario2d + lnames(1,ilev)=1 + names(1,ilev)=fv3lam_ens_io_dynmetvars2d_nouv(i) + ilev=ilev+1 + enddo + - ! - ! INITIALIZE ENSEMBLE MEAN ACCUMULATORS - en_bar(m)%values=zero - - do n=1,n_ens - en_perts(n,m)%valuesr4 = zero - enddo - - mm1=mype+1 - kap1=rd_over_cp+one - kapr=one/rd_over_cp - ! - ! LOOP OVER ENSEMBLE MEMBERS - do n=1,n_ens - write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n + call general_sub2grid_create_info(grd_fv3lam_ens_dynvar_io_nouv,inner_vars,grd_ens%nlat,& + grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + + inner_vars=1 + numfields=inner_vars*(ntracerio3d*grd_ens%nsig+ntracerio2d) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,ntracerio3d + do k=1,grd_ens%nsig + lnames(1,ilev)=k + names(1,ilev)=fv3lam_ens_io_tracermetvars3d_nouv(i) + ilev=ilev+1 + enddo + enddo + do i=1,ntracerio2d + lnames(1,ilev)=1 + names(1,ilev)=fv3lam_ens_io_tracermetvars2d_nouv(i) + ilev=ilev+1 + enddo + + + call general_sub2grid_create_info(grd_fv3lam_ens_tracer_io_nouv,inner_vars,grd_ens%nlat,& + grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + + + + + + numfields=grd_ens%nsig + inner_vars=2 + allocate(uvlnames(inner_vars,numfields),uvnames(inner_vars,numfields)) + ilev=1 + do k=1,grd_ens%nsig + uvlnames(1,ilev)=k + uvlnames(2,ilev)=k + uvnames(1,ilev)='u' + uvnames(2,ilev)='v' + ilev=ilev+1 + enddo + call general_sub2grid_create_info(grd_fv3lam_ens_uv,inner_vars,grd_ens%nlat,& + grd_ens%nlon,grd_ens%nsig,numfields,regional,names=uvnames,lnames=uvlnames) + + deallocate(lnames,names,uvlnames,uvnames) + + + do m=1,ntlevs_ens + call gsi_bundlecreate(en_bar(m),grid_ens,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind) + if(istatus/=0) then + write(6,*)' get_fv3_regional_ensperts_netcdf: trouble creating en_bar bundle' + call stop2(9991) + endif + enddo ! for m + + ! print info message for dirZDA + if(mype==0)then + if (l_use_cvpqx) then + if ( cvpqx_pval == 0._r_kind ) then ! CVlogq + write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg to log transform.' + else if ( cvpqx_pval > 0._r_kind ) then ! CVpq + write(6,*) 'general_read_fv3_regional_dirZDA: ', & + 'reset minimum of qr/qs/qg to specified values before power transform.' + write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg with power transform.' + end if + end if + if ( cld_nt_updt > 0 .and. l_cvpnr) then + write(6,*) 'general_read_fv3_regional_dirZDA: ', & + 'reset minimum of qnr to a specified value before power transform' + write(6,*) 'general_read_fv3_regional_dirZDA: convert qnr with power transform .' + end if + end if + + + do m=1,ntlevs_ens + + + + ! + ! INITIALIZE ENSEMBLE MEAN ACCUMULATORS + en_bar(m)%values=zero + + do n=1,n_ens + en_perts(n,m)%valuesr4 = zero + enddo + + mm1=mype+1 + kap1=rd_over_cp+one + kapr=one/rd_over_cp + ! + ! LOOP OVER ENSEMBLE MEMBERS + do n=1,n_ens + write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n 22 format(a,'fv3SAR',i2.2,'_ens_mem',i3.3) - ! DEFINE INPUT FILE NAME - fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' !exmaple thinktobe - fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk' - fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars' - fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer" - fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" - fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" - ! - ! READ ENEMBLE MEMBERS DATA - if (mype == 0) write(6,'(a,a)') 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) - - if ( l_use_dbz_directDA ) then ! Read additional hydrometers and w for dirZDA - call this%general_read_fv3_regional_dirZDA(fv3_filename,ql,qi,qr,qs,qg,qnr,w) - end if - - ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR - do ic3=1,nc3d - - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n - call stop2(9992) - end if - call gsi_bundlegetpointer(en_bar(m),trim(cvars3d(ic3)),x3,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for en_bar' - call stop2(9993) - end if - - select case (trim(cvars3d(ic3))) - - case('sf','SF') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = u(j,i,k) - x3(j,i,k)=x3(j,i,k)+u(j,i,k) - end do - end do - end do - - case('vp','VP') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = v(j,i,k) - x3(j,i,k)=x3(j,i,k)+v(j,i,k) - end do - end do - end do - - case('t','T') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = tv(j,i,k) - x3(j,i,k)=x3(j,i,k)+tv(j,i,k) - end do - end do - end do - - case('q','Q') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = rh(j,i,k) - x3(j,i,k)=x3(j,i,k)+rh(j,i,k) - end do - end do - end do - - case('oz','OZ') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = oz(j,i,k) - x3(j,i,k)=x3(j,i,k)+oz(j,i,k) - end do - end do - end do - + ! DEFINE INPUT FILE NAME + fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' !exmaple thinktobe + fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk' + fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars' + fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer" + fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" + fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" + ! + ! READ ENEMBLE MEMBERS DATA + if (mype == 0) write(6,'(a,a)') & + 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) + if (.not.l_use_dbz_directDA ) then ! Read additional hydrometers and w for dirZDA + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) + else + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + end if + + ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR + do ic3=1,nc3d + + call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n + call stop2(9992) + end if + call gsi_bundlegetpointer(en_bar(m),trim(cvars3d(ic3)),x3,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for en_bar' + call stop2(9993) + end if + + select case (trim(cvars3d(ic3))) + + case('sf','SF') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = u(j,i,k) + x3(j,i,k)=x3(j,i,k)+u(j,i,k) + end do + end do + end do + + case('vp','VP') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = v(j,i,k) + x3(j,i,k)=x3(j,i,k)+v(j,i,k) + end do + end do + end do + + case('t','T') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = tv(j,i,k) + x3(j,i,k)=x3(j,i,k)+tv(j,i,k) + end do + end do + end do + + case('q','Q') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = rh(j,i,k) + x3(j,i,k)=x3(j,i,k)+rh(j,i,k) + end do + end do + end do + + case('oz','OZ') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = oz(j,i,k) + x3(j,i,k)=x3(j,i,k)+oz(j,i,k) + end do + end do + end do + ! save additional ensemble varaible data for direct reflectivity DA - case('ql','QL') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = ql(j,i,k) - x3(j,i,k)=x3(j,i,k)+ql(j,i,k) - end do - end do - end do - - case('qi','QI') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = qi(j,i,k) - x3(j,i,k)=x3(j,i,k)+qi(j,i,k) - end do - end do - end do - - case('qr','QR') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = qr(j,i,k) - x3(j,i,k)=x3(j,i,k)+qr(j,i,k) - end do - end do - end do - - case('qs','QS') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = qs(j,i,k) - x3(j,i,k)=x3(j,i,k)+qs(j,i,k) - end do - end do - end do - - case('qg','QG') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = qg(j,i,k) - x3(j,i,k)=x3(j,i,k)+qg(j,i,k) - end do - end do - end do - - case('qnr','QNR') - - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - if ( l_use_dbz_directDA ) then ! direct reflectivity DA - if ( cld_nt_updt > 0 ) then ! Update Nc - w3(j,i,k) = qnr(j,i,k) - x3(j,i,k)=x3(j,i,k)+qnr(j,i,k) - end if - else ! .not. l_use_dbz_directDA - w3(j,i,k) = qnr(j,i,k) - x3(j,i,k)=x3(j,i,k)+qnr(j,i,k) - end if - end do + case('ql','QL') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = ql(j,i,k) + x3(j,i,k)=x3(j,i,k)+ql(j,i,k) end do end do + end do - case('w','W') - do k=1,grd_ens%nsig - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w3(j,i,k) = w(j,i,k) - x3(j,i,k)=x3(j,i,k)+w(j,i,k) - end do - end do - end do - end select + case('qi','QI') - + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = qi(j,i,k) + x3(j,i,k)=x3(j,i,k)+qi(j,i,k) + end do + end do + end do - end do - - do ic2=1,nc2d - - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n - call stop2(9994) - end if - call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar' - call stop2(9995) - end if - - select case (trim(cvars2d(ic2))) - - case('ps','PS') - - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w2(j,i) = ps(j,i) - x2(j,i)=x2(j,i)+ps(j,i) - end do - end do - - case('sst','SST') - ! IGNORE SST IN HYBRID for now - - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - w2(j,i) = zero - x2(j,i)=zero - end do - end do - - end select - end do - enddo - ! - ! CALCULATE ENSEMBLE MEAN - bar_norm = one/float(n_ens) - en_bar(m)%values=en_bar(m)%values*bar_norm - - ! Copy pbar to module array. ps_bar may be needed for vertical localization - ! in terms of scale heights/normalized p/p - do ic2=1,nc2d - - if(trim(cvars2d(ic2)) == 'ps'.or.trim(cvars2d(ic2)) == 'PS') then - - call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar to get ps_bar' - call stop2(9996) - end if - - do i=1,grd_ens%lon2 - do j=1,grd_ens%lat2 - ps_bar(j,i,1)=x2(j,i) - end do - end do - exit - end if - end do - - call mpi_barrier(mpi_comm_world,ierror) - ! - ! - ! CONVERT ENSEMBLE MEMBERS TO ENSEMBLE PERTURBATIONS - sig_norm=sqrt(one/max(one,n_ens-one)) - - do n=1,n_ens - do i=1,nelen - en_perts(n,m)%valuesr4(i)=(en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm - end do - end do + case('qr','QR') - enddo ! it 4d loop - ! CALCULATE ENSEMBLE SPREAD - write_ens_sprd=.true. - if(write_ens_sprd ) then - call this%ens_spread_dualres_regional(mype,en_perts,nelen) - call mpi_barrier(mpi_comm_world,ierror) ! do we need this mpi_barrier here? - endif - do m=1,ntlevs_ens + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = qr(j,i,k) + x3(j,i,k)=x3(j,i,k)+qr(j,i,k) + end do + end do + end do + + case('qs','QS') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = qs(j,i,k) + x3(j,i,k)=x3(j,i,k)+qs(j,i,k) + end do + end do + end do + + case('qg','QG') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = qg(j,i,k) + x3(j,i,k)=x3(j,i,k)+qg(j,i,k) + end do + end do + end do + + case('qnr','QNR') + + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + if ( l_use_dbz_directDA ) then ! direct reflectivity DA + if ( cld_nt_updt > 0 ) then ! Update Nc + w3(j,i,k) = qnr(j,i,k) + x3(j,i,k)=x3(j,i,k)+qnr(j,i,k) + end if + else ! .not. l_use_dbz_directDA + w3(j,i,k) = qnr(j,i,k) + x3(j,i,k)=x3(j,i,k)+qnr(j,i,k) + end if + end do + end do + end do + + case('w','W') + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = w(j,i,k) + x3(j,i,k)=x3(j,i,k)+w(j,i,k) + end do + end do + end do + end select + + + + end do + + do ic2=1,nc2d + + call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n + call stop2(9994) + end if + call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar' + call stop2(9995) + end if + + select case (trim(cvars2d(ic2))) + + case('ps','PS') + + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w2(j,i) = ps(j,i) + x2(j,i)=x2(j,i)+ps(j,i) + end do + end do + + case('sst','SST') + ! IGNORE SST IN HYBRID for now + + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w2(j,i) = zero + x2(j,i)=zero + end do + end do + + end select + end do + enddo + ! + ! CALCULATE ENSEMBLE MEAN + bar_norm = one/float(n_ens) + en_bar(m)%values=en_bar(m)%values*bar_norm + + ! Copy pbar to module array. ps_bar may be needed for vertical localization + ! in terms of scale heights/normalized p/p + do ic2=1,nc2d + + if(trim(cvars2d(ic2)) == 'ps'.or.trim(cvars2d(ic2)) == 'PS') then + + call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar to get ps_bar' + call stop2(9996) + end if + + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + ps_bar(j,i,1)=x2(j,i) + end do + end do + exit + end if + end do + + call mpi_barrier(mpi_comm_world,ierror) + ! + ! + ! CONVERT ENSEMBLE MEMBERS TO ENSEMBLE PERTURBATIONS + sig_norm=sqrt(one/max(one,n_ens-one)) + + do n=1,n_ens + do i=1,nelen + en_perts(n,m)%valuesr4(i)=(en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm + end do + end do + + enddo ! it 4d loop + ! CALCULATE ENSEMBLE SPREAD + write_ens_sprd=.true. + if(write_ens_sprd ) then + call this%ens_spread_dualres_regional(mype,en_perts,nelen) + call mpi_barrier(mpi_comm_world,ierror) ! do we need this mpi_barrier here? + endif + do m=1,ntlevs_ens call gsi_bundledestroy(en_bar(m),istatus) if(istatus/=0) then write(6,*)' in get_fv3_regional_ensperts_netcdf: trouble destroying en_bar bundle' - call stop2(9997) + call stop2(9997) endif - end do + end do - deallocate(en_bar) + deallocate(en_bar) ! - return + return 30 write(6,*) 'get_fv3_regional_ensperts_netcdf: open filelist failed ' call stop2(555) @@ -408,7 +523,8 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end subroutine get_fv3_regional_ensperts_run - subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz) + subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w) !$$$ subprogram documentation block ! first compied from general_read_arw_regional . . . . ! subprogram: general_read_fv3_regional read fv3sar model ensemble members @@ -430,65 +546,84 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g ! !$$$ end documentation block - use netcdf, only: nf90_nowrite - use netcdf, only: nf90_open,nf90_close - use netcdf, only: nf90_inq_dimid,nf90_inquire_dimension - use netcdf, only: nf90_inq_varid,nf90_inquire_variable,nf90_get_var - use kinds, only: r_kind,r_single,i_kind - use gridmod, only: eta1_ll,eta2_ll - use constants, only: zero,one,fv,zero_single,one_tenth,h300 - use hybrid_ensemble_parameters, only: grd_ens,q_hyb_ens - use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt - - use mpimod, only: mpi_comm_world,mpi_rtype - use netcdf_mod, only: nc_check - use gsi_rfv3io_mod,only: type_fv3regfilenameg - use gsi_rfv3io_mod,only:n2d - use gsi_rfv3io_mod,only:mype_t,mype_p ,mype_q,mype_oz - use constants, only: half,zero - use gsi_rfv3io_mod, only: gsi_fv3ncdf_read - use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_v1 - use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv - use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv_v1 - use gsi_rfv3io_mod, only: gsi_fv3ncdf2d_read_v1 - - implicit none - ! - ! Declare passed variables - class(get_fv3_regional_ensperts_class), intent(inout) :: this - type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) ::g_tsen, g_q,g_prsl - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig+1) ::g_prsi - ! - ! Declare local parameters - real(r_kind),parameter:: r0_01 = 0.01_r_kind - real(r_kind),parameter:: r10 = 10.0_r_kind - real(r_kind),parameter:: r100 = 100.0_r_kind - ! - ! Declare local variables - - integer(i_kind):: i,j,k,kp - integer(i_kind) iderivative - - - logical ice + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_open,nf90_close + use netcdf, only: nf90_inq_dimid,nf90_inquire_dimension + use netcdf, only: nf90_inq_varid,nf90_inquire_variable,nf90_get_var + use kinds, only: r_kind,r_single,i_kind + use gridmod, only: eta1_ll,eta2_ll + use constants, only: zero,one,fv,zero_single,one_tenth,h300 + use hybrid_ensemble_parameters, only: grd_ens,q_hyb_ens + use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt + + use mpimod, only: mpi_comm_world,mpi_rtype + use netcdf_mod, only: nc_check + use gsi_rfv3io_mod,only: type_fv3regfilenameg + use gsi_rfv3io_mod,only:n2d + use constants, only: half,zero + use gsi_rfv3io_mod, only: gsi_fv3ncdf_read + use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_v1 + use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv + use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv_v1 + use gsi_rfv3io_mod, only: gsi_fv3ncdf2d_read_v1 + use directDA_radaruse_mod, only: l_use_dbz_directDA + use gsi_bundlemod, only: gsi_gridcreate + use gsi_bundlemod, only: gsi_grid + use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundledestroy + use gsi_bundlemod, only: gsi_bundlegetvar + use hybrid_ensemble_parameters, only: grd_ens + use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt + use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval + + + + implicit none +! +! Declare passed variables + class(get_fv3_regional_ensperts_class), intent(inout) :: this + type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w + + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig+1) ::g_prsi + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) ::g_prsl ,g_tsen,g_q,g_delp +! +! Declare local parameters + real(r_kind),parameter:: r0_01 = 0.01_r_kind + real(r_kind),parameter:: r10 = 10.0_r_kind + real(r_kind),parameter:: r100 = 100.0_r_kind + ! +! Declare local variables + + integer(i_kind):: i,j,k,kp + + integer(i_kind) iderivative + + + logical ice + + character(len=24),parameter :: myname_ = 'general_read_fv3_regional' + type(gsi_bundle) :: gsibundle_fv3lam_ens_dynvar_nouv + type(gsi_bundle) :: gsibundle_fv3lam_ens_tracer_nouv + type(gsi_grid):: grid_ens - character(len=24),parameter :: myname_ = 'general_read_fv3_regional' + 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' + integer (i_kind) ier,istatus - 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' - - associate( this => this ) ! eliminates warning for unused dummy argument needed for binding - end associate + + associate( this => this ) ! eliminates warning for unused dummy argument needed for binding + end associate + call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) grid_spec=fv3_filenameginput%grid_spec ak_bk=fv3_filenameginput%ak_bk dynvars=fv3_filenameginput%dynvars @@ -496,46 +631,77 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g sfcdata=fv3_filenameginput%sfcdata couplerres=fv3_filenameginput%couplerres - -!cltthinktobe should be contained in variable like grd_ens + + + if (allocated(fv3lam_ens_io_dynmetvars2d_nouv) ) then + call gsi_bundlecreate(gsibundle_fv3lam_ens_dynvar_nouv,grid_ens,'gsibundle_fv3lam_ens_dynvar_nouv', istatus,& + names2d=fv3lam_ens_io_dynmetvars2d_nouv,names3d=fv3lam_ens_io_dynmetvars3d_nouv) + else + call gsi_bundlecreate(gsibundle_fv3lam_ens_dynvar_nouv,grid_ens,'gsibundle_fv3lam_ens_dynvar_nouv',istatus, & + names3d=fv3lam_ens_io_dynmetvars3d_nouv) + endif + + + if (allocated(fv3lam_ens_io_tracermetvars2d_nouv) ) then + call gsi_bundlecreate(gsibundle_fv3lam_ens_tracer_nouv,grid_ens,'gsibundle_fv3lam_ens_tracer_nouv', istatus,& + names2d=fv3lam_ens_io_tracermetvars2d_nouv,names3d=fv3lam_ens_io_tracermetvars3d_nouv) + else + call gsi_bundlecreate(gsibundle_fv3lam_ens_tracer_nouv,grid_ens,'gsibundle_fv3lam_ens_tracer_nouv',istatus, & + names3d=fv3lam_ens_io_tracermetvars3d_nouv) + endif + + + if(fv3sar_ensemble_opt == 0 ) then - call gsi_fv3ncdf_readuv(dynvars,g_u,g_v) + call gsi_fv3ncdf_readuv(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput) else - call gsi_fv3ncdf_readuv_v1(dynvars,g_u,g_v) + call gsi_fv3ncdf_readuv_v1(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput) endif if(fv3sar_ensemble_opt == 0) then - call gsi_fv3ncdf_read(dynvars,'T','t',g_tsen,mype_t) + call gsi_fv3ncdf_read(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& + fv3_filenameginput%dynvars,fv3_filenameginput) + call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& + fv3_filenameginput%tracers,fv3_filenameginput) else - call gsi_fv3ncdf_read_v1(dynvars,'t','T',g_tsen,mype_t) + call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& + fv3_filenameginput%dynvars,fv3_filenameginput) + call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& + fv3_filenameginput%tracers,fv3_filenameginput) endif + ier=0 + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q' ,g_q ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz' ,g_oz ,istatus );ier=ier+istatus + if (l_use_dbz_directDA) then + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'ql' ,g_ql ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qi' ,g_qi ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qr' ,g_qr ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qs' ,g_qs ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qg' ,g_qg ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qnr',g_qnr ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'w' , g_w ,istatus );ier=ier+istatus + end if + + if (fv3sar_ensemble_opt == 0) then - call gsi_fv3ncdf_read(dynvars,'DELP','delp',g_prsi,mype_p) + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'delp' ,g_delp ,istatus );ier=ier+istatus g_prsi(:,:,grd_ens%nsig+1)=eta1_ll(grd_ens%nsig+1) !thinkto be done , should use eta1_ll from ensemble grid do i=grd_ens%nsig,1,-1 - g_prsi(:,:,i)=g_prsi(:,:,i)*0.001_r_kind+g_prsi(:,:,i+1) + g_prsi(:,:,i)=g_delp(:,:,i)*0.001_r_kind+g_prsi(:,:,i+1) enddo - g_ps(:,:)=g_prsi(:,:,1) + g_ps(:,:)=g_prsi(:,:,1) else ! for the ensemble processed frm CHGRES - call gsi_fv3ncdf2d_read_v1(dynvars,'ps','PS',g_ps,mype_p) + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'ps' ,g_ps ,istatus );ier=ier+istatus g_ps=g_ps*0.001_r_kind do k=1,grd_ens%nsig+1 g_prsi(:,:,k)=eta1_ll(k)+eta2_ll(k)*g_ps enddo - endif - if(fv3sar_ensemble_opt == 0) then - call gsi_fv3ncdf_read(tracers,'SPHUM','sphum',g_q,mype_q) - call gsi_fv3ncdf_read(tracers,'O3MR','o3mr',g_oz,mype_oz) - else - call gsi_fv3ncdf_read_v1(tracers,'sphum','SPHUM',g_q,mype_q) - call gsi_fv3ncdf_read_v1(tracers,'o3mr','O3MR',g_oz,mype_oz) - endif - !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,grd_ens%nsig do j=1,grd_ens%lon2 @@ -544,136 +710,36 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g enddo enddo enddo - if (.not.q_hyb_ens) then - ice=.true. - iderivative=0 - do k=1,grd_ens%nsig - kp=k+1 - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - g_prsl(i,j,k)=(g_prsi(i,j,k)+g_prsi(i,j,kp))*half - end do - end do - end do - call genqsat(g_rh,g_tsen(1,1,1),g_prsl(1,1,1),grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative) - do k=1,grd_ens%nsig - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - g_rh(i,j,k) = g_q(i,j,k)/g_rh(i,j,k) - end do - end do + if (.not.q_hyb_ens) then + ice=.true. + iderivative=0 + do k=1,grd_ens%nsig + kp=k+1 + do j=1,grd_ens%lon2 + do i=1,grd_ens%lat2 + g_prsl(i,j,k)=(g_prsi(i,j,k)+g_prsi(i,j,kp))*half + end do + end do + end do + call genqsat(g_rh,g_tsen(1,1,1),g_prsl(1,1,1),grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative) + do k=1,grd_ens%nsig + do j=1,grd_ens%lon2 + do i=1,grd_ens%lat2 + g_rh(i,j,k) = g_q(i,j,k)/g_rh(i,j,k) + end do + end do + end do + else + do k=1,grd_ens%nsig + do j=1,grd_ens%lon2 + do i=1,grd_ens%lat2 + g_rh(i,j,k) = g_q(i,j,k) + end do end do - else - do k=1,grd_ens%nsig - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - g_rh(i,j,k) = g_q(i,j,k) - end do - end do - end do - end if - - - - - - return - end subroutine general_read_fv3_regional - - subroutine general_read_fv3_regional_dirZDA(this,fv3_filenameginput,g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w) - !$$$ subprogram documentation block - ! firstly copied from general_read_arw_regional . . . . - ! subprogram: general_read_fv3_regional_dirZda read fv3sar model ensemble members - ! prgmmr: Ting org: emc/ncep date: 2018 - ! - ! abstract: read ensemble members from the fv3 model netcdf format (dynamic and tracer) - ! for use with hybrid ensemble option. - ! - ! program history log: - ! 2018- Ting - intial versions - ! 2019-03-13 CAPS(C. Tong) - Port direct radar DA capabilities - ! 2021-05-05 CAPS(J. Park) - Modified to read hydrometeors and w only - ! - ! input argument list: - ! - ! output argument list: - ! - ! attributes: - ! language: f90 - ! machine: ibm RS/6000 SP - ! - !$$$ end documentation block - - use netcdf, only: nf90_nowrite - use netcdf, only: nf90_open,nf90_close - use netcdf, only: nf90_inq_dimid,nf90_inquire_dimension - use netcdf, only: nf90_inq_varid,nf90_inquire_variable,nf90_get_var - use kinds, only: r_kind,r_single,i_kind - use constants, only: one - use hybrid_ensemble_parameters, only: grd_ens - use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt - - use mpimod, only: mpi_comm_world,mpi_rtype - use netcdf_mod, only: nc_check - use gsi_rfv3io_mod,only: type_fv3regfilenameg - use gsi_rfv3io_mod,only:n2d - use gsi_rfv3io_mod,only:mype_qr,mype_qs,mype_qg,mype_qnr,mype_w - use gsi_rfv3io_mod,only:mype_ql,mype_qi - use gsi_rfv3io_mod, only: gsi_fv3ncdf_read - use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_v1 - use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt - use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval - - implicit none - ! - ! Declare passed variables - class(get_fv3_regional_ensperts_class), intent(inout) :: this - type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_ql,g_qi,g_qr - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_qs,g_qg,g_qnr,g_w - ! - ! Declare local parameters - real(r_kind),parameter:: r0_01 = 0.01_r_kind - real(r_kind),parameter:: r10 = 10.0_r_kind - real(r_kind),parameter:: r100 = 100.0_r_kind - ! - ! Declare local variables - - integer(i_kind):: i,j,k - - character(len=24),parameter :: myname_ = 'general_read_fv3_regional_dirZDA' - - 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' - - associate( this => this ) ! eliminates warning for unused dummy argument needed for binding - end associate + end do + end if - grid_spec=fv3_filenameginput%grid_spec - ak_bk=fv3_filenameginput%ak_bk - dynvars=fv3_filenameginput%dynvars - tracers=fv3_filenameginput%tracers - sfcdata=fv3_filenameginput%sfcdata - couplerres=fv3_filenameginput%couplerres - if(fv3sar_ensemble_opt == 0) then -! variables for direct reflectivity DA - call gsi_fv3ncdf_read(dynvars,'W','w',g_w,mype_w) - call gsi_fv3ncdf_read(tracers,'LIQ_WAT','liq_wat',g_ql,mype_ql) - call gsi_fv3ncdf_read(tracers,'ICE_WAT','ice_wat',g_qi,mype_qi) - call gsi_fv3ncdf_read(tracers,'RAINWAT','rainwat',g_qr,mype_qr) - call gsi_fv3ncdf_read(tracers,'SNOWWAT','snowwat',g_qs,mype_qs) - call gsi_fv3ncdf_read(tracers,'GRAUPEL','graupel',g_qg,mype_qg) - call gsi_fv3ncdf_read(tracers,'RAIN_NC','rain_nc',g_qnr,mype_qnr) - else - write(6,*) "get_fv3_regional_ensperts for 'fv3sar_bg_opt == 0 is only available & - for now in direct relfectivity DA" - stop - end if ! CV transform do k=1,grd_ens%nsig do i=1,grd_ens%lon2 @@ -717,10 +783,12 @@ subroutine general_read_fv3_regional_dirZDA(this,fv3_filenameginput,g_ql,g_qi,g_ enddo enddo enddo + call gsi_bundledestroy(gsibundle_fv3lam_ens_dynvar_nouv) + call gsi_bundledestroy(gsibundle_fv3lam_ens_tracer_nouv) return - end subroutine general_read_fv3_regional_dirZDA + end subroutine general_read_fv3_regional subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) !$$$ subprogram documentation block diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 44a5344297..a61433d653 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -13,6 +13,7 @@ module gsi_rfv3io_mod ! 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 ! subroutines included: ! sub gsi_rfv3io_get_grid_specs ! sub read_fv3_files @@ -37,6 +38,9 @@ module gsi_rfv3io_mod 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 implicit none public type_fv3regfilenameg public bg_fv3regfilenameg @@ -44,7 +48,7 @@ module gsi_rfv3io_mod ! directory names (hardwired for now) type type_fv3regfilenameg - character(len=:),allocatable :: grid_spec !='fv3_grid_spec' + 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' @@ -55,6 +59,7 @@ module gsi_rfv3io_mod 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(:,:) @@ -62,6 +67,24 @@ module gsi_rfv3io_mod 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), parameter :: vardynvars(ndynvarslist) =(/"u","v","u_w","u_s", & + "v_w","v_s","t","tv","tsen","w","delp","ps","delzinc"/) + character(len=max_varname_length), parameter :: vartracers(ntracerslist) =(/'q','oz', & + 'ql','qi','qr','qs','qg','qnr'/) + character(len=max_varname_length), parameter :: varfv3name(15) =(/'u','v','W','T','delp','sphum','o3mr', & + 'liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ps','DZ'/) + character(len=max_varname_length), parameter :: vgsiname(15) =(/'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 @@ -80,6 +103,8 @@ module gsi_rfv3io_mod 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 @@ -100,6 +125,20 @@ module gsi_rfv3io_mod 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, & @@ -144,8 +183,7 @@ subroutine fv3regfilename_init(this,grid_spec_input,ak_bk_input,dynvars_input, & else this%couplerres='coupler.res' endif - - end subroutine fv3regfilename_init +end subroutine fv3regfilename_init subroutine gsi_rfv3io_get_grid_specs(fv3filenamegin,ierr) @@ -200,14 +238,14 @@ subroutine gsi_rfv3io_get_grid_specs(fv3filenamegin,ierr) 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 + 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 - coupler_res_filenam=fv3filenamegin%couplerres - grid_spec=fv3filenamegin%grid_spec - ak_bk=fv3filenamegin%ak_bk + 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') @@ -289,7 +327,7 @@ subroutine gsi_rfv3io_get_grid_specs(fv3filenamegin,ierr) 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 - elseif(grid_type_fv3_regional == 2) then + 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 @@ -622,16 +660,24 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) !$$$ end documentation block use kinds, only: r_kind,i_kind use mpimod, only: npe - use guess_grids, only: ges_tsen,ges_prsi + 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_bundlegetpointer + 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 @@ -640,13 +686,15 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) 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_ql=>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() @@ -656,47 +704,24 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) 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=:),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=mod(1,npe) - mype_t=mod(2,npe) - mype_p=mod(3,npe) - mype_q=mod(4,npe) - mype_ql=mod(5,npe) - mype_oz=mod(6,npe) - mype_2d=mod(7,npe) - mype_delz=mod(8,npe) - - if (l_use_dbz_directDA) then ! direct reflectivity DA - if(npe< 15) then - call die('read_fv3_netcdf_guess','not enough PEs to read in fv3 fields for direct reflectivity DA' ) - endif - mype_qi=mod(9,npe) - mype_qr=mod(10,npe) - mype_qs=mod(11,npe) - mype_qg=mod(12,npe) - mype_qnr=mod(13,npe) - mype_w=mod(14,npe) - end if - + 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 + ijns2d(i)=ijn_s(i)*n2d enddo displss(1)=0 displsz_g(1)=0 @@ -707,102 +732,314 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) 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 (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), 'qr' ,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(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 + it=ntguessig - 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) - ges_ps=ges_ps*0.001_r_kind - do k=1,nsig - ges_prsi(:,:,k,it)=eta1_ll(k)+eta2_ll(k)*ges_ps - enddo - endif + 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) + - 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) + 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 (l_use_dbz_directDA ) then - if( fv3sar_bg_opt == 0) then - call gsi_fv3ncdf_read(dynvars,'W','w',ges_w,mype_w) - call gsi_fv3ncdf_read(tracers,'LIQ_WAT','liq_wat',ges_ql,mype_ql) - call gsi_fv3ncdf_read(tracers,'ICE_WAT','ice_wat',ges_qi,mype_qi) - call gsi_fv3ncdf_read(tracers,'RAINWAT','rainwat',ges_qr,mype_qr) - ges_iqr=ges_qr - call gsi_fv3ncdf_read(tracers,'SNOWWAT','snowwat',ges_qs,mype_qs) - call gsi_fv3ncdf_read(tracers,'GRAUPEL','graupel',ges_qg,mype_qg) - call gsi_fv3ncdf_read(tracers,'RAIN_NC','rain_nc',ges_qnr,mype_qnr) + 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 - write(6,*) "FV3 IO READ for 'fv3sar_bg_opt == 0' is only available for now in direct reflectivity DA" - stop - end if + 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 - 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 + + 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 @@ -844,7 +1081,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) 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 + 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 @@ -865,52 +1102,109 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) 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 + 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 ) + 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)) + else + write(*,*) "wrong dimension number ndim =",ndim + call stop2(119) + endif + + call fv3_h_to_ll(sfc(:,:,1),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 + iret=nf90_close(gfile_loc) + + !!!! read in orog from dynam !!!!!!!!!!!! + 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 - 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,grid_reverse_flag) + + 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,grid_reverse_flag) kk=0 do n=1,npe @@ -923,57 +1217,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z) 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,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) - endif ! mype + deallocate (dim_id,sfc,sfc1,dim) + endif ! mype !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1051,7 +1297,7 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io) allocate (work(itotsub)) if(mype==mype_io ) then - iret=nf90_open(trim(filenamein),nf90_nowrite,gfile_loc) + 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) @@ -1096,14 +1342,14 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io) return end subroutine gsi_fv3ncdf2d_read_v1 -subroutine gsi_fv3ncdf_read(filenamein,varname,varname2,work_sub,mype_io) +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 -! -! abstract: read in a field from a netcdf FV3 file in mype_io -! then scatter the field to each PE +! 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: @@ -1123,95 +1369,79 @@ subroutine gsi_fv3ncdf_read(filenamein,varname,varname2,work_sub,mype_io) 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 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 general_commvars_mod, only: ltosi_s,ltosj_s + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub 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 + 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 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,grid_reverse_flag) - 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) + 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: 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) + 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/) + + iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) + 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 ! ilevtot + call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) + iret=nf90_close(gfile_loc) + deallocate (uu2d) - deallocate (work) + + return end subroutine gsi_fv3ncdf_read -subroutine gsi_fv3ncdf_read_v1(filenamein,varname,varname2,work_sub,mype_io) +subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . @@ -1240,98 +1470,86 @@ subroutine gsi_fv3ncdf_read_v1(filenamein,varname,varname2,work_sub,mype_io) 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 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 general_commvars_mod, only: ltosi_s,ltosj_s + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub 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 + 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) nztmp,nzp1,kk,j,mm1,i,ir,ii,jj - integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid + integer(i_kind) nzp1,mm1 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 + 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 - 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 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) - allocate(uu(nx,ny,nsig)) - allocate(temp0(nx,ny,nsig+1)) + 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) - 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,grid_reverse_flag) - 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 (uu2d) - 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) +subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_readuv @@ -1353,68 +1571,80 @@ subroutine gsi_fv3ncdf_readuv(dynvarsfile,ges_u,ges_v) ! !$$$ 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 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,nya,nxa,fv3uv2earth - use general_commvars_mod, only: ltosi_s,ltosj_s + use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub 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 + 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) 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 + integer(i_kind) nz,nzp1,mm1 - 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)) + 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 + 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 + 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/) + + 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) + 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) - 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) ! 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. @@ -1430,60 +1660,24 @@ subroutine gsi_fv3ncdf_readuv(dynvarsfile,ges_u,ges_v) ! ! and the last input parameter for fv3_h_to_ll is alway true: ! -! call fv3_h_to_ll(u,a,nx,ny,nxa,nya,.true.) ! - if(trim(name)=='u'.or.trim(name)=='U') then - allocate(temp1(dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3)))) - if(.not.grid_reverse_flag) call reverse_grid_r_uv(uu,dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3))) - temp1=uu - else if(trim(name)=='v'.or.trim(name)=='V') then - if(.not.grid_reverse_flag) call reverse_grid_r_uv(uu,dim(dim_id(1)),dim(dim_id(2)),dim(dim_id(3))) - 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,.true.) - else - call fv3_h_to_ll(v,a,nx,ny,nxa,nya,.true.) - 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) + 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 + + call general_grid2sub(grd_uv,hwork,worksub) + ges_u=worksub(1,:,:,:) + ges_v=worksub(2,:,:,:) + iret=nf90_close(gfile_loc) + + deallocate(u2d,v2d,uc2d,vc2d,worksub) end subroutine gsi_fv3ncdf_readuv -subroutine gsi_fv3ncdf_readuv_v1(dynvarsfile,ges_u,ges_v) +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 . . . +! 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: @@ -1501,113 +1695,96 @@ subroutine gsi_fv3ncdf_readuv_v1(dynvarsfile,ges_u,ges_v) !$$$ 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 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,nya,nxa,fv3uv2earth - use general_commvars_mod, only: ltosi_s,ltosj_s + use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub 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 + 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) nztmp,nzp1,kk,j,mm1,i,ir,ii,jj - integer(i_kind) ndimensions,nvariables,nattributes,unlimiteddimid + 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 (work(itotsub*nsig)) + allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%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) + 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/) - 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,grid_reverse_flag) - 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,grid_reverse_flag) - 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) + 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) @@ -1635,16 +1812,17 @@ subroutine wrfv3_netcdf(fv3filenamegin) 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 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 gridmod, only: lat2,lon2,nsig 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 @@ -1652,8 +1830,6 @@ subroutine wrfv3_netcdf(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() @@ -1661,8 +1837,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_v =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL() - real(r_kind),allocatable,dimension(:,:,:)::ges_delzinc - integer(i_kind) k + integer(i_kind) i,k real(r_kind),pointer,dimension(:,:,:):: ges_ql =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qi =>NULL() @@ -1671,12 +1846,14 @@ subroutine wrfv3_netcdf(fv3filenamegin) 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(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 - dynvars=fv3filenamegin%dynvars - tracers=fv3filenamegin%tracers it=ntguessig ier=0 @@ -1693,36 +1870,18 @@ subroutine wrfv3_netcdf(fv3filenamegin) 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('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) - 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_fv3ncdf_write_fv3_dz(dynvars,"DZ",ges_delzinc,mype_delz,add_saved) - deallocate(ges_delzinc) - endif - - else - call gsi_fv3ncdf_write_v1(dynvars,'t',ges_tsen(1,1,1,it),mype_t,add_saved) - call gsi_fv3ncdf_write_v1(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) - + 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 @@ -1731,26 +1890,65 @@ subroutine wrfv3_netcdf(fv3filenamegin) 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 + ges_ps=ges_ps*1000.0_r_kind + call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'ps',ges_ps,istatus) + endif ! write out - if( fv3sar_bg_opt == 0) then - add_saved=.true. - call gsi_fv3ncdf_write(dynvars,'W',ges_w,mype_w,add_saved) - call gsi_fv3ncdf_write(tracers,'liq_wat',ges_ql,mype_ql,add_saved) - call gsi_fv3ncdf_write(tracers,'ice_wat',ges_qi,mype_qi,add_saved) - call gsi_fv3ncdf_write(tracers,'rainwat',io_arr_qr,mype_qr,add_saved) - call gsi_fv3ncdf_write(tracers,'snowwat',io_arr_qs,mype_qs,add_saved) - call gsi_fv3ncdf_write(tracers,'graupel',io_arr_qg,mype_qg,add_saved) - call gsi_fv3ncdf_write(tracers,'rain_nc',io_arr_qnr,mype_qnr,add_saved) - else - write(6,*) "FV3 IO Write for 'fv3sar_bg_opt == 0 is only available for now in direct relfectivity DA" - stop - end if - end if + 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) + +! additional I/O for direct reflectivity DA capabilities end subroutine wrfv3_netcdf -subroutine gsi_fv3ncdf_writeuv(dynvars,varu,varv,mype_io,add_saved) +subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeuv @@ -1775,254 +1973,125 @@ subroutine gsi_fv3ncdf_writeuv(dynvars,varu,varv,mype_io,add_saved) ! !$$$ 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 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 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 + use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none - character(len=*),intent(in) :: dynvars !='fv3_dynvars' + 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 - 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 + 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,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 + 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 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) + + 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 - 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)) + end do + call general_sub2grid(grd_uv,worksub,hwork) + 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 + 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 !!!!!!!!!!!!!!!! - call check( nf90_get_var(gfile_loc,ugrd_VarId,work_bu) ) - call check( nf90_get_var(gfile_loc,vgrd_VarId,work_bv) ) - if(.not.grid_reverse_flag) then - call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,nsig) - call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,nsig) - endif - 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,.true.) - call fv3_h_to_ll(v,workav2,nlon_regional,nlat_regional,nlon,nlat,.true.) + 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) ) + 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(:,:,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) + 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(:,:,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 - if(.not.grid_reverse_flag) then - call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,nsig) - call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,nsig) - 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 + 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 + + 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) ) + enddo !ilevltot + call check( nf90_close(gfile_loc) ) + deallocate(work_bu,work_bv,u2d,v2d) + deallocate(work_au,work_av) - 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) ) - allocate( workb2(nlon_regional,nlat_regional)) - - if(add_saved)then - 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,grid_reverse_flag) -!!!!!!! analysis_inc Psfc: work_a - work_a(:,:)=work_a(:,:)-worka2(:,:) - call fv3_ll_to_h(work_a,workb2,nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) - 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,grid_reverse_flag) - 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) ) - if (allocated(worka2)) deallocate(worka2) - if (allocated(workb2)) deallocate(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) +subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_writeuv @@ -2050,305 +2119,196 @@ subroutine gsi_fv3ncdf_writeuv_v1(dynvars,varu,varv,mype_io,add_saved) !$$$ 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 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 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 - + use general_sub2grid_mod, only: sub2grid_info,general_sub2grid 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 + 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,n,k,ns,kr,m - integer(i_kind) ilev0 - 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 + 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 - ilev0=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) + 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 - 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) -! u and v would contain winds at either D-grid or A-grid -! 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) ) + 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/) - allocate( u(nlon_regional,nlat_regional)) - allocate( v(nlon_regional,nlat_regional)) - allocate( work_bu_s(nlon_regional,nlat_regional+1,nsig+1)) - allocate( work_bv_s(nlon_regional,nlat_regional+1,nsig+1)) - allocate( work_bu_w(nlon_regional+1,nlat_regional,nsig+1)) - allocate( work_bv_w(nlon_regional+1,nlat_regional,nsig+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) ) + 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) ) - 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) ) - - if(add_saved)then - allocate( workau2(nlat,nlon),workav2(nlat,nlon)) - do k=1,nsig - do j=1,nlat_regional - u(:,j)=half * (work_bu_s(:,j,ilev0+k)+ work_bu_s(:,j+1,ilev0+k)) - enddo - do i=1,nlon_regional - v(i,:)=half*(work_bv_w(i,:,ilev0+k)+work_bv_w(i+1,:,ilev0+k)) - enddo - call fv3_h_to_ll(u,workau2,nlon_regional,nlat_regional,nlon,nlat,grid_reverse_flag) - call fv3_h_to_ll(v,workav2,nlon_regional,nlat_regional,nlon,nlat,grid_reverse_flag) + 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(:,:,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,grid_reverse_flag) - call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) + 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*(u(i-1,:)+u(i,:)) - workbv_w2(i,:)=half*(v(i-1,:)+v(i,:)) - 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 - workbu_s2(:,j)=half*(u(:,j-1)+u(:,j)) - workbv_s2(:,j)=half*(v(:,j-1)+v(:,j)) - 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(:,:,ilev0+k)=work_bu_w(:,:,ilev0+k)+workbu_w2(:,:) - work_bu_s(:,:,ilev0+k)=work_bu_s(:,:,ilev0+k)+workbu_s2(:,:) - work_bv_w(:,:,ilev0+k)=work_bv_w(:,:,ilev0+k)+workbv_w2(:,:) - work_bv_s(:,:,ilev0+k)=work_bv_s(:,:,ilev0+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,grid_reverse_flag) - call fv3_ll_to_h(work_av(:,:,k),v,nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) - - do i=2,nlon_regional - work_bu_w(i,:,k)=half*(u(i-1,:)+u(i,:)) - work_bv_w(i,:,k)=half*(v(i-1,:)+v(i,:)) - enddo - work_bu_w(1,:,ilev0+k)=u(1,:) - work_bv_w(1,:,ilev0+k)=v(1,:) - work_bu_w(nlon_regional+1,:,ilev0+k)=u(nlon_regional,:) - work_bv_w(nlon_regional+1,:,ilev0+k)=v(nlon_regional,:) - - do j=2,nlat_regional - work_bu_s(:,j,ilev0+k)=half*(u(:,j-1)+u(:,j)) - work_bv_s(:,j,ilev0+k)=half*(v(:,j-1)+v(:,j)) - enddo - work_bu_s(:,1,ilev0+k)=u(:,1) - work_bv_s(:,1,ilev0+k)=v(:,1) - work_bu_s(:,nlat_regional+1,ilev0+k)=u(:,nlat_regional) - work_bv_s(:,nlat_regional+1,ilev0+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 + 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(work))deallocate(work) - if(allocated(work_sub))deallocate(work_sub) + if(allocated(worksub))deallocate(worksub) 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,grid_reverse_flag) -!!!!!!! analysis_inc Psfc: work_a - work_a(:,:)=work_a(:,:)-worka2(:,:) - call fv3_ll_to_h(work_a,workb2,nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) - work_b(:,:)=work_b(:,:)+workb2(:,:) - else - call fv3_ll_to_h(work_a,work_b,nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) - - endif - - call check( nf90_put_var(gfile_loc,VarId,work_b) ) - call check( nf90_close(gfile_loc) ) - if (allocated(worka2)) deallocate(worka2) - if ( allocated(workb2)) deallocate(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) +subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write @@ -2372,193 +2332,98 @@ subroutine gsi_fv3ncdf_write(filename,varname,var,mype_io,add_saved) ! !$$$ 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 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 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 + 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 - 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 + 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) 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 + 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 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)) + call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) ) + 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 + call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + 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 + call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + - 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,grid_reverse_flag) + call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! 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,grid_reverse_flag) - 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,grid_reverse_flag) - enddo - endif + 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 + call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) - 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 + enddo !ilevtotl loop + call check(nf90_close(gfile_loc)) + deallocate(work_b,work_a) - deallocate(work,work_sub) end subroutine gsi_fv3ncdf_write -subroutine gsi_fv3ncdf_write_fv3_dz(filename,varname,varinc,mype_io,add_saved) -!$$$ subprogram documentation block -! . . . . -! subprogram: gsi_nemsio_write_fv3_dz from gsi_nemsio_write -! pgrmmr: lei -! -! abstract: -! -! program history log: -! -! input argument list: -! varin -! 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 ) :: varinc(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_ainc - real(r_kind),allocatable,dimension(:,:,:):: work_b - real(r_kind),allocatable,dimension(:,:):: workb2 - - 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)=varinc(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_ainc(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_ainc(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(.not. add_saved)then - write(6,*)'here the input is increments to be added to the read-in background, & - hence, add_saved has to be true' - endif - allocate( workb2(nlon_regional,nlat_regional)) - call check( nf90_get_var(gfile_loc,VarId,work_b) ) - - do k=1,nsig - call fv3_ll_to_h(work_ainc(1,1,k),workb2(:,:),nlon,nlat,nlon_regional,nlat_regional,.true.) -!!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! - work_b(:,:,k)=workb2(:,:)+work_b(:,:,k) - enddo - deallocate(workb2) - - 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_ainc) - end if !mype_io - - deallocate(work,work_sub) - -end subroutine gsi_fv3ncdf_write_fv3_dz subroutine check(status) use kinds, only: i_kind use netcdf, only: nf90_noerr,nf90_strerror @@ -2569,7 +2434,7 @@ subroutine check(status) stop end if end subroutine check -subroutine gsi_fv3ncdf_write_v1(filename,varname,var,mype_io,add_saved) +subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3filenamegin) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write @@ -2593,91 +2458,90 @@ subroutine gsi_fv3ncdf_write_v1(filename,varname,var,mype_io,add_saved) ! !$$$ 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 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 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 + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none - real(r_kind) ,intent(in ) :: var(lat2,lon2,nsig) - integer(i_kind),intent(in ) :: mype_io + type(sub2grid_info), intent(in):: grd_ionouv + type(gsi_bundle),intent(inout) :: cstate_nouv logical ,intent(in ) :: add_saved - character(*) ,intent(in ) :: varname,filename + 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) :: ilev0 - 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 + 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 - - 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+1)) - - call check( nf90_open(trim(filename),nf90_write,gfile_loc) ) - call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) - call check( nf90_get_var(gfile_loc,VarId,work_b) ) - ilev0=1 - - if(add_saved)then - allocate( workb2(nlon_regional,nlat_regional)) - allocate( worka2(nlat,nlon)) + 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 - do k=1,nsig - call fv3_h_to_ll(work_b(:,:,ilev0+k),worka2,nlon_regional,nlat_regional,nlon,nlat,grid_reverse_flag) + call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! 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,grid_reverse_flag) - work_b(:,:,ilev0+k)=work_b(:,:,ilev0+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,ilev0+k),nlon,nlat,nlon_regional,nlat_regional,grid_reverse_flag) - 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) + 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 @@ -3163,25 +3027,25 @@ subroutine convert_cvpnx_to_nx(qnx_arr,cvpnr,cvpnr_pvalue,cloud_nt_updt,q_arr,qr ! 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) + 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 + 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 if + tmparr_qnr(j,i,k)=qnr_tmp end do end do @@ -3190,5 +3054,99 @@ subroutine convert_cvpnx_to_nx(qnx_arr,cvpnr,cvpnr_pvalue,cloud_nt_updt,q_arr,qr 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 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 1e12efdd39..1c8b73841f 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -1678,6 +1678,10 @@ subroutine gsimain_initialize end if end if if (fv3sar_bg_opt /= 0) l_reg_update_hydro_delz=.false. + if(regional_ensemble_option == 5 .and. (fv3sar_ensemble_opt /= fv3sar_bg_opt)) then + write(6,*)'this setup doesn"t work, stop' + call stop2(137) + endif if (anisotropic) then call init_fgrid2agrid(pf2aP1) call init_fgrid2agrid(pf2aP2) diff --git a/src/gsi/mod_fv3_lola.f90 b/src/gsi/mod_fv3_lola.f90 index f1726fdde8..84f8144968 100644 --- a/src/gsi/mod_fv3_lola.f90 +++ b/src/gsi/mod_fv3_lola.f90 @@ -361,6 +361,13 @@ subroutine generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt) endif jb2=jb1 jb1=gya + if(ib1+1 > nx)then !this block( 6 lines) is copied from GSL gsi repository + ib1=ib1-1 + endif + if(jb1+1 > ny)then + jb1=jb1-1 + endif + if((ib1 == ib2) .and. (jb1 == jb2)) exit if(n==3 ) then