diff --git a/ci/spack.yaml b/ci/spack.yaml index a831de16ad..4374d8dd85 100644 --- a/ci/spack.yaml +++ b/ci/spack.yaml @@ -19,7 +19,7 @@ spack: - wrf-io@1.2.0 - ncio@1.1.2 - crtm@2.4.0 - - gsi-ncdiag@1.0.0 + - gsi-ncdiag@1.1.0 view: true concretizer: unify: true diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index bb2421c89c..4aa2613c63 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -131,7 +131,7 @@ subroutine init_controlvec() cvars3d(nc3d) = trim(adjustl(var)) clevels(nc3d) = ilev + clevels(nc3d-1) else - if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev + if (nproc .eq. 0) print *,'Error controlvec: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev call stop2(503) endif enddo @@ -212,7 +212,10 @@ subroutine read_control() ! read in whole control vector on i/o procs - keep in memory ! (needed in write_ensemble) allocate(grdin(npts,ncdim,nbackgrounds,nanals_per_iotask)) -allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask)) +! if only updating the sfc fields, qsat will not be calculated in readgriddata +! only allocate if needed. +q_ind = getindex(cvars3d, 'q') +if (q_ind > 0) allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask)) if (paranc) then if (nproc == 0) t1 = mpi_wtime() call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & @@ -225,7 +228,8 @@ subroutine read_control() fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) end if !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat) - if (use_qsatensmean) then + q_ind = getindex(cvars3d, 'q') + if (use_qsatensmean .and. q_ind>0 ) then allocate(qsatmean(npts,nlevs,nbackgrounds)) allocate(qsat_tmp(npts)) ! compute ensemble mean qsat @@ -257,7 +261,6 @@ subroutine read_control() ! print *,'min/max qsatmean proc',nproc,'=',& ! minval(qsatmean(:,:,nbackgrounds/2+1)),maxval(qsatmean(:,:,nbackgrounds/2+1)) !endif - q_ind = getindex(cvars3d, 'q') if (pseudo_rh .and. q_ind > 0) then if (use_qsatensmean) then do ne=1,nanals_per_iotask diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90 index c2e2b10f57..317ca2221c 100644 --- a/src/enkf/gridinfo_gfs.f90 +++ b/src/enkf/gridinfo_gfs.f90 @@ -66,7 +66,7 @@ module gridinfo ! supported variable names in anavinfo character(len=max_varname_length),public, dimension(13) :: vars3d_supported = (/'u ', 'v ', 'tv ', 'q ', 'oz ', 'cw ', 'tsen', 'prse', & 'ql ', 'qi ', 'qr ', 'qs ', 'qg '/) -character(len=max_varname_length),public, dimension(3) :: vars2d_supported = (/'ps ', 'pst', 'sst' /) +character(len=max_varname_length),public, dimension(13) :: vars2d_supported = (/'ps ', 'pst', 'sst', 't2m', 'q2m', 'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /) ! supported variable names in anavinfo contains diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 456f8126e2..5d560be3a8 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -35,6 +35,7 @@ module gridio ! a required input for EFSO calculations ! 2019-03-13 Add precipitation components ! 2019-07-10 Add convective clouds +! 2022-07-21 Draper: added read/write for sfc file for nc io (writeincrements, and readgridata) ! ! attributes: ! language: f95 @@ -100,12 +101,21 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & logical ice logical use_full_hydro integer(i_kind), allocatable, dimension(:) :: mem_pe, lev_pe1, lev_pe2, iocomms - integer(i_kind) :: iope, ionumproc, iolevs, krev + integer(i_kind) :: iope, ionumproc, iolevs, krev, ierr integer(i_kind) :: ncstart(3), nccount(3) ! mpi gatherv things integer(i_kind), allocatable, dimension(:) :: recvcounts, displs real(r_single), dimension(nlons,nlats,nlevs) :: ug3d_0, vg3d_0 + logical :: read_sfc_file, read_atm_file + + 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)) @@ -554,7 +564,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, type(sigio_head) :: sighead type(sigio_data) :: sigdata type(nemsio_gfile) :: gfile - type(Dataset) :: dset + type(Dataset) :: dset, dset_sfc type(Dimension) :: londim,latdim,levdim type(nemsio_gfile) :: gfilesfc @@ -562,14 +572,20 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, 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 + 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,iunitsig,iret,nb,i,idvc,nlonsin,nlatsin,nlevsin,ne,nanal integer(i_kind) :: nlonsin_sfc,nlatsin_sfc logical ice logical use_full_hydro + logical read_sfc_file, read_atm_file use_full_hydro = .false. + ! determine which files will be read in + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, read_sfc_file, read_atm_file) + ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -662,6 +678,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ! print *, 'ql: ', ql_ind, ', prse: ', prse_ind ! print *, 'ps: ', ps_ind, ', pst: ', pst_ind, ', sst: ', sst_ind ! endif + if (read_atm_file) then if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) @@ -731,7 +748,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*psg if (nanal .eq. 1) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) enddo - deallocate(ak,bk,values_2d) + deallocate(ak,bk) else vrtspec = sigdata%ps call sptez_s(vrtspec,psg,1) @@ -1100,10 +1117,134 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, deallocate(pressi,pslg) deallocate(psg) if (pst_ind > 0) deallocate(vmassdiv,pstend) + endif ! read_atm_file + if (use_gfs_nemsio) call nemsio_close(gfile,iret=iret) if (use_gfs_ncio) call close_dataset(dset) if (use_gfs_nemsio) call nemsio_close(gfilesfc,iret=iret) + if ( read_sfc_file ) then + + if ( .not. use_gfs_ncio ) then + write(6,*) 'griddio/griddata for sfc update vars only coded for nc io' + call stop2(23) + endif + if ( reducedgrid ) then + write(6,*) "reducedgrid=T interpolation not valid for writing sfc files" + call stop2(22) + 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') + + dset_sfc = open_dataset(filenamesfc) + ! 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 *,'error reading tmp2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + 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 *,'error reading spfh2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + 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 *,'error reading soilt1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + 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 *,'error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + 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 *,'error reading soilt3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + 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 *,'error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne)) + endif + if (slc1_ind > 0) then + call read_vardata(dset_sfc, 'slc1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne)) + endif + if (slc2_ind > 0) then + call read_vardata(dset_sfc, 'slc2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne)) + endif + if (slc3_ind > 0) then + call read_vardata(dset_sfc, 'slc3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne)) + endif + if (slc4_ind > 0) then + call read_vardata(dset_sfc, 'slc4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc4_ind,nb,ne)) + endif + + call close_dataset(dset_sfc) + + endif ! sfc read + + if ( allocated(values_2d) ) deallocate(values_2d) + end do backgroundloop ! loop over backgrounds to read in end do ensmemloop ! loop over ens members to read in @@ -1998,6 +2139,15 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n integer k,krev,nt,ierr,iunitsig,nb,i,ne,nanal logical :: nocompress + logical :: write_sfc_file, write_atm_file + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + + if (write_sfc_file .and. nproc==0 ) then + ! adding the sfc increments requires adjusting several other variables. This is done is a separate + ! program. + write(6,*)'gridio/writegriddata: not coded to write sfc analysis, use separate add_incr program instead' + endif nocompress = .true. if (nccompress) nocompress = .false. @@ -3402,7 +3552,7 @@ end subroutine writegriddata subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& - datestring,nhr_anal,write_ensmean + datestring,nhr_anal,write_ensmean, fgsfcfileprefixes,incsfcfileprefixes use constants, only: grav use mpi use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -3438,7 +3588,12 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, 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 + tvarid, sphumvarid, liqwatvarid, o3varid, icvarid, & + 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 @@ -3450,10 +3605,17 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! 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_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 + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + if ( write_atm_file) then use_full_hydro = .false. clip = tiny_r_kind read(datestring,*) iadateout @@ -3774,14 +3936,267 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & start = ncstart, count = nccount)) + call close_dataset(dsfg,errcode=iret) + if (iret/=0) then + write(6,*)'gridio/writeincrement gfs model: problem closing netcdf fg dataset, iret=',iret + call stop2(23) + endif ! deallocate things deallocate(inc3d,inc3d2,inc3dout) deallocate(tmp,tv,q,tmpanl,tvanl,qanl) - deallocate(delzb,psges) + if (allocated(delzb)) deallocate(delzb) + if (allocated(psges)) deallocate(psges) end do backgroundloop ! loop over backgrounds to read in end do ensmemloop ! loop over ens members to read in + endif ! write_atm_file + + if (write_sfc_file) then + + ne = 0 + sfcensmemloop: do nanal=nanal1,nanal2 + ne = ne + 1 + write(charnanal,'(i3.3)') nanal + sfcbackgroundloop: do nb=1,nbackgrounds + + if (nanal == 0 .and. write_ensmean) then + filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"ensmean" + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"ensmean" + else + 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 + endif + + ! 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)) + ! 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, dimids3(1:2), tmp2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids3(1:2), spfh2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids3(1:2), soilt1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids3(1:2), soilt2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids3(1:2), soilt3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids3(1:2), soilt4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids3(1:2), slc1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids3(1:2), slc2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids3(1:2), slc3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids3(1:2), slc4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids3(1:2), 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)) + + 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') + + dsfg = open_dataset(filenamein) + + ! 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, 'slc1', 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 gfs model: problem closing netcdf sfc fg dataset, iret=',iret + call stop2(23) + endif + ! deallocate things + deallocate(inc2d,inc2dout) + + end do sfcbackgroundloop ! loop over backgrounds to read in + end do sfcensmemloop ! loop over ens members to read in + + endif ! write_sfc_file + return contains @@ -4402,6 +4817,64 @@ end subroutine copyfromgrdin end subroutine writeincrement_pnc + subroutine set_ncio_file_flags(vars3d, n3d, vars2d, n2d, sfc_file, atm_file) + ! determine if variables are in sfc and/or atm file, for ncio case. + character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d + character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d + integer, intent(in) :: n2d, n3d + logical, intent(out) :: sfc_file, atm_file + + 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 + 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 + + ! atmos file variables + u_ind = getindex(vars3d, 'u') !< indices in the state or control var arrays + v_ind = getindex(vars3d, 'v') ! U and V (3D) + tv_ind = getindex(vars3d, 'tv') ! Tv (3D) + q_ind = getindex(vars3d, 'q') ! Q (3D) + oz_ind = getindex(vars3d, 'oz') ! Oz (3D) + cw_ind = getindex(vars3d, 'cw') ! CW (3D) + tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) + ql_ind = getindex(vars3d, 'ql') ! QL (3D) + qi_ind = getindex(vars3d, 'qi') ! QI (3D) + prse_ind = getindex(vars3d, 'prse') + qr_ind = getindex(vars3d, 'qr') ! QR (3D) + qs_ind = getindex(vars3d, 'qs') ! QS (3D) + qg_ind = getindex(vars3d, 'qg') ! QG (3D) + ps_ind = getindex(vars2d, 'ps') ! Ps (2D) + pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of + ! old logical massbal_adjust, if non-zero + sst_ind = getindex(vars2d, 'sst') ! is this really in the atmos file? + + ! for nc gfs io determine if requested variables are in sfc and/or atmos file + atm_file = ( u_ind>0 .or. v_ind>0 .or. tv_ind>0 .or. q_ind>0 .or. sst_ind>0 .or. & + oz_ind>0 .or. cw_ind>0 .or. tsen_ind>0 .or. ql_ind>0 .or. & + qi_ind>0 .or. prse_ind>0 .or. qr_ind>0 .or. qs_ind>0 .or. qg_ind>0 ) + + ! sfc file 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') + + sfc_file = ( tmp2m_ind > 0 .or. spfh2m_ind > 0 .or. soilt1_ind > 0 .or. & + slc1_ind > 0 .or. soilt2_ind > 0 .or. slc2_ind > 0 .or. & + soilt3_ind > 0 .or. slc3_ind > 0 .or. soilt4_ind > 0 .or. & + slc4_ind > 0 ) + + end subroutine set_ncio_file_flags + + logical function checkfield(field,fields,nrec) result(hasfield) use nemsio_module, only: nemsio_charkind integer, intent(in) :: nrec diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 593e5a5ec4..b21d88abd0 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -85,7 +85,9 @@ module params character(len=120),dimension(7),public :: statefileprefixes character(len=120),dimension(7),public :: statesfcfileprefixes character(len=120),dimension(7),public :: anlfileprefixes +character(len=120),dimension(7),public :: anlsfcfileprefixes character(len=120),dimension(7),public :: incfileprefixes +character(len=120),dimension(7),public :: incsfcfileprefixes ! analysis date string (YYYYMMDDHH) character(len=10), public :: datestring ! Hour for datestring @@ -266,7 +268,7 @@ module params lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,& lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh,& fgfileprefixes,fgsfcfileprefixes,anlfileprefixes, & - incfileprefixes, & + anlsfcfileprefixes,incfileprefixes,incsfcfileprefixes,& statefileprefixes,statesfcfileprefixes, & covl_minfact,covl_efold,lupd_obspace_serial,letkf_novlocal,& analpertwtnh,analpertwtsh,analpertwttr,sprd_tol,& @@ -460,8 +462,8 @@ subroutine read_namelist() ! Initialize first-guess and analysis file name prefixes. ! (blank means use default names) fgfileprefixes = ''; anlfileprefixes=''; statefileprefixes='' -fgsfcfileprefixes = ''; statesfcfileprefixes='' -incfileprefixes = '' +anlsfcfileprefixes=''; fgsfcfileprefixes = ''; statesfcfileprefixes='' +incfileprefixes = ''; incsfcfileprefixes = '' ! option for including convective clouds in the all-sky cnvw_option=.false. @@ -720,7 +722,7 @@ subroutine read_namelist() endif endif if (trim(fgsfcfileprefixes(nbackgrounds+1)) .eq. "") then - fgsfcfileprefixes(nbackgrounds+1)="sfgsfc_"//datestring//"_fhr"//charfhr_anal(nbackgrounds+1)//"_" + fgsfcfileprefixes(nbackgrounds+1)="bfg_"//datestring//"_fhr"//charfhr_anal(nbackgrounds+1)//"_" end if nbackgrounds = nbackgrounds+1 end do @@ -742,7 +744,7 @@ subroutine read_namelist() endif endif if (trim(statesfcfileprefixes(nstatefields+1)) .eq. "") then - statesfcfileprefixes(nstatefields+1)="sfgsfc_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_" + statesfcfileprefixes(nstatefields+1)="bfg_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_" end if nstatefields = nstatefields+1 end do @@ -762,6 +764,23 @@ subroutine read_namelist() incfileprefixes(nb)="incr_"//datestring//"_fhr"//charfhr_anal(nb)//"_" ! else ! anlfileprefixes(nb)="sanl_"//datestring//"_" +! endif + endif + endif + if (trim(anlsfcfileprefixes(nb)) .eq. "") then + ! default analysis file prefix + if (regional) then + if (nbackgrounds > 1) then + anlsfcfileprefixes(nb)="sfc_analysis_fhr"//charfhr_anal(nb)//"." + else + anlsfcfileprefixes(nb)="sfc_analysis." + endif + else ! global +! if (nbackgrounds > 1) then + anlsfcfileprefixes(nb)="banl_"//datestring//"_fhr"//charfhr_anal(nb)//"_" + incsfcfileprefixes(nb)="sfcincr_"//datestring//"_fhr"//charfhr_anal(nb)//"_" +! else +! anlfileprefixes(nb)="sanl_"//datestring//"_" ! endif endif endif diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index f3b1e38f04..d1f4ec3ff8 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -24,6 +24,7 @@ module readconvobs ! reflectivity and radial velocity assimilation. POC: xuguang.wang@ou.edu ! 2017-12-13 shlyaeva - added netcdf diag read/write capability ! 2019-03-21 CAPS(C. Tong) - added direct reflectivity DA capability +! 2022-03-23 draper - added option to not scale qobs by forecast qsat. ! ! attributes: ! language: f95 @@ -32,7 +33,8 @@ module readconvobs use kinds, only: r_kind,i_kind,r_single,r_double use constants, only: one,zero,deg2rad -use params, only: npefiles, netcdf_diag, modelspace_vloc, l_use_enkf_directZDA +use params, only: npefiles, netcdf_diag, modelspace_vloc, & + l_use_enkf_directZDA implicit none private @@ -329,7 +331,6 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) call nc_diag_read_close(obsfile) - num_obs_totdiag = num_obs_totdiag + nobs_curr do i = 1, nobs_curr diff --git a/src/enkf/statevec.f90 b/src/enkf/statevec.f90 index d1be91af3c..5ad70346aa 100644 --- a/src/enkf/statevec.f90 +++ b/src/enkf/statevec.f90 @@ -14,7 +14,7 @@ module statevec ! ! Public Variables: ! nanals: (integer scalar) number of ensemble members (from module params) -! nlevs: number of analysis vertical levels (from module params). +! nlevs: number of analysis atmos vertical levels (from module params). ! ns3d: number of 3D variables ! ns2d: number of 2D variables ! svars3d: names of 3D variables @@ -120,7 +120,7 @@ subroutine init_statevec() svars3d(ns3d)=trim(adjustl(var)) slevels(ns3d)=ilev + slevels(ns3d-1) else - if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev + if (nproc .eq. 0) print *,'Error statevec: - only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev call stop2(503) endif enddo diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index 0216b95fb6..ffdcd90c79 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -411,6 +411,77 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & return end subroutine general_reload2 + +! 2m reload +subroutine general_reload_sfc(grd,g_t2m, g_q2m,g_ps,icount,iflag,work) +! !USES: + use kinds, only: r_kind,i_kind + use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype + use general_sub2grid_mod, only: sub2grid_info + + implicit none +! !INPUT PARAMETERS: + + type(sub2grid_info), intent(in ) :: grd + integer(i_kind), intent(inout) :: icount + integer(i_kind),dimension(npe), intent(inout) :: iflag + real(r_kind),dimension(grd%itotsub),intent(in ) :: work + +! !OUTPUT PARAMETERS: + + real(r_kind),dimension(grd%lat2,grd%lon2), intent( out) :: g_t2m,& + g_q2m, g_ps + +! !DESCRIPTION: version of general_reload, for 2m variables. +! +! !REVISION HISTORY: +! 2023-03-2 Draper +!------------------------------------------------------------------------- + + integer(i_kind) i,j,ij,k + real(r_kind),dimension(grd%lat2*grd%lon2,npe):: sub + + call mpi_alltoallv(work,grd%sendcounts_s,grd%sdispls_s,mpi_rtype,& + sub,grd%recvcounts_s,grd%rdispls_s,mpi_rtype,& + mpi_comm_world,ierror) + +!$omp parallel do schedule(dynamic,1) private(k,i,j,ij) + + do k=1,icount + if ( iflag(k) == 2 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_t2m(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 3 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_q2m(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 4 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ps(i,j)=sub(ij,k) + enddo + enddo + endif + enddo ! do k=1,icount + + icount=0 + iflag=0 + + return + +end subroutine general_reload_sfc + end module gfsreadmod subroutine general_read_gfsatm(grd,sp_a,sp_b,filename,uvflag,vordivflag,zflag, & @@ -431,6 +502,7 @@ subroutine general_read_gfsatm(grd,sp_a,sp_b,filename,uvflag,vordivflag,zflag, & ! 2014-11-30 todling - genelize interface to handle bundle instead of fields; ! internal code should be generalized ! 2014-12-03 derber - introduce vordivflag, zflag and optimize routines +! 2023-03-23 draper - added option to read sfc files (for 2m variables) ! ! input argument list: ! grd - structure variable containing information about grid @@ -1892,7 +1964,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & use gsi_bundlemod, only: gsi_bundlegetpointer use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& close_dataset, get_dim, read_vardata,get_idate_from_time_units - use gfsreadmod, only: general_reload + use gfsreadmod, only: general_reload, general_reload_sfc implicit none @@ -1910,6 +1982,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & real(r_kind),pointer,dimension(:,:) :: ptr2d real(r_kind),pointer,dimension(:,:,:) :: ptr3d real(r_kind),pointer,dimension(:,:) :: g_ps + real(r_kind),pointer,dimension(:,:) :: g_t2m, g_q2m real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& g_cwmr,g_q,g_oz,g_tv @@ -1942,10 +2015,9 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & logical :: procuse,diff_res,eqspace type(egrid2agrid_parm) :: p_high logical,dimension(1) :: vector - type(Dataset) :: atmges + type(Dataset) :: filges type(Dimension) :: ncdim - - + logical :: read_2m, read_z !****************************************************************************** ! Initialize variables used below @@ -1959,6 +2031,19 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & mype_use=-1 icount=0 procuse=.false. + + if (filename(1:3) == 'sfc') then + read_2m = .true. + read_z = .false. + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading 2m variables from ', trim(filename) + else + read_2m = .false. + read_z = zflag + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading atmos variables from ', trim(filename) + endif + if ( mype == 0 ) procuse = .true. do i=1,npe if ( grd%recvcounts_s(i-1) > 0 ) then @@ -1992,20 +2077,21 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & if ( procuse ) then - atmges = open_dataset(filename, paropen=.true., mpicomm=mpi_comm_read) + filges = open_dataset(filename, paropen=.true., mpicomm=mpi_comm_read) ! get dimension sizes - ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len - ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len - ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len + ncdim = get_dim(filges, 'grid_xt'); lonb = ncdim%len + ncdim = get_dim(filges, 'grid_yt'); latb = ncdim%len + if (.not. read_2m) & + ncdim = get_dim(filges, 'pfull'); levs = ncdim%len ! get time information - idate = get_idate_from_time_units(atmges) + idate = get_idate_from_time_units(filges) odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year - call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later + call read_vardata(filges, 'time', fhour) ! might need to change this to attribute later ! depends on model changes from ! Jeff Whitaker fhour = float(nint(fhour)) @@ -2030,11 +2116,13 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & trim(my_name),grd%nlon,lonb !call stop2(101) endif - if ( levs /= grd%nsig ) then - if ( mype == 0 ) write(6, & - '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & - trim(my_name),grd%nsig,levs - call stop2(101) + if (.not. read_2m) then + if ( levs /= grd%nsig ) then + if ( mype == 0 ) write(6, & + '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & + trim(my_name),grd%nsig,levs + call stop2(101) + endif endif allocate( spec_vor(sp_a%nc), spec_div(sp_a%nc) ) @@ -2047,8 +2135,8 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & allocate(rwork3d1(lonb,latb,1)) allocate(rwork2d(lonb,latb)) allocate(rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb)) - call read_vardata(atmges, 'grid_xt', rlons_tmp) - call read_vardata(atmges, 'grid_yt', rlats_tmp) + call read_vardata(filges, 'grid_xt', rlons_tmp) + call read_vardata(filges, 'grid_yt', rlats_tmp) do j=1,latb rlats(latb+2-j)=deg2rad*rlats_tmp(j) end do @@ -2073,57 +2161,73 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & endif ! if ( procuse ) ! Get pointer to relevant variables (this should be made flexible and general) - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both sf and div' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both vp and vor' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'cannot handle having both t and tv' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - istatus=0 - call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier - if ( istatus /= 0 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_nems: ERROR' - write(6,*) 'Missing some of the required fields' - write(6,*) 'Aborting ... ' - endif - call stop2(999) + if (.not. read_2m) then + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nc: ERROR' + write(6,*) 'cannot handle having both sf and div' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nc: ERROR' + write(6,*) 'cannot handle having both vp and vor' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nc: ERROR' + write(6,*) 'cannot handle having both t and tv' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nc: ERROR' + write(6,*) 'Missing some of the required fields' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + else ! read 2m vars + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'t2m',g_t2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q2m',g_q2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_nc: ERROR' + write(6,*) 'Missing 2m required variables' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif endif allocate(g_u(grd%lat2,grd%lon2,grd%nsig),g_v(grd%lat2,grd%lon2,grd%nsig)) allocate(g_z(grd%lat2,grd%lon2)) @@ -2135,8 +2239,8 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & ! Once on the grid, fields need to be scattered from the full domain to ! sub-domains. - ! Only read Terrain when zflag is true. - if ( zflag ) then + ! Only read Terrain when read_z is true. + if ( read_z ) then icount=icount+1 iflag(icount)=1 @@ -2145,7 +2249,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & ! Terrain: spectral --> grid transform, scatter to all mpi tasks if (mype==mype_use(icount)) then ! read hs - call read_vardata(atmges, 'hgtsfc', rwork2d) + call read_vardata(filges, 'hgtsfc', rwork2d) if ( diff_res ) then grid_b=rwork2d vector(1)=.false. @@ -2161,415 +2265,498 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & call general_fill_ns(grd,grid,work) endif endif - if ( icount == icm ) then + if ( icount == icm ) then call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & icount,iflag,ilev,work,uvflag,vordivflag) endif endif - icount=icount+1 - iflag(icount)=2 - ilev(icount)=1 + if (.not. read_2m) then + + icount=icount+1 + iflag(icount)=2 + ilev(icount)=1 + + ! Surface pressure: same procedure as terrain + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(filges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + + ! Thermodynamic variable: s-->g transform, communicate to all tasks + ! For multilevel fields, each task handles a given level. Periodic + ! mpi_alltoallv calls communicate the grids to all mpi tasks. + ! Finally, the grids are loaded into guess arrays used later in the + ! code. + + do k=1,nlevs + + icount=icount+1 + iflag(icount)=3 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'spfh', rwork3d1, nslice=kr, slicedim=3) + call read_vardata(filges, 'tmp', rwork3d0, nslice=kr, slicedim=3) + rwork2d = rwork3d0(:,:,1) * (one+fv*rwork3d1(:,:,1)) + if ( diff_res ) then + grid_b=rwork2d + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + if ( vordivflag .or. .not. uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=4 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Vorticity + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_vor(grd%nlon,nlatm2)) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_vor,work) + deallocate(grid_vor) + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + + end do + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=5 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Divergence + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_div(grd%nlon,nlatm2) ) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_div,work) + deallocate(grid_div) + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + + end do + endif ! if ( vordivflag .or. .not. uvflag ) + if ( uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=6 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + + icount=icount+1 + iflag(icount)=7 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! V + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + ! Note work_v and work are switched because output must be in work. + call general_filluv_ns(grd,slons,clons,grid,grid_v,work_v,work) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + endif ! if ( uvflag ) + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=8 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! Specific humidity + call read_vardata(filges, 'spfh', rwork3d0, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=9 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'o3mr', rwork3d0, nslice=kr, slicedim=3) + ! Ozone mixing ratio + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + do k=1,nlevs + icount=icount+1 + iflag(icount)=10 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'clwmr', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'icmr', rwork3d1, nslice=kr, slicedim=3) + ! Cloud condensate mixing ratio. + rwork2d = rwork3d0(:,:,1)+rwork3d1(:,:,1) + if ( diff_res ) then + grid_b=rwork2d + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + + endif + + if ( icount == icm .or. k == nlevs ) then + call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & + icount,iflag,ilev,work,uvflag,vordivflag) + endif + + enddo ! do k=1,nlevs + else ! read_2m + + icount=icount+1 + iflag(icount)=2 + + ! 2m temperature from sfc file + if (mype==mype_use(icount)) then + call read_vardata(filges, 'tmp2m', rwork2d) + + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + icount=icount + 1 + iflag(icount)=3 + + ! 2m humidity from sfc file + if (mype==mype_use(icount)) then + call read_vardata(filges, 'spfh2m', rwork2d) + + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + icount=icount + 1 + iflag(icount)=4 + + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(filges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + ! not using all procs. doesn't trigger. todo: figure out trigger + ! for when reading fewer vars. + !if ( icount == icm ) then + call general_reload_sfc(grd,g_t2m, g_q2m, g_ps, icount,iflag,work) + !endif + + endif ! read_2m - ! Surface pressure: same procedure as terrain - if (mype==mype_use(icount)) then - ! read ps - call read_vardata(atmges, 'pressfc', rwork2d) - rwork2d = r0_001*rwork2d ! convert Pa to cb - if ( diff_res ) then - vector(1)=.false. - grid_b=rwork2d - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork2d - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) + if ( procuse ) then + if ( diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) + call destroy_egrid2agrid(p_high) + deallocate(spec_div,spec_vor) + deallocate(rwork3d1,rwork3d0,clons,slons) + deallocate(rwork2d) + deallocate(grid,grid_v) + call close_dataset(filges) endif + deallocate(work) - ! Thermodynamic variable: s-->g transform, communicate to all tasks - ! For multilevel fields, each task handles a given level. Periodic - ! mpi_alltoallv calls communicate the grids to all mpi tasks. - ! Finally, the grids are loaded into guess arrays used later in the - ! code. - - do k=1,nlevs + ! Convert dry temperature to virtual temperature + !do k=1,grd%nsig + ! do j=1,grd%lon2 + ! do i=1,grd%lat2 + ! g_tv(i,j,k) = g_tv(i,j,k)*(one+fv*g_q(i,j,k)) + ! enddo + ! enddo + !enddo - icount=icount+1 - iflag(icount)=3 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip + ! Load u->div and v->vor slot when uv are used instead + if ( .not. read_2m ) then + if ( uvflag ) then + call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif ! read_2m + if (read_z) then + call gsi_bundlegetpointer(gfs_bundle,'z' ,ptr2d,ier) + if ( ier == 0 ) ptr2d=g_z + endif - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'spfh', rwork3d1, nslice=kr, slicedim=3) - call read_vardata(atmges, 'tmp', rwork3d0, nslice=kr, slicedim=3) - rwork2d = rwork3d0(:,:,1) * (one+fv*rwork3d1(:,:,1)) - if ( diff_res ) then - grid_b=rwork2d - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork2d - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do + ! Clean up + deallocate(g_z) + deallocate(g_u,g_v) - if ( vordivflag .or. .not. uvflag ) then - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=4 - ilev(icount)=k + ! Print date/time stamp + if ( mype == 0 ) then + write(6,700) lonb,latb,nlevs,grd%nlon,nlatm2,& + fhour,odate,trim(filename) +700 format('GENERAL_READ_GFSATM_NC: read lonb,latb,levs=',& + 3i6,', scatter nlon,nlat=',2i6,', hour=',f6.1,', idate=',4i5,1x,a) + endif - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - ! Vorticity - ! Convert grid u,v to div and vor - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work_v(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid_v(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - allocate( grid_vor(grd%nlon,nlatm2)) - call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) - call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) - ! Load values into rows for south and north pole - call general_fill_ns(grd,grid_vor,work) - deallocate(grid_vor) - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif + return - end do - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip +end subroutine general_read_gfsatm_nc - icount=icount+1 - iflag(icount)=5 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - ! Divergence - ! Convert grid u,v to div and vor - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work_v(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid_v(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - allocate( grid_div(grd%nlon,nlatm2) ) - call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) - call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) - ! Load values into rows for south and north pole - call general_fill_ns(grd,grid_div,work) - deallocate(grid_div) - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - - end do - endif ! if ( vordivflag .or. .not. uvflag ) - if ( uvflag ) then - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=6 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - - icount=icount+1 - iflag(icount)=7 - ilev(icount)=k - - if (mype==mype_use(icount)) then - ! V - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - ! Note work_v and work are switched because output must be in work. - call general_filluv_ns(grd,slons,clons,grid,grid_v,work_v,work) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - endif ! if ( uvflag ) - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=8 - ilev(icount)=k - - if (mype==mype_use(icount)) then - ! Specific humidity - call read_vardata(atmges, 'spfh', rwork3d0, nslice=kr, slicedim=3) - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - - icount=icount+1 - iflag(icount)=9 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'o3mr', rwork3d0, nslice=kr, slicedim=3) - ! Ozone mixing ratio - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - - do k=1,nlevs - icount=icount+1 - iflag(icount)=10 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'clwmr', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'icmr', rwork3d1, nslice=kr, slicedim=3) - ! Cloud condensate mixing ratio. - rwork2d = rwork3d0(:,:,1)+rwork3d1(:,:,1) - if ( diff_res ) then - grid_b=rwork2d - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork2d - call general_fill_ns(grd,grid,work) - endif - - endif - - if ( icount == icm .or. k == nlevs ) then - call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & - icount,iflag,ilev,work,uvflag,vordivflag) - endif - - enddo ! do k=1,nlevs - - if ( procuse ) then - if ( diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) - call destroy_egrid2agrid(p_high) - deallocate(spec_div,spec_vor) - deallocate(rwork3d1,rwork3d0,clons,slons) - deallocate(rwork2d) - deallocate(grid,grid_v) - call close_dataset(atmges) - endif - deallocate(work) - - ! Convert dry temperature to virtual temperature - !do k=1,grd%nsig - ! do j=1,grd%lon2 - ! do i=1,grd%lat2 - ! g_tv(i,j,k) = g_tv(i,j,k)*(one+fv*g_q(i,j,k)) - ! enddo - ! enddo - !enddo - - ! Load u->div and v->vor slot when uv are used instead - if ( uvflag ) then - call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - endif - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - if (zflag) then - call gsi_bundlegetpointer(gfs_bundle,'z' ,ptr2d,ier) - if ( ier == 0 ) ptr2d=g_z - endif - - ! Clean up - deallocate(g_z) - deallocate(g_u,g_v) - - ! Print date/time stamp - if ( mype == 0 ) then - write(6,700) lonb,latb,nlevs,grd%nlon,nlatm2,& - fhour,odate,trim(filename) -700 format('GENERAL_READ_GFSATM_NC: read lonb,latb,levs=',& - 3i6,', scatter nlon,nlat=',2i6,', hour=',f6.1,', idate=',4i5,1x,a) - endif - - return - -end subroutine general_read_gfsatm_nc subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & gfs_bundle,iret_read) !$$$ subprogram documentation block @@ -2618,7 +2805,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z use gsi_bundlemod, only: gsi_bundlegetpointer use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& close_dataset, get_dim, read_vardata,get_idate_from_time_units - use gfsreadmod, only: general_reload2 + use gfsreadmod, only: general_reload2, general_reload_sfc use ncepnems_io, only: imp_physics implicit none @@ -2637,6 +2824,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z real(r_kind),pointer,dimension(:,:) :: ptr2d real(r_kind),pointer,dimension(:,:,:) :: ptr3d real(r_kind),pointer,dimension(:,:) :: g_ps + real(r_kind),pointer,dimension(:,:) :: g_t2m, g_q2m real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& g_q,g_oz,g_tv real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg @@ -2668,8 +2856,9 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z logical :: procuse,diff_res,eqspace type(egrid2agrid_parm) :: p_high logical,dimension(1) :: vector - type(Dataset) :: atmges + type(Dataset) :: filges type(Dimension) :: ncdim + logical :: read_2m, read_z @@ -2685,6 +2874,19 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z mype_use=-1 icount=0 procuse=.false. + + if (filename(1:3) == 'sfc') then + read_2m = .true. + read_z = .false. + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading 2m variables from ', trim(filename) + else + read_2m = .false. + read_z = zflag + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading atmos variables from ', trim(filename) + endif + if ( mype == 0 ) procuse = .true. do i=1,npe if ( grd%recvcounts_s(i-1) > 0 ) then @@ -2694,26 +2896,26 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif enddo icm=icount - allocate( work(grd%itotsub),work_v(grd%itotsub) ) + allocate( work(grd%itotsub)) work=zero - work_v=zero if ( procuse ) then - atmges = open_dataset(filename, paropen=.true.) + filges = open_dataset(filename, paropen=.true.) ! get dimension sizes - ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len - ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len - ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len + ncdim = get_dim(filges, 'grid_xt'); lonb = ncdim%len + ncdim = get_dim(filges, 'grid_yt'); latb = ncdim%len + if (.not. read_2m) & + ncdim = get_dim(filges, 'pfull'); levs = ncdim%len ! get time information - idate = get_idate_from_time_units(atmges) + idate = get_idate_from_time_units(filges) odate(1) = idate(4) !hour odate(2) = idate(2) !month odate(3) = idate(3) !day odate(4) = idate(1) !year - call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later + call read_vardata(filges, 'time', fhour) ! might need to change this to attribute later ! depends on model changes from ! Jeff Whitaker fhour = float(nint(fhour)) @@ -2738,11 +2940,13 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z trim(my_name),grd%nlon,lonb !call stop2(101) endif - if ( levs /= grd%nsig ) then - if ( mype == 0 ) write(6, & - '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & - trim(my_name),grd%nsig,levs - call stop2(101) + if (.not. read_2m) then + if ( levs /= grd%nsig ) then + if ( mype == 0 ) write(6, & + '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & + trim(my_name),grd%nsig,levs + call stop2(101) + endif endif allocate( spec_vor(sp_a%nc), spec_div(sp_a%nc) ) @@ -2755,8 +2959,8 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z allocate(rwork3d1(lonb,latb,1)) allocate(rwork2d(lonb,latb)) allocate(rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb)) - call read_vardata(atmges, 'grid_xt', rlons_tmp) - call read_vardata(atmges, 'grid_yt', rlats_tmp) + call read_vardata(filges, 'grid_xt', rlons_tmp) + call read_vardata(filges, 'grid_yt', rlats_tmp) do j=1,latb rlats(latb+2-j)=deg2rad*rlats_tmp(j) end do @@ -2781,64 +2985,79 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif ! if ( procuse ) ! Get pointer to relevant variables (this should be made flexible and general) - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' - write(6,*) 'cannot handle having both sf and div' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' - write(6,*) 'cannot handle having both vp and vor' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - iredundant=0 - call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) - if ( ier == 0 ) iredundant = iredundant + 1 - if ( iredundant==2 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' - write(6,*) 'cannot handle having both t and tv' - write(6,*) 'Aborting ... ' - endif - call stop2(999) - endif - istatus=0 - call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier - call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier -! call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier - istatus1=0 - call gsi_bundlegetpointer(gfs_bundle,'ql',g_ql ,ier);istatus1=istatus1+ier - call gsi_bundlegetpointer(gfs_bundle,'qi',g_qi ,ier);istatus1=istatus1+ier - call gsi_bundlegetpointer(gfs_bundle,'qr',g_qr ,ier);istatus1=istatus1+ier - call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs ,ier);istatus1=istatus1+ier - call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg ,ier);istatus1=istatus1+ier -! call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf ,ier);istatus1=istatus1+ier - if ( istatus1 /= 0 ) then - if ( mype == 0 ) then - write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' - write(6,*) 'Missing some of the required hydrometeor fields for imp_physics = ', imp_physics - write(6,*) 'Aborting ... ' - endif - call stop2(999) + if (.not. read_2m) then + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both sf and div' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both vp and vor' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both t and tv' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier + ! call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier + istatus1=0 + call gsi_bundlegetpointer(gfs_bundle,'ql',g_ql ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qi',g_qi ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qr',g_qr ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg ,ier);istatus1=istatus1+ier + ! call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf ,ier);istatus1=istatus1+ier + if ( istatus1 /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'Missing some of the required hydrometeor fields for imp_physics = ', imp_physics + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + else ! read 2m vars + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'t2m',g_t2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q2m',g_q2m ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + if ( istatus /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'Missing 2m required variables' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif endif allocate(g_u(grd%lat2,grd%lon2,grd%nsig),g_v(grd%lat2,grd%lon2,grd%nsig)) @@ -2853,7 +3072,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z ! Only read Terrain when zflag is true. - if ( zflag ) then + if ( read_z ) then icount=icount+1 iflag(icount)=1 @@ -2862,7 +3081,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z ! Terrain: spectral --> grid transform, scatter to all mpi tasks if (mype==mype_use(icount)) then ! read hs - call read_vardata(atmges, 'hgtsfc', rwork2d) + call read_vardata(filges, 'hgtsfc', rwork2d) if ( diff_res ) then grid_b=rwork2d vector(1)=.false. @@ -2884,465 +3103,470 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif endif - icount=icount+1 - iflag(icount)=2 - ilev(icount)=1 - - ! Surface pressure: same procedure as terrain - if (mype==mype_use(icount)) then - ! read ps - call read_vardata(atmges, 'pressfc', rwork2d) - rwork2d = r0_001*rwork2d ! convert Pa to cb - if ( diff_res ) then - vector(1)=.false. - grid_b=rwork2d - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork2d - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - - ! Thermodynamic variable: s-->g transform, communicate to all tasks - ! For multilevel fields, each task handles a given level. Periodic - ! mpi_alltoallv calls communicate the grids to all mpi tasks. - ! Finally, the grids are loaded into guess arrays used later in the - ! code. - - do k=1,nlevs - - icount=icount+1 - iflag(icount)=3 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'spfh', rwork3d1, nslice=kr, slicedim=3) - call read_vardata(atmges, 'tmp', rwork3d0, nslice=kr, slicedim=3) - rwork2d = rwork3d0(:,:,1) * (one+fv*rwork3d1(:,:,1)) - if ( diff_res ) then - grid_b=rwork2d - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork2d - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - - if ( vordivflag .or. .not. uvflag ) then - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=4 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - ! Vorticity - ! Convert grid u,v to div and vor - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work_v(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid_v(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - allocate( grid_vor(grd%nlon,nlatm2)) - call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) - call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) - ! Load values into rows for south and north pole - call general_fill_ns(grd,grid_vor,work) - deallocate(grid_vor) - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - - icount=icount+1 - iflag(icount)=5 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - ! Divergence - ! Convert grid u,v to div and vor - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work_v(kk)=grid2(i,j,1) - enddo - do j=1,grd%nlon - do i=2,grd%nlat-1 - grid_v(j,grd%nlat-i)=grid2(i,j,1) - enddo - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - allocate( grid_div(grd%nlon,nlatm2) ) - call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) - call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) - ! Load values into rows for south and north pole - call general_fill_ns(grd,grid_div,work) - deallocate(grid_div) - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - - end do - endif ! if ( vordivflag .or. .not. uvflag ) - - if ( uvflag ) then - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=6 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - - icount=icount+1 - iflag(icount)=7 - ilev(icount)=k - - if (mype==mype_use(icount)) then - ! V - call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) - call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) - if ( diff_res ) then - grid_b = rwork3d0(:,:,1) - grid_b2 = rwork3d1(:,:,1) - vector(1)=.true. - call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) - call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - grid_v = rwork3d1(:,:,1) - ! Note work_v and work are switched because output must be in work. - call general_filluv_ns(grd,slons,clons,grid,grid_v,work_v,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - endif ! if ( uvflag ) - - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - icount=icount+1 - iflag(icount)=8 - ilev(icount)=k - - if (mype==mype_use(icount)) then - ! Specific humidity - call read_vardata(atmges, 'spfh', rwork3d0, nslice=kr, slicedim=3) - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid = rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - - do k=1,nlevs - kr = levs+1-k ! netcdf is top to bottom, need to flip - - icount=icount+1 - iflag(icount)=9 - ilev(icount)=k - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'o3mr', rwork3d0, nslice=kr, slicedim=3) - ! Ozone mixing ratio - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - end do - - do k=1,nlevs - icount=icount+1 - iflag(icount)=10 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'clwmr', rwork3d0, nslice=kr, slicedim=3) - ! Cloud liquid water mixing ratio. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - enddo ! do k=1,nlevs - - do k=1,nlevs - icount=icount+1 - iflag(icount)=11 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'icmr', rwork3d0, nslice=kr, slicedim=3) - ! Cloud ice water mixing ratio. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - enddo ! do k=1,nlevs - - do k=1,nlevs - icount=icount+1 - iflag(icount)=12 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'rwmr', rwork3d0, nslice=kr, slicedim=3) - ! Rain water mixing ratio. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - enddo ! do k=1,nlevs - - do k=1,nlevs - icount=icount+1 - iflag(icount)=13 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'snmr', rwork3d0, nslice=kr, slicedim=3) - ! Snow water mixing ratio. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - enddo ! do k=1,nlevs - - do k=1,nlevs - icount=icount+1 - iflag(icount)=14 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip - - if (mype==mype_use(icount)) then - call read_vardata(atmges, 'grle', rwork3d0, nslice=kr, slicedim=3) - ! Graupel mixing ratio. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm .or. k==nlevs) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) - endif - enddo ! do k=1,nlevs + if (.not. read_2m) then + + allocate( work_v(grd%itotsub) ) + work_v=zero + + icount=icount+1 + iflag(icount)=2 + ilev(icount)=1 + + ! Surface pressure: same procedure as terrain + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(filges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + ! Thermodynamic variable: s-->g transform, communicate to all tasks + ! For multilevel fields, each task handles a given level. Periodic + ! mpi_alltoallv calls communicate the grids to all mpi tasks. + ! Finally, the grids are loaded into guess arrays used later in the + ! code. + + do k=1,nlevs + + icount=icount+1 + iflag(icount)=3 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'spfh', rwork3d1, nslice=kr, slicedim=3) + call read_vardata(filges, 'tmp', rwork3d0, nslice=kr, slicedim=3) + rwork2d = rwork3d0(:,:,1) * (one+fv*rwork3d1(:,:,1)) + if ( diff_res ) then + grid_b=rwork2d + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + if ( vordivflag .or. .not. uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=4 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Vorticity + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_vor(grd%nlon,nlatm2)) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_vor,work) + deallocate(grid_vor) + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=5 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Divergence + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_div(grd%nlon,nlatm2) ) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_div,work) + deallocate(grid_div) + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + end do + endif ! if ( vordivflag .or. .not. uvflag ) + + if ( uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=6 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + icount=icount+1 + iflag(icount)=7 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! V + call read_vardata(filges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(filges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + ! Note work_v and work are switched because output must be in work. + call general_filluv_ns(grd,slons,clons,grid,grid_v,work_v,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + endif ! if ( uvflag ) + + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=8 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! Specific humidity + call read_vardata(filges, 'spfh', rwork3d0, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=9 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'o3mr', rwork3d0, nslice=kr, slicedim=3) + ! Ozone mixing ratio + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + do k=1,nlevs + icount=icount+1 + iflag(icount)=10 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'clwmr', rwork3d0, nslice=kr, slicedim=3) + ! Cloud liquid water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=11 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'icmr', rwork3d0, nslice=kr, slicedim=3) + ! Cloud ice water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=12 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'rwmr', rwork3d0, nslice=kr, slicedim=3) + ! Rain water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=13 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'snmr', rwork3d0, nslice=kr, slicedim=3) + ! Snow water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=14 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'grle', rwork3d0, nslice=kr, slicedim=3) + ! Graupel mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs ! do k=1,nlevs ! icount=icount+1 @@ -3351,7 +3575,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z ! kr = levs+1-k ! netcdf is top to bottom, need to flip ! ! if (mype==mype_use(icount)) then -! call read_vardata(atmges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3) +! call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3) ! ! Cloud amount (cloud fraction). ! if ( diff_res ) then ! grid_b=rwork3d0(:,:,1) @@ -3375,6 +3599,87 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z ! endif ! enddo ! do k=1,nlevs + else ! read_2m + + icount=icount+1 + iflag(icount)=2 + + ! 2m temperature from sfc file + if (mype==mype_use(icount)) then + call read_vardata(filges, 'tmp2m', rwork2d) + + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + icount=icount + 1 + iflag(icount)=3 + + ! 2m humidity from sfc file + if (mype==mype_use(icount)) then + call read_vardata(filges, 'spfh2m', rwork2d) + + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + icount=icount + 1 + iflag(icount)=4 + + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(filges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + + ! not necessarily using all assigned tasks (fewer vars), so below doesn't trigger. + ! todo: figure out what icm should be here. + !if ( icount == icm ) then + call general_reload_sfc(grd,g_t2m, g_q2m, g_ps, icount,iflag,work) + !endif + + endif ! read_2m + + + if ( procuse ) then if ( diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) call destroy_egrid2agrid(p_high) @@ -3382,9 +3687,10 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z deallocate(rwork3d1,rwork3d0,clons,slons) deallocate(rwork2d) deallocate(grid,grid_v) - call close_dataset(atmges) + call close_dataset(filges) endif - deallocate(work, work_v) + deallocate(work) + if (allocated(work_v)) deallocate(work_v) ! Convert dry temperature to virtual temperature !do k=1,grd%nsig @@ -3396,27 +3702,29 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z !enddo ! Load u->div and v->vor slot when uv are used instead - if ( uvflag ) then - call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) then - ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - endif - else ! in this case, overload: return u/v in sf/vp slot - call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_u - call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) - if ( ier == 0 ) ptr3d=g_v - endif - if (zflag) then + if ( .not. read_2m ) then + if ( uvflag ) then + call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif !read_2m + if (read_z) then call gsi_bundlegetpointer(gfs_bundle,'z' ,ptr2d,ier) if ( ier == 0 ) ptr2d=g_z endif @@ -3429,7 +3737,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z if ( mype == 0 ) then write(6,700) lonb,latb,nlevs,grd%nlon,nlatm2,& fhour,odate,trim(filename) -700 format('GENERAL_READ_GFSATM_NC: read lonb,latb,levs=',& +700 format('GENERAL_READ_GFSATM_ALLHYDRO_NC: read lonb,latb,levs=',& 3i6,', scatter nlon,nlat=',2i6,', hour=',f6.1,', idate=',4i5,1x,a) endif diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index ed4bbb1e55..cf885c2b64 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -99,7 +99,7 @@ module gsimod factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& bcoption,diurnalbc,print_diag_pcg,tsensible,diag_precon,step_start,pseudo_q2,& - clip_supersaturation,cnvw_option + clip_supersaturation,cnvw_option,hofx_2m_sfcfile use state_vectors, only: init_anasv,final_anasv use control_vectors, only: init_anacv,final_anacv,nrf,nvars,nrf_3d,cvars3d,cvars2d,& nrf_var,lcalc_gfdl_cfrac,incvars_to_zero,incvars_zero_strat,incvars_efold @@ -1047,7 +1047,7 @@ module gsimod ! l_foreaft_thin - separate TDR fore/aft scan for thinning namelist/obs_input/dmesh,time_window_max,time_window_rad, & - ext_sonde,l_foreaft_thin + ext_sonde,l_foreaft_thin,hofx_2m_sfcfile ! SINGLEOB_TEST (one observation test case setup): ! maginnov - magnitude of innovation for one ob diff --git a/src/gsi/guess_grids.F90 b/src/gsi/guess_grids.F90 index b8cb188d01..bf493a0628 100644 --- a/src/gsi/guess_grids.F90 +++ b/src/gsi/guess_grids.F90 @@ -981,26 +981,25 @@ subroutine create_gesfinfo hrdifsig_all(nfldsig_all), & stat=istatus) if (istatus/=0) & - write(6,*)'CREATE_GESFINFO(hrdifsig,..): allocate error, istatus=',& - istatus + call die('CREATE_GESFINFO', '(hrdifsig,..): allocate error, istatus=', istatus) if(nfldsfc>0) allocate(hrdifsfc(nfldsfc),ifilesfc(nfldsfc), & hrdifsfc_all(nfldsfc_all), & stat=istatus) if (istatus/=0) & - write(6,*)'CREATE_GESFINFO(hrdifsfc,..): allocate error, istatus=',& - istatus + call die('CREATE_GESFINFO', '(hrdifsfc,..): allocate error, istatus=',& + istatus) if(nfldnst>0) allocate(hrdifnst(nfldnst),ifilenst(nfldnst), & hrdifnst_all(nfldnst_all), & - stat=istatus) + stat=istatus) if (istatus/=0) & - write(6,*)'CREATE_GESFINFO(hrdifnst,..): allocate error, istatus=',& - istatus + call die('CREATE_GESFINFO', '(hrdifnst,..): allocate error, istatus=',& + istatus) if(nfldnst>0) allocate(hrdifaer(nfldaer),ifileaer(nfldaer), & hrdifaer_all(nfldaer_all), & - stat=istatus) + stat=istatus) if (istatus/=0) & - write(6,*)'CREATE_GESFINFO(hrdifaer,..): allocate error, istatus=',& - istatus + call die('CREATE_GESFINFO', '(hrdifaer,..): allocate error, istatus=',& + istatus) #endif /* HAVE_ESMF */ return @@ -1044,16 +1043,16 @@ subroutine destroy_gesfinfo #ifndef HAVE_ESMF if(nfldsig>0) deallocate(hrdifsig,ifilesig,hrdifsig_all,stat=istatus) if (istatus/=0) & - write(6,*)'DESTROY_GESFINFO: deallocate error, istatus=',istatus + call die('DESTROY_GESFINFO', 'deallocate error, istatus=',istatus) if(nfldsfc>0) deallocate(hrdifsfc,ifilesfc,hrdifsfc_all,stat=istatus) if (istatus/=0) & - write(6,*)'DESTROY_GESFINFO: deallocate error, istatus=',istatus + call die('DESTROY_GESFINFO', 'deallocate error, istatus=',istatus) if(nfldnst>0) deallocate(hrdifnst,ifilenst,hrdifnst_all,stat=istatus) if (istatus/=0) & - write(6,*)'DESTROY_GESFINFO: deallocate error, istatus=',istatus + call die('DESTROY_GESFINFO', 'deallocate error, istatus=',istatus) if(nfldnst>0) deallocate(hrdifaer,ifileaer,hrdifaer_all,stat=istatus) if (istatus/=0) & - write(6,*)'DESTROY_GESFINFO: deallocate error, istatus=',istatus + call die('DESTROY_GESFINFO', 'deallocate error, istatus=',istatus) nfldsfc_all=0 nfldnst_all=0 diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index 1b92ad2e94..616f835218 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -136,10 +136,12 @@ module jfunc public :: pseudo_q2 public :: varq public :: cnvw_option + public :: hofx_2m_sfcfile logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option logical pseudo_q2,limitqobs + logical hofx_2m_sfcfile logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -249,6 +251,9 @@ subroutine init_jfunc ! option for including convective clouds in the all-sky assimilation cnvw_option=.false. +! option to calculate hofx for T2m and q2m by interpolating from 2m vars in sfc file + hofx_2m_sfcfile=.false. + return end subroutine init_jfunc diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index e1255b773c..41e8f33e03 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -105,6 +105,7 @@ subroutine read_ ! ! program history log: ! 2019-09-24 Martin - create routine based on read_nems +! 2022-03-23 Draper - add option to include T2m and q2m in MetGuess ! ! input argument list: ! @@ -129,6 +130,7 @@ subroutine read_ use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info use mpimod, only: npe,mype use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound + use jfunc, only: hofx_2m_sfcfile use gridmod, only: fv3_full_hydro implicit none @@ -141,6 +143,8 @@ subroutine read_ real(r_kind),pointer,dimension(:,: ):: ges_ps_it =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_z_it =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_t2m_it =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_q2m_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_u_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_v_it =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_div_it =>NULL() @@ -164,8 +168,10 @@ subroutine read_ type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 ! integer(i_kind),parameter :: n3d=8 + integer(i_kind),parameter :: n2d_2m=4 integer(i_kind),parameter :: n3d=14 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) + character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z ', 'ps ','t2m ','q2m ' /) ! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & ! 'vor ', 'div ', & ! 'tv ', 'q ', & @@ -189,8 +195,11 @@ subroutine read_ ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) - - call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) + if (hofx_2m_sfcfile) then + call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d_with2m,names3d=vars3d) + else + call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) + endif if(istatus/=0) then write(6,*) myname_,': trouble creating atm_bundle' call stop2(999) @@ -198,9 +207,15 @@ subroutine read_ do it=1,nfldsig - write(filename,'(''sigf'',i2.2)') ifilesig(it) - ! Read background fields into bundle + if (hofx_2m_sfcfile) then + if (mype==0) write(*,*) 'calling general_read_gfsatm_nc for 2m data', it + write(filename,'(''sfcf'',i2.2)') ifilesig(it) + call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& + atm_bundle,.true.,istatus) + if (mype==0) write(*,*) 'done with general_read_gfsatm_nc for 2m data', it + end if + write(filename,'(''sigf'',i2.2)') ifilesig(it) if (fv3_full_hydro) then if (mype==0) write(*,*) 'calling general_read_gfsatm_allhydro_nc', it call general_read_gfsatm_allhydro_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& @@ -273,6 +288,16 @@ subroutine set_guess_ call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,istatus) if(istatus==0) ges_z_it = ptr2d endif + call gsi_bundlegetpointer (atm_bundle,'t2m',ptr2d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'t2m' ,ges_t2m_it ,istatus) + if(istatus==0) ges_t2m_it = ptr2d + endif + call gsi_bundlegetpointer (atm_bundle,'q2m',ptr2d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q2m' ,ges_q2m_it ,istatus) + if(istatus==0) ges_q2m_it = ptr2d + endif call gsi_bundlegetpointer (atm_bundle,'u',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,istatus) diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index e8d7ec756f..771d2b1443 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -148,6 +148,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! ! 2020-05-04 wu - no rotate_wind for fv3_regional ! 2020-09-05 CAPS(C. Tong) - add flag for new vadwind obs to assimilate around the analysis time only +! 2023-03-23 draper - add code for processing T2m and q2m for global system ! input argument list: ! infile - unit from which to read BUFR data @@ -212,7 +213,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& use hilbertcurve,only: init_hilbertcurve, accum_hilbertcurve, & apply_hilbertcurve,destroy_hilbertcurve use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error - use jfunc, only: tsensible + use jfunc, only: tsensible, hofx_2m_sfcfile use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_deter @@ -263,7 +264,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& logical tob,qob,uvob,spdob,sstob,pwob,psob,gustob,visob,tdob,mxtmob,mitmob,pmob,howvob,cldchob logical metarcldobs,goesctpobs,tcamtob,lcbasob logical outside,driftl,convobs,inflate_error - logical sfctype + logical sfctype, global_2m_land logical luse,ithinp,windcorr logical patch_fog logical aircraftset,aircraftobs,aircraftobst,aircrafttype @@ -1614,10 +1615,15 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& pmq(k)=idnint(qcmark(8,k)) end do +! 181, 183, 187, and 188 are the screen-level obs over land + global_2m_land = ( (kx==181 .or. kx==183 .or. kx==188 .or. kx==188 ) .and. hofx_2m_sfcfile ) + ! If temperature ob, extract information regarding virtual ! versus sensible temperature if(tob) then - if (.not. twodvar_regional .or. .not.tsensible) then + ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode + if ( (.not. tsensible) .and. .not. (twodvar_regional .or. global_2m_land) ) then + do k=1,levs tvflg(k)=one ! initialize as sensible do j=1,20 @@ -1914,6 +1920,26 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Missing Values ==> Cycling! In this case for howv only. #ww3 if (howvob .and. owave(1,k) > r0_1_bmiss) cycle LOOP_K_LEVS +! Over-ride QM=9 and hard-wire errors for land obs and hofx_sfc_file option +! Can be deleted once prepbufr processing updated. + if ( global_2m_land ) then + if (tob .and. qm==9 ) then + pqm(k)=2 ! otherwise, type 183 will be discarded. + qm=2 + tqm(k)=2 + if (kx==187) obserr(3,k)=2.2 + if (kx==181) obserr(3,k)=1.5 + if (kx==183) obserr(3,k)=2.6 + endif + if (qob .and. qm == 9 ) then + qm = 2 + ! qob err specified as fraction of qsat, multiplied by 10. + if (kx==187) obserr(2,k)=1.0 + if (kx==181) obserr(2,k)=1.0 + if (kx==183) obserr(2,k)=1.0 + endif + + endif ! Set usage variable usage = zero if(icuse(nc) <= 0)usage=100._r_kind @@ -1957,6 +1983,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if(obsdat(12,k) > 32.2_r_kind) usage=118._r_kind ! > 90F endif endif + ! to-do: should we add qob checks from above for landsfctype too? if ((kx>129.and.kx<140).or.(kx>229.and.kx<240) ) then call get_aircraft_usagerj(kx,obstype,c_station_id,usage) diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index cc06d169a7..aa557b72c2 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -111,6 +111,8 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! information in diagonostic file, which is used ! in offline observation quality control program (AutoObsQC) ! for 3D-RTMA (if l_obsprvdiag is true). +! 2023-03-09 Draper added option to interpolate screen-level q from model 2m output. +! (hofx_2m_sfcfile) ! ! ! input argument list: @@ -160,7 +162,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use constants, only: huge_single,wgtlim,three use constants, only: tiny_r_kind,five,half,two,huge_r_kind,r0_01 use qcmod, only: npres_print,ptopq,pbotq,dfact,dfact1,njqc,vqc,nvqc - use jfunc, only: jiter,last,jiterstart,miter,superfact,limitqobs + use jfunc, only: jiter,last,jiterstart,miter,superfact,limitqobs,hofx_2m_sfcfile use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype use convinfo, only: ibeta,ikapa use convinfo, only: icsubtype @@ -217,7 +219,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Declare local variables real(r_double) rstation_id - real(r_kind) qob,qges,qsges,q2mges,q2mges_water + real(r_kind) qob,qges,qsges,q2mges,q2mges_water,qsges_o real(r_kind) ratio_errors,dlat,dlon,dtime,dpres,rmaxerr,error real(r_kind) rsig,dprpx,rlow,rhgh,presq,tfact,ramp real(r_kind) psges,sfcchk,ddiff,errorx @@ -231,6 +233,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_kind),dimension(nobs):: dup real(r_kind),dimension(lat2,lon2,nsig,nfldsig):: qg real(r_kind),dimension(lat2,lon2,nfldsig):: qg2m + real(r_kind),dimension(lat2,lon2,nfldsig):: qg2m_o real(r_kind),dimension(nsig):: prsltmp real(r_kind),dimension(34):: ptablq real(r_single),allocatable,dimension(:,:)::rdiagbuf @@ -277,10 +280,22 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q real(r_kind),allocatable,dimension(:,:,: ) :: ges_q2m + real(r_kind),allocatable,dimension(:,:,: ) :: ges_t2m logical:: l_pbl_pseudo_itype integer(i_kind):: ich0 type(obsLList),pointer,dimension(:):: qhead + + logical :: landsfctype + + real(r_kind) :: delta_z, lapse_error, q_delta_terrain + real(r_kind), parameter :: T_lapse = -0.0045 ! standard lapse rate, K/m +! use 4.5 K/km, in place of more standard 6.5 K/km, following +! https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019EA000984 +! lapse_error_frac around 0.5 ~ 2K/km, from Figure 2 of above. + real(r_kind), parameter :: lapse_error_frac = 0.5 ! inflation factor for obs error when vertically interpolating + real(r_kind), parameter :: max_delta_z = 300. ! max. vertical mismatch allowed (later: relax this) + qhead => obsLL(:) save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf @@ -359,8 +374,11 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav hr_offset=min_offset/60.0_r_kind dup=one do k=1,nobs + ikx=nint(data(ikxx,k)) + itype=ictype(ikx) + landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188 ) do l=k+1,nobs - if (twodvar_regional) then + if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & data(ilon,k) == data(ilon,l) .and. & data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & @@ -425,9 +443,15 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ice=.false. ! get larger (in rh) q obs error for mixed and ice phases iderivative=0 + + ! calculate qsat and 2m qsat do jj=1,nfldsig - call genqsat(qg(1,1,1,jj),ges_tsen(1,1,1,jj),ges_prsl(1,1,1,jj),lat2,lon2,nsig,ice,iderivative) - qg2m(:,:,jj)=qg(:,:,1,jj) + call genqsat(qg(:,:,:,jj),ges_tsen(:,:,:,jj),ges_prsl(:,:,:,jj),lat2,lon2,nsig,ice,iderivative) + if (i_use_2mq4b > 0) then ! use lowest model level + qg2m(:,:,jj)=qg(:,:,1,jj) + elseif ( hofx_2m_sfcfile ) then ! calculate from 2m model output + call genqsat(qg2m(:,:,jj),ges_t2m(:,:,jj),ges_ps(:,:,jj),lat2,lon2,1,ice,iderivative) + endif end do @@ -440,10 +464,10 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call dtime_check(dtime, in_curbin, in_anybin) if(.not.in_anybin) cycle + landsfctype =( itype==181 .or. itype==183 .or. itype==187 ) ! Flag static conditions to create PBL_pseudo_surfobsq obs. - l_pbl_pseudo_itype = l_pbl_pseudo_surfobsq .and. & - ( itype==181 .or. itype==183 .or.itype==187 ) + l_pbl_pseudo_itype = l_pbl_pseudo_surfobsq .and. landsfctype if(in_curbin) then ! Convert obs lats and lons to grid coordinates @@ -509,24 +533,28 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav presq=r10*exp(dpres) itype=ictype(ikx) dprpx=zero - if(((itype > 179 .and. itype < 190) .or. itype == 199) & + + if ( hofx_2m_sfcfile .and. landsfctype) then + dpres = one ! put obs on surface + else + if(((itype > 179 .and. itype < 190) .or. itype == 199) & .and. .not.twodvar_regional)then - dprpx=abs(one-exp(dpres-log(psges)))*r10 - end if + dprpx=abs(one-exp(dpres-log(psges)))*r10 + endif ! Put obs pressure in correct units to get grid coord. number - call grdcrd1(dpres,prsltmp(1),nsig,-1) + call grdcrd1(dpres,prsltmp(1),nsig,-1) ! Get approximate k value of surface by using surface pressure - sfcchk=log(psges) - call grdcrd1(sfcchk,prsltmp(1),nsig,-1) + sfcchk=log(psges) + call grdcrd1(sfcchk,prsltmp(1),nsig,-1) ! Check to see if observations is above the top of the model (regional mode) - if( dpres>=nsig+1)dprpx=1.e6_r_kind - if((itype > 179 .and. itype < 186) .or. itype == 199) dpres=one + if( dpres>=nsig+1)dprpx=1.e6_r_kind + if((itype > 179 .and. itype < 186) .or. itype == 199) dpres=one + + endif -! Scale errors by guess saturation q - qob = data(iqob,i) if(limitqobs) then call tintrp31(ges_qsat,qsges,dlat,dlon,dpres,dtime,hrdifsig,& @@ -534,11 +562,13 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav qob=min(qob,superfact*qsges) end if +! get qsges, to be used to scale the obs error call tintrp31(qg,qsges,dlat,dlon,dpres,dtime,hrdifsig,& mype,nfldsig) -! Interpolate 2-m qs to obs locations/times - if((i_use_2mq4b > 0) .and. ((itype > 179 .and. itype < 190) .or. itype == 199) & - .and. .not.twodvar_regional)then + +! overwrite qsges with 2-m qs if sfc obs scheme + if( ( (i_use_2mq4b > 0) .and. ((itype > 179 .and. itype < 190) .or. itype == 199) & + .and. .not.twodvar_regional) .or. (hofx_2m_sfcfile .and. landsfctype) )then call tintrp2a11(qg2m,qsges,dlat,dlon,dtime,hrdifsig,mype,nfldsig) endif @@ -549,10 +579,36 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav rmaxerr=max(small2,rmaxerr) errorx =(data(ier,i)+dprpx)*qsges -! Interpolate guess moisture to observation location and time - call tintrp31(ges_q,qges,dlat,dlon,dpres,dtime, & - hrdifsig,mype,nfldsig) - +! qges: Interpolate guess moisture to observation location and time + + if (.not. ( hofx_2m_sfcfile .and. landsfctype) ) then + call tintrp31(ges_q,qges,dlat,dlon,dpres,dtime, & + hrdifsig,mype,nfldsig) + else + ! only use land locations + if (int(data(idomsfc,i)) .NE. 1 ) muse(i) = .false. + + call tintrp2a11(ges_q2m,qges,dlat,dlon,dtime,hrdifsig,mype,nfldsig) + + ! terrain correction: assume RH_zo = RH_zm, and correct T with + ! same lapse rate as used for T2m terrain correction + + delta_z = data(istnelv,i) - data(izz,i) ! obs -model + + do jj=1,nfldsig + ! qsat in model at height of obs + call genqsat(qg2m_o(:,:,jj),ges_t2m(:,:,jj)+delta_z*T_lapse,ges_ps(:,:,jj),lat2,lon2,1,ice,iderivative) + enddo + + call tintrp2a11(qg2m_o,qsges_o,dlat,dlon,dtime,hrdifsig,mype,nfldsig) + q_delta_terrain = (qsges/qsges_o - 1)*qob + qob = qob * ( qsges/qsges_o) + + !update the station elevation + data(istnelv,i) = data(izz,i) + + endif + ddiff=qob-qges ! Setup dynamic ob error specification for aircraft recon in hurricanes @@ -572,18 +628,22 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif errorx =max(small1,errorx) - ! Adjust observation error to reflect the size of the residual. ! If extrapolation occurred, then further adjust error according to ! amount of extrapolation. - rlow=max(sfcchk-dpres,zero) + if (.not. (hofx_2m_sfcfile .and. landsfctype) ) then + rlow=max(sfcchk-dpres,zero) ! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] - if(l_sfcobserror_ramp_q) then - ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind)*0.001_r_kind + if(l_sfcobserror_ramp_q) then + ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind)*0.001_r_kind + else + ramp=rlow + endif else - ramp=rlow + rlow = zero + ramp = zero endif rhgh=max(dpres-r0_001-rsig,zero) @@ -594,7 +654,20 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(rhgh/=zero) awork(3) = awork(3) + one end if - ratio_errors=error*qsges/(errorx+1.0e6_r_kind*rhgh+r8*ramp) +! inflate error for uncertainty in the terrain adjustment + lapse_error = 0. + if ( hofx_2m_sfcfile .and. landsfctype) then + if (abs(delta_z)max_delta_z do not assim. + ! inflate obs error to account for error in lapse_rate + ! also include some representativity error here (assuming + ! delta_z ~ heterogeneity) + lapse_error = abs(lapse_error_frac*q_delta_terrain) + else + muse(i)=.false. + endif + endif + + ratio_errors=error*qsges/(errorx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) ! Check to see if observations is above the top of the model (regional mode) if (dpres > rsig) ratio_errors=zero @@ -618,7 +691,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dhx_dx%val(2) = delz ! weight for iz+1's level endif -! Interpolate 2-m q to obs locations/times +! i_use_2mq4b: Interpolate 2-m q to obs locations/times if(i_use_2mq4b>0 .and. itype > 179 .and. itype < 190 .and. .not.twodvar_regional)then if(i_coastline==2 .or. i_coastline==3) then @@ -643,7 +716,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call stop2(100) endif ddiff=qob-qges - endif + endif ! i_use_2mq4b ! If requested, setup for single obs test. @@ -943,7 +1016,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav my_head => null() ENDDO - endif ! 181,183,187 + endif ! l_pbl_pseudo_itype !!!!!!!!!!!!!!!!!! PBL pseudo surface obs !!!!!!!!!!!!!!!!!!!!!!! ! End of loop over observations @@ -1025,7 +1098,7 @@ subroutine init_vars_ call stop2(999) endif ! get q2m ... - if (i_use_2mq4b>0) then + if (i_use_2mq4b>0 .or. hofx_2m_sfcfile) then varname='q2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) if (istatus==0) then @@ -1044,6 +1117,25 @@ subroutine init_vars_ call stop2(999) endif endif ! i_use_2mq4b + if (hofx_2m_sfcfile) then + varname='t2m' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_t2m))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_t2m(size(rank2,1),size(rank2,2),nfldsig)) + ges_t2m(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_t2m(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + endif ! hofx_2m_sfcfile ! get q ... varname='q' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) @@ -1272,8 +1364,10 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) call nc_diag_metadata_to_single("Latitude", data(ilate,i) ) call nc_diag_metadata_to_single("Longitude", data(ilone,i) ) +! this is the obs height after being interpolated to the model (=model height) call nc_diag_metadata_to_single("Station_Elevation",data(istnelv,i) ) call nc_diag_metadata_to_single("Pressure", presq ) +! this is the original obs height (= stn elevation, before being interpolated) call nc_diag_metadata_to_single("Height", data(iobshgt,i) ) call nc_diag_metadata_to_single("Time", dtime,time_offset,'-' ) call nc_diag_metadata_to_single("Prep_QC_Mark", data(iqc,i) ) @@ -1392,6 +1486,7 @@ end subroutine contents_netcdf_diagp_ subroutine final_vars_ if(allocated(ges_q2m)) deallocate(ges_q2m) + if(allocated(ges_t2m)) deallocate(ges_t2m) if(allocated(ges_q )) deallocate(ges_q ) if(allocated(ges_ps)) deallocate(ges_ps) end subroutine final_vars_ diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index f33dcf69df..5467a6dec9 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -54,11 +54,11 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use gridmod, only: nsig,twodvar_regional,regional use gridmod, only: get_ijk,pt_ll - use jfunc, only: jiter,last,jiterstart,miter + use jfunc, only: jiter,last,jiterstart,miter,hofx_2m_sfcfile use guess_grids, only: nfldsig, hrdifsig,ges_lnprsl,& geop_hgtl,ges_tsen,pbl_height - use state_vectors, only: svars3d, levels + use state_vectors, only: svars3d, levels, ns3d, svars2d use constants, only: zero, one, four,t0c,rd_over_cp,three,rd_over_cp_mass,ten use constants, only: tiny_r_kind,half,two @@ -228,6 +228,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! for 3D-RTMA (if l_obsprvdiag is true). ! 2022-03-15 Hu change all th2 to t2m to indicate that 2m temperature ! is sensible instead of potentionl temperature +! 2023-03-21 Draper added option to interpolate screen-level T from model 2m output. +! (hofx_2m_sfcfile) ! ! !REMARKS: ! language: f90 @@ -309,7 +311,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav logical,dimension(nobs):: luse,muse integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID - logical sfctype + logical sfctype, landsfctype logical iqtflg logical aircraftobst logical duplogic @@ -342,6 +344,18 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav integer(i_kind):: ich0 type(obsLList),pointer,dimension(:):: thead + + real(r_kind) :: delta_z, lapse_error + real(r_kind), parameter :: T_lapse = -0.0045 ! standard lapse rate, K/m +! use 4.5 K/km, in place of more standard 6.5 K/km, following +! https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019EA000984 +! lapse_error_frac around 0.5 ~ 2K/km, from Figure 2 of above. + real(r_kind), parameter :: lapse_error_frac = 0.5 ! inflation factor for obs error when vertically interpolating + real(r_kind), parameter :: max_delta_z = 300. ! max. vertical mismatch allowed + +! CSD - move this to where the namelists are read in. + if (i_use_2mt4b>0) hofx_2m_sfcfile=.false. + thead => obsLL(:) save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf @@ -432,8 +446,11 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav hr_offset=min_offset/60.0_r_kind dup=one do k=1,nobs + ikx=nint(data(ikxx,k)) + itype=ictype(ikx) + landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188 ) do l=k+1,nobs - if (twodvar_regional) then + if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & data(ilon,k) == data(ilon,l) .and. & data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & @@ -465,6 +482,10 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end do ! Run a buddy-check +! Note: buddy check crashes for hofx_2m_sfcfile option. +! Ccurrent params have buddy radius of 108 km, max diff of 8 K. +! The gross error check removes O-F > 7., so this is probably removing +! most obs that fail the buddy check already if (twodvar_regional .and. buddycheck_t) call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) ! If requested, save select data for output to diagnostic file @@ -521,6 +542,10 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav rstation_id = data(id,i) prest=r10*exp(dpres) ! in mb sfctype=(itype>179.and.itype<190).or.(itype>=192.and.itype<=199) +! hofx_2m_sfcfile option to calculate hofx from 2m model output (rather than LML) +! is restricted to landsfctype only. GDAS assimilates 180 and 182 over ocean, +! should we also use 2m model output for the over-ocean obs? + landsfctype =( itype==181 .or. itype==183 .or. itype==187 ) iqtflg=nint(data(iqt,i)) == 0 var_jb=data(ijb,i) @@ -654,17 +679,22 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) - drpx=zero - if(sfctype .and. .not.twodvar_regional) then - drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c - end if + drpx = zero + if ( hofx_2m_sfcfile .and. landsfctype) then + dpres = one ! put obs at surface + else + if(sfctype .and. .not.twodvar_regional) then + drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c + end if -! Put obs pressure in correct units to get grid coord. number - call grdcrd1(dpres,prsltmp(1),nsig,-1) +! Put obs pressure in correct units to get grid coord. number + call grdcrd1(dpres,prsltmp(1),nsig,-1) + endif ! Implementation of forward model ---------- - if(sfctype.and.sfcmodel) then +! SCENARIO 1: If obs is sfctype, and sfcmodel is requested. Outdated. + if(sfctype .and. sfcmodel) then tgges=data(iskint,i) roges=data(isfcr,i) @@ -694,8 +724,47 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges +! SCENARIO 2: obs is sfctype, and hofx_2m_sfcfile scheme is on. +! 2m forecast has been read from the sfc guess files + elseif (landsfctype .and. hofx_2m_sfcfile ) then + +! mask: 0 - sea, 1 - land, 2-ice, >= 3 mixed +! for now, use only pure land + if (int(data(idomsfc,i)) .NE. 1 ) muse(i) = .false. + + call tintrp2a11(ges_t2m,tges2m,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + +! correct obs to model terrain height using a standard lapse rate. +! Later: look into updating with lapse-rate from the model (similar to gsd_terrain_match) + + delta_z = data(izz,i) - data(istnelv,i) + tob = tob + delta_z*T_lapse + !update the station elevation + data(istnelv,i) = data(izz,i) + + if(save_jacobian) then + t_ind = getindex(svars2d, 't2m') + if (t_ind < 0) then + print *, 'Error: no variable t2m in state vector.Exiting.' + call stop2(1300) + endif + dhx_dx%st_ind(1) = sum(levels(1:ns3d)) + t_ind + dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t_ind + dhx_dx%val(1) = one + dhx_dx%val(2) = zero ! in this case, there is no vertical interp + ! and nnz (=dim(dhx_dx%val)) should be one, + ! but nnz is a file attribute, so need to use + ! same value as for vertical profile obs. Get + ! around this by setting val(2) to zero. + endif + +! SCENARIO 3: obs is sfctype, and neither sfcmodel nor hofx_2m_sfcfile is chosen +! .or. obs is not sfctype. Interpoate hofx from model levels. else + if(iqtflg)then +! SCENARIO 3a: obs is a virtual temp. ! Interpolate guess tv to observation location and time call tintrp31(ges_tv,tges,dlat,dlon,dpres,dtime, & hrdifsig,mype,nfldsig) @@ -717,6 +786,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dhx_dx%val(2) = delz ! weight for iz+1's level endif else +! SCENARIO 3b: obs is a sensible temp. ! Interpolate guess tsen to observation location and time call tintrp31(ges_tsen,tges,dlat,dlon,dpres,dtime, & hrdifsig,mype,nfldsig) @@ -739,6 +809,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif end if + +! SCENARIO 4: obs is sfctype, and i_use_2mt4b flag is on (turns on regional sfc DA) if(i_use_2mt4b>0 .and. sfctype) then if(i_coastline==1 .or. i_coastline==3) then @@ -773,17 +845,23 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call grdcrd1(sfcchk,prsltmp(1),nsig,-1) ! Check to see if observations is above the top of the model (regional mode) - if(sfctype)then + if(sfctype .and. .not. (hofx_2m_sfcfile .and. landsfctype) )then if(abs(dpres)>four) drpx=1.0e10_r_kind pres_diff=prest-r10*psges if (twodvar_regional .and. abs(pres_diff)>=r1000) drpx=1.0e10_r_kind end if - rlow=max(sfcchk-dpres,zero) -! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] - if(l_sfcobserror_ramp_t) then - ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind) + + if (.not. (hofx_2m_sfcfile .and. landsfctype) ) then + rlow=max(sfcchk-dpres,zero) +! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] + if(l_sfcobserror_ramp_t) then + ramp=min(max(((rlow-1.0_r_kind)/(15.0_r_kind-1.0_r_kind)),0.0_r_kind),1.0_r_kind) + else + ramp=rlow + endif else - ramp=rlow + rlow = zero + ramp = zero endif rhgh=max(zero,dpres-rsigp-r0_001) @@ -795,12 +873,26 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(rlow/=zero) awork(2) = awork(2) + one if(rhgh/=zero) awork(3) = awork(3) + one end if - - ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp) + +! inflate error for uncertainty in the terrain adjustment + lapse_error = 0. + if ( hofx_2m_sfcfile .and. landsfctype) then + if (abs(delta_z)max_delta_z do not assim. + ! inflate obs error to account for error in lapse_rate + ! also include some representativity error here (assuming + ! delta_z ~ heterogeneity) + lapse_error = abs(lapse_error_frac*T_lapse*delta_z) + else + muse(i)=.false. + endif + endif + + ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) ! Compute innovation - if(i_use_2mt4b>0 .and. sfctype) then + if( (sfctype .and. i_use_2mt4b>0) .or. (hofx_2m_sfcfile .and. landsfctype) ) then ddiff = tob-tges2m + if (hofx_2m_sfcfile) tges=tges2m else ddiff = tob-tges endif @@ -1411,7 +1503,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif - if(i_use_2mt4b>0) then + if(i_use_2mt4b>0 .or. hofx_2m_sfcfile) then ! get t2m ... varname='t2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1430,6 +1522,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif + ! get q2m ... varname='q2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1676,8 +1769,10 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) call nc_diag_metadata_to_single("Latitude",data(ilate,i)) call nc_diag_metadata_to_single("Longitude",data(ilone,i)) +! this is the obs height after being interpolated to the model (=model height) call nc_diag_metadata_to_single("Station_Elevation",data(istnelv,i)) call nc_diag_metadata_to_single("Pressure",prest) +! this is the original obs height (= stn elevation, before being interpolated) call nc_diag_metadata_to_single("Height",data(iobshgt,i)) call nc_diag_metadata_to_single("Time",dtime,time_offset,'-') call nc_diag_metadata_to_single("Prep_QC_Mark",data(iqc,i)) @@ -1693,7 +1788,11 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata_to_single("Errinv_Input",errinv_input ) call nc_diag_metadata_to_single("Errinv_Adjust",errinv_adjst ) call nc_diag_metadata_to_single("Errinv_Final",errinv_final ) - call nc_diag_metadata_to_single("Observation",data(itob,i)) + if (hofx_2m_sfcfile ) then + call nc_diag_metadata_to_single("Observation", tob ) + else + call nc_diag_metadata_to_single("Observation", data(itob,i) ) + endif call nc_diag_metadata_to_single("Obs_Minus_Forecast_adjusted",ddiff ) call nc_diag_metadata_to_single("Obs_Minus_Forecast_unadjusted",tob,tges,'-') diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index e5a0f64245..1885542a1b 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,hofx_2m_sfcfile use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -458,7 +458,7 @@ subroutine update_guess(sval,sbias) endif call gsd_update_soil_tq(tinc_1st,is_t,qinc_1st,is_q,it) endif ! l_gsd_soilTQ_nudge - if (i_use_2mt4b > 0 .and. is_t>0) then + if ( (i_use_2mt4b > 0.or. hofx_2m_sfcfile) .and. is_t>0) then do j=1,lon2 do i=1,lat2 tinc_1st(i,j)=p_tv(i,j,1) @@ -466,7 +466,7 @@ subroutine update_guess(sval,sbias) end do call gsd_update_t2m(tinc_1st,it) endif ! l_gsd_t2m_adjust - if (i_use_2mq4b > 0 .and. is_q>0) then + if ( (i_use_2mq4b > 0.or. hofx_2m_sfcfile) .and. is_q>0) then do j=1,lon2 do i=1,lat2 qinc_1st(i,j)=p_q(i,j,1)