From 4e8107c0ae054bc9dcc7b9a2684479d97a1e4261 Mon Sep 17 00:00:00 2001 From: "Tseganeh Z. Gichamo" Date: Tue, 19 Mar 2024 12:44:32 -0600 Subject: [PATCH] Add parallel netcdf read/write from EnKF for sfc files (paranc option) (#709) --- src/enkf/controlvec.f90 | 16 +- src/enkf/gridio_gfs.f90 | 453 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 452 insertions(+), 17 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 808eae2e28..0961549634 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -191,7 +191,7 @@ subroutine read_control() ! read ensemble members on IO tasks implicit none real(r_double) :: t1,t2 -integer(i_kind) :: nb,ne +integer(i_kind) :: nb,nlev,ne integer(i_kind) :: q_ind integer(i_kind) :: ierr @@ -218,19 +218,23 @@ subroutine read_control() if (nproc == 0) t1 = mpi_wtime() call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in readgrid_pnc on root',t2-t1,'secs' + end if end if if (nproc <= ntasks_io-1) then if (.not. paranc) then if (nproc == 0) t1 = mpi_wtime() call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in readgrid on root',t2-t1,'secs' + end if end if !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat) q_ind = getindex(cvars3d, 'q') - if (nproc == 0) then - t2 = mpi_wtime() - print *,'time in readgridata on root',t2-t1,'secs' - end if if (pseudo_rh .and. q_ind > 0) then do ne=1,nanals_per_iotask do nb=1,nbackgrounds @@ -357,7 +361,7 @@ subroutine write_control(no_inflate_flag) endif deallocate(grdin_mean) t2 = mpi_wtime() - print *,'time in write_control on root',t2-t1,'secs' + print *,'time in write_control paranc on root',t2-t1,'secs' endif end if diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 8afb012fcf..35a0c3fbe4 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -89,15 +89,21 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & real(r_kind), dimension(ndimspec) :: vrtspec,divspec real(r_kind), allocatable, dimension(:) :: psg,pstend,ak,bk real(r_single),allocatable,dimension(:,:,:) :: ug3d,vg3d - type(Dataset) :: dset + type(Dataset) :: dset, dset_sfc type(Dimension) :: londim,latdim,levdim integer(i_kind) :: u_ind, v_ind, tv_ind, q_ind, oz_ind, cw_ind integer(i_kind) :: qr_ind, qs_ind, qg_ind integer(i_kind) :: tsen_ind, ql_ind, qi_ind, prse_ind integer(i_kind) :: ps_ind, pst_ind, sst_ind + ! surface + integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind, soilt3_ind + integer(i_kind) :: soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind integer(i_kind) :: k,iret,nb,i,imem,idvc,nlonsin,nlatsin,nlevsin,ne,nanal + ! surface + integer(i_kind) :: nlonsin_sfc,nlatsin_sfc + logical ice logical use_full_hydro integer(i_kind), allocatable, dimension(:) :: mem_pe, lev_pe1, lev_pe2, iocomms @@ -111,12 +117,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, read_sfc_file, read_atm_file) - if (read_sfc_file) then - print *,'paranc not supported for reading surface files' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - endif - ! figure out what member to read and do MPI sub-communicator things allocate(mem_pe(0:numproc-1)) allocate(iocomms(nanals)) @@ -152,6 +152,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & displs(i+1) = ((lev_pe1(i)-1)*nlons*nlats) end do + if (read_atm_file) then ! loop through times and do the read ne = 1 @@ -159,7 +160,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & write(charnanal,'(a3, i3.3)') 'mem', nanal filename = trim(adjustl(datapath))//trim(adjustl(fileprefixes(nb)))//trim(charnanal) - sfcfilename = trim(adjustl(datapath))//trim(adjustl(filesfcprefixes(nb)))//trim(charnanal) if (use_gfs_ncio) then dset = open_dataset(filename, paropen=.true., mpicomm=iocomms(mem_pe(nproc))) londim = get_dim(dset,'grid_xt'); nlonsin = londim%len @@ -496,6 +496,141 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end do backgroundloop ! loop over backgrounds to read in + end if !read_atm_file + + if (read_sfc_file) then + ! loop through times and do the read + ne = 1 + sfcbackgroundloop: do nb=1,ntimes + + write(charnanal,'(a3, i3.3)') 'mem', nanal + sfcfilename = trim(adjustl(datapath))//trim(adjustl(filesfcprefixes(nb)))//trim(charnanal) + if (use_gfs_ncio) then + dset_sfc = open_dataset(sfcfilename, paropen=.true., mpicomm=iocomms(mem_pe(nproc))) + else + write(6,*)'READGRIDDATA_PNC sfc: ***FATAL ERROR*** parallel read only supported for netCDF' , ' PROGRAM STOPS' + call mpi_barrier(mpi_comm_world,ierr) + call mpi_finalize(ierr) + end if + if ( reducedgrid ) then + write(6,*) "READGRIDDATA_PNC sfc: reducedgrid=T interpolation not valid for writing sfc files" + call mpi_barrier(mpi_comm_world,ierr) + call mpi_finalize(ierr) + endif + + ! land sfc DA variables + tmp2m_ind = getindex(vars2d, 't2m') + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'st1') + slc1_ind = getindex(vars2d, 'sl1') + soilt2_ind = getindex(vars2d, 'st2') + slc2_ind = getindex(vars2d, 'sl2') + soilt3_ind = getindex(vars2d, 'st3') + slc3_ind = getindex(vars2d, 'sl3') + soilt4_ind = getindex(vars2d, 'st4') + slc4_ind = getindex(vars2d, 'sl4') + + ! read in sfc vars, if requested + if (tmp2m_ind > 0) then + call read_vardata(dset_sfc, 'tmp2m', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading tmp2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + tmp2m_ind,nb,ne)) + endif + if (spfh2m_ind > 0) then + call read_vardata(dset_sfc, 'spfh2m', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading spfh2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + spfh2m_ind,nb,ne)) + endif + if (soilt1_ind > 0) then + call read_vardata(dset_sfc, 'soilt1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt1_ind,nb,ne)) + endif + if (soilt2_ind > 0) then + call read_vardata(dset_sfc, 'soilt2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt2_ind,nb,ne)) + endif + if (soilt3_ind > 0) then + call read_vardata(dset_sfc, 'soilt3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt3_ind,nb,ne)) + endif + if (soilt4_ind > 0) then + call read_vardata(dset_sfc, 'soilt4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne)) + endif + if (slc1_ind > 0) then + call read_vardata(dset_sfc, 'soill1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne)) + endif + if (slc2_ind > 0) then + call read_vardata(dset_sfc, 'soill2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne)) + endif + if (slc3_ind > 0) then + call read_vardata(dset_sfc, 'soill3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne)) + endif + if (slc4_ind > 0) then + call read_vardata(dset_sfc, 'soill4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill4' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc4_ind,nb,ne)) + endif + + ! bring all the subdomains back to the main PE + call mpi_barrier(iocomms(mem_pe(nproc)), iret) + if (allocated(values_2d)) deallocate(values_2d) + call close_dataset(dset_sfc) + call mpi_barrier(iocomms(mem_pe(nproc)), iret) + + end do sfcbackgroundloop ! loop over backgrounds to read in + end if !if (read_sfc_file) + ! remove the sub communicators call mpi_barrier(iocomms(mem_pe(nproc)), iret) call mpi_comm_free(iocomms(mem_pe(nproc)), iret) @@ -1322,6 +1457,20 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ integer(i_kind) :: ncstart(4), nccount(4) logical :: nocompress + logical :: write_sfc_file, write_atm_file + character(len=max_varname_length), dimension(n3d) :: no_vars3d + character(len=max_varname_length), dimension(n2d) :: no_vars2d + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + + if (write_sfc_file ) then + ! adding the sfc increments requires adjusting several other variables. + ! This is done is a separate program. + if (nproc == 0) write(6,*) 'gridio/writegriddata_pnc: not coded to write sfc analysis, will write increment for sfc fields' + no_vars3d='' + call writeincrement_pnc(no_vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) + endif + nocompress = .true. if (nccompress) nocompress = .false. @@ -3616,6 +3765,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! soil / snow mask (not fixed) integer(i_kind), dimension(nlons,nlats) :: mask logical :: write_sfc_file, write_atm_file + real(r_double) :: t1,t2 call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) @@ -3627,6 +3777,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, nccount = (/nlons, nlats, nlevs/) if ( write_atm_file) then + if (nproc == 0) t1 = mpi_wtime() ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -3954,10 +4105,15 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, end do backgroundloop ! loop over backgrounds to read in end do ensmemloop ! loop over ens members to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement atm_file on root',t2-t1,'secs' + endif endif ! write_atm_file if (write_sfc_file) then + if (nproc == 0) t1 = mpi_wtime() ne = 0 sfcensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -4199,6 +4355,10 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, end do sfcbackgroundloop ! loop over backgrounds to read in end do sfcensmemloop ! loop over ens members to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement sfc_file on root',t2-t1,'secs' + endif endif ! write_sfc_file @@ -4225,7 +4385,8 @@ end subroutine writeincrement subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& - datestring,nhr_anal + datestring,nhr_anal, & + incsfcfileprefixes,fgsfcfileprefixes use constants, only: grav,qcmin use mpi use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -4257,12 +4418,17 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate integer :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind ! netcdf things - integer(i_kind) :: dimids3(3),nccount(3),ncstart(3) + integer(i_kind) :: dimids3(3),nccount(3),ncstart(3), dimids2(2) integer(i_kind) :: ncid_out, lon_dimid, lat_dimid, lev_dimid, ilev_dimid integer(i_kind) :: lonvarid, latvarid, levvarid, pfullvarid, ilevvarid, & hyaivarid, hybivarid, uvarid, vvarid, delpvarid, delzvarid, & tvarid, sphumvarid, liqwatvarid, o3varid, icvarid, & - rwmrvarid, snmrvarid, grlevarid + rwmrvarid, snmrvarid, grlevarid, & + tmp2mvarid, spfh2mvarid, soilt1varid, soilt2varid, & + soilt3varid, soilt4varid, slc1varid, slc2varid, & + slc3varid, slc4varid, maskvarid + integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind,soilt3_ind, & + soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind integer(i_kind) :: iadateout ! fixed fields such as lat, lon, levs @@ -4274,11 +4440,20 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ! increment real(r_kind), dimension(nlons*nlats) :: psinc, inc, ug, vg, work real(r_single), allocatable, dimension(:,:,:) :: inc3d, inc3d2, inc3dout + real(r_single), allocatable, dimension(:,:) :: inc2d, inc2dout real(r_single), allocatable, dimension(:,:,:) :: tv, tvanl, tmp, tmpanl, q, qanl real(r_single), allocatable, dimension(:,:,:) :: q2, qanl2 real(r_kind), allocatable, dimension(:,:) :: values_2d real(r_kind), allocatable, dimension(:) :: psges, delzb, values_1d + ! soil / snow mask (not fixed) + integer(i_kind), dimension(nlons,nlats) :: mask + + logical :: write_sfc_file, write_atm_file + real(r_double) :: t1,t2 + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + use_full_hydro = .false. clip = tiny_r_kind read(datestring,*) iadateout @@ -4317,6 +4492,9 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) enddo + if (write_atm_file ) then + + if (nproc == 0) t1 = mpi_wtime() ! loop through times and do the read ne = 1 backgroundloop: do nb=1,nbackgrounds @@ -4772,8 +4950,261 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate if (allocated(delzb)) deallocate(delzb) if (allocated(psges)) deallocate(psges) + !closing file + call nccheck_incr(nf90_close(ncid_out)) end do backgroundloop ! loop over backgrounds to write out + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement_pnc atm_file on root',t2-t1,'secs' + endif + end if ! if (write_atm_file) + + if (write_sfc_file ) then + + if (nproc == 0) t1 = mpi_wtime() + + tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'st1') + slc1_ind = getindex(vars2d, 'sl1') + soilt2_ind = getindex(vars2d, 'st2') + slc2_ind = getindex(vars2d, 'sl2') + soilt3_ind = getindex(vars2d, 'st3') + slc3_ind = getindex(vars2d, 'sl3') + soilt4_ind = getindex(vars2d, 'st4') + slc4_ind = getindex(vars2d, 'sl4') + + ! loop through times and do the read + ne = 1 + write(charnanal,'(i3.3)') nanal + sfcbackgroundloop: do nb=1,nbackgrounds + + if(no_inflate_flag) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"nimem"//charnanal + else + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"mem"//charnanal + end if + filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"mem"//charnanal + + !! note: only iope=0 is writing the outputs. Having all pes in iocomm write to a file slows it down. + !! + if (iope==0) then + dsfg = open_dataset(filenamein) + ! create the output netCDF increment file + call nccheck_incr(nf90_create(path=trim(filenameout), cmode=nf90_netcdf4,ncid=ncid_out)) + + ! create dimensions based on analysis resolution, not guess + call nccheck_incr(nf90_def_dim(ncid_out, "longitude", nlons, lon_dimid)) + call nccheck_incr(nf90_def_dim(ncid_out, "latitude", nlats, lat_dimid)) + dimids2 = (/ lon_dimid, lat_dimid /) + ! create variables + call nccheck_incr(nf90_def_var(ncid_out, "longitude", nf90_real,(/lon_dimid/), lonvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "latitude", nf90_real,(/lat_dimid/), latvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids2,tmp2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids2,spfh2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids2,soilt1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids2,soilt2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids2,soilt3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids2,soilt4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids2,slc1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids2,slc2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids2,slc3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids2,slc4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int,dimids2, maskvarid)) + ! place global attributes to serial calc_increment output + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI EnKF")) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", & + "global landsfc anal increment from writeincrement")) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "analysis_time",iadateout)) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global,"IAU_hour_from_guess", nhr_anal(nb))) + ! add units to lat/lon because that's what the calc_increment utility has + call nccheck_incr(nf90_put_att(ncid_out, lonvarid, "units","degrees_east")) + call nccheck_incr(nf90_put_att(ncid_out, latvarid, "units","degrees_north")) + ! end the netCDF file definition + call nccheck_incr(nf90_enddef(ncid_out)) + + ! longitudes + call read_vardata(dsfg, 'grid_xt', values_1d, errcode=iret) + deglons(:) = values_1d + call nccheck_incr(nf90_put_var(ncid_out, lonvarid, deglons, & + start = (/1/), count = (/nlons/))) + + call read_vardata(dsfg, 'grid_yt', values_1d, errcode=iret) + ! latitudes + do j=1,nlats + deglats(nlats-j+1) = values_1d(j) + end do + call nccheck_incr(nf90_put_var(ncid_out, latvarid, deglats, & + start = (/1/), count = (/nlats/))) + ! construct mask (1 - soil, 2 - snow, 0 - not snow) + ! note: same logic/threshold used in global_cycle to produce + ! mask on model grid. + call read_vardata(dsfg, 'soill1', values_2d, errcode=iret) + mask = 0 + do j=1,nlats + do i = 1, nlons + if (values_2d(i,j) .LT. 1.0) then + mask(i,nlats-j+1) = 1 + endif + enddo + end do + call read_vardata(dsfg, 'weasd', values_2d, errcode=iret) + do j=1,nlats + do i = 1, nlons + if (values_2d(i,j) .GT. 0.001) then + mask(i,nlats-j+1) = 2 + endif + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, maskvarid, mask, & + start = ncstart(1:2), count = nccount(1:2))) + + allocate(inc2d(nlons,nlats)) + allocate(inc2dout(nlons,nlats)) + + ! tmp2m increment + inc(:) = zero + if (tmp2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d) + tmp2m_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, tmp2mvarid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! spfh2m increment + inc(:) = zero + if (spfh2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+spfh2m_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, spfh2mvarid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt1 increment + inc(:) = zero + if (soilt1_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt1_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt1varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt2 increment + inc(:) = zero + if (soilt2_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt2_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt2varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt3 increment + inc(:) = zero + if (soilt3_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt3_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt3varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt4 increment + inc(:) = zero + if (soilt4_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt4_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt4varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc1 increment + inc(:) = zero + if (slc1_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc1_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc1varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc2 increment + inc(:) = zero + if (slc2_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc2_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc2varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc3 increment + inc(:) = zero + if (slc3_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc3_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc3varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc4 increment + inc(:) = zero + if (slc4_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc4_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc4varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + + call close_dataset(dsfg,errcode=iret) + if (iret/=0) then + write(6,*)'gridio/writeincrement_par: problem closing netcdf sfc fg dataset, iret=',iret + call stop2(23) + endif + ! deallocate things + deallocate(inc2d,inc2dout) + + call nccheck_incr(nf90_close(ncid_out)) + + end if + + end do sfcbackgroundloop ! loop over backgrounds to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement_pnc sfc_file on root',t2-t1,'secs' + endif + endif !write_Sfc + ! remove the sub communicators call mpi_barrier(iocomms(mem_pe(nproc)), iret) call mpi_comm_free(iocomms(mem_pe(nproc)), iret)