From 3e0143da09eec387930b41ff01e6717fcc30c8d6 Mon Sep 17 00:00:00 2001 From: Ricardo Todling Date: Thu, 24 Oct 2024 11:58:41 -0400 Subject: [PATCH 1/2] Allow for consistency of NCEP-sfc fields with those used in JEDI --- GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt | 1 + .../GSI_GridComp/GSI_GridCompMod.F90 | 139 +++++-- GEOSaana_GridComp/GSI_GridComp/analyzer | 34 +- .../GSI_GridComp/m_nc_ncepsfc.f90 | 349 ++++++++++++++++++ 4 files changed, 474 insertions(+), 49 deletions(-) create mode 100644 GEOSaana_GridComp/GSI_GridComp/m_nc_ncepsfc.f90 diff --git a/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt b/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt index cd308b5..0424f0d 100644 --- a/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt +++ b/GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt @@ -411,6 +411,7 @@ set (SRCS_SOLVER logcldch_to_cldch.f90 logvis_to_vis.f90 m_nc_berror.f90 + m_nc_ncepsfc.f90 m_berror_stats.f90 m_berror_stats_reg.f90 m_cvgridLookup.F90 diff --git a/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 b/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 index e7a921a..95bb5f7 100644 --- a/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 +++ b/GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90 @@ -130,6 +130,7 @@ module GSI_GridCompMod ! (buffer points on ends) ak5,bk5,ck5, & ! coefficients for hybrid vertical coordinate idvc5, & ! vertical coordinate identifier + nvege_type, & ! veg-type choice ! spectral transform grid info sp_a @@ -1823,7 +1824,6 @@ subroutine Scale_Import_() where ( du002p /= MAPL_UNDEF ) du002p = du002p * KGpKG2PPBVaero ! convert from kg/kg to ppbv endwhere - print *, 'DEBUG doing the right thing ...' endif if(associated(du003p)) then where ( du003p /= MAPL_UNDEF ) @@ -2566,6 +2566,9 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) use sfcio_module use mpeu_util, only: die + use m_nc_ncepsfc, only: nc_ncepsfc_read + use m_nc_ncepsfc, only: nc_ncepsfc_vars + use m_nc_ncepsfc, only: nc_ncepsfc_vars_final implicit none integer, intent(in) :: lit ! logical index of first guess time level to copy @@ -2575,6 +2578,12 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) ! 17Mar2016 Todling Revisit handling of extra surface fields: use sfcio; ! add option to set veg-type to constant for sensitivity ! evalulation. +! 22Oct2024 Todling In an attempt to minimize differences between GSI and JEDI +! this routine provides the means to read an nc4 file with +! the three surface fields in question that are consistent +! with the data used in JEDI. Notice: this implies that +! a lat-lon file (at GSI''s resolution) has been created +! from the cubed files used in JEDI. ! !EOPI !------------------------------------------------------------------------- @@ -2594,14 +2603,19 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) type(sfcio_head):: head type(sfcio_data):: data + type(nc_ncepsfc_vars) :: sfcvars + character(len=*), parameter :: & IAm='GSI_GridCompGetNCEPsfcFromFile_' character(len=ESMF_MAXSTR) :: opt integer(i_kind) :: it,nn,iret integer(i_kind) :: iset_veg_type + real(r_single) :: rvege logical :: wrt_ncep_sfc_grd logical :: use_sfcio + logical :: use_sfcnc + logical :: bypass_interp ! start @@ -2618,12 +2632,13 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) endif call ESMF_ConfigGetAttribute( CF, opt, label ='USE_SFCIO4NCEP_SFC:', default="YES", rc = STATUS ) VERIFY_(STATUS) - use_sfcio=.true. - if ( trim(opt) == "YES" ) then - if(IamRoot) print *,trim(Iam),': using SFCIO to read NCEP surface fields' + use_sfcnc=nvege_type==20 ! wired for now + if ( use_sfcnc ) then + if(IamRoot) print *,trim(Iam),': read NCEP surface fields from nc4 file with 20-veg' else - use_sfcio=.false. - if(IamRoot) print *,trim(Iam),': NOT using SFCIO to read NCEP surface fields' + use_sfcio=.false. + if (trim(opt) == 'YES') use_sfcio=.true. + if(IamRoot) print *,trim(Iam),': read NCEP surface fields using SFCIO: ', use_sfcio endif if(lit > nfldsfc) then @@ -2640,6 +2655,7 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) ! if(IamRoot) then + bypass_interp=.false. if (wrt_ncep_sfc_grd) then call baopenwt(36,'ncepsfc.grd',iret) endif @@ -2660,11 +2676,26 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) print*, trim(Iam), ', Veg-type reset by request to: ',iset_veg_type endif else - open(34,file='ncepsfc', form='unformatted') - read (34) - read(34) yhour,igdate,glon,glat2,version - close(34) - glat=glat2+2 + if (use_sfcnc) then + yhour=0 + igdate=1776074 + call nc_ncepsfc_read('ncepsfc',sfcvars,STATUS) + VERIFY_(status) + if(sfcvars%nveg/=nvege_type) then + status=99 + VERIFY_(STATUS) + endif + glon =size(sfcvars%vfrac,1) + glat2=size(sfcvars%vfrac,2) + bypass_interp= (glat2==nlat).and.(glon==nlon) + print*, trim(Iam), ', sfc-file-handling bypassing interpolation: ', bypass_interp + else + open(34,file='ncepsfc', form='unformatted') + read (34) + read(34) yhour,igdate,glon,glat2,version + close(34) + glat=glat2+2 + endif endif print *,trim(Iam),': ncepsfc hour/date : ',yhour,' / ',igdate print *,trim(Iam),': ncepsfc fields dims, version: ',glon,' x ',glat2, ' v', version @@ -2680,9 +2711,13 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) VERIFY_(STATUS) if (.not.use_sfcio ) then - call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & - STATUS, jrec=12, fld=sfcbuf) - VERIFY_(STATUS) + if (use_sfcnc) then + vfrbuf=sfcvars%vfrac + else + call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & + STATUS, jrec=12, fld=sfcbuf) + VERIFY_(STATUS) + endif else if (associated(data%vfrac)) then sfcbuf=data%vfrac @@ -2690,12 +2725,15 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) call die(Iam,' vfrac not found in NCEP input file', 99) endif endif - call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) - call GSI_GridCompFlipLons_(sfcbufpp) - call hinterp2 ( sfcbufpp ,glon, glat , & - vfrbuf, nlon, nlat, 1, & - UNDEF_SSI_) + if (.not.bypass_interp) then + call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) + call GSI_GridCompFlipLons_(sfcbufpp) + call hinterp2 ( sfcbufpp ,glon, glat , & + vfrbuf, nlon, nlat, 1, & + UNDEF_SSI_) + endif call GSI_GridCompFlipLons_(vfrbuf) + where(vfrbuf<0.0) vfrbuf=0.0 end where @@ -2707,9 +2745,13 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) endif if (.not.use_sfcio ) then - call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & - STATUS, jrec=15, fld=sfcbuf) - VERIFY_(STATUS) + if (use_sfcnc) then + vtybuf=sfcvars%vtype + else + call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & + STATUS, jrec=15, fld=sfcbuf) + VERIFY_(STATUS) + endif else if (associated(data%vtype)) then sfcbuf=data%vtype @@ -2717,23 +2759,30 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) call die(Iam,' vtype not found in NCEP input file', 99) endif endif - call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) - call GSI_GridCompFlipLons_(sfcbufpp) - call hinterp2 ( sfcbufpp ,glon, glat , & - vtybuf, nlon, nlat, 1, & - UNDEF_SSI_) + if (.not.bypass_interp) then + call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) + call GSI_GridCompFlipLons_(sfcbufpp) + call hinterp2 ( sfcbufpp ,glon, glat , & + vtybuf, nlon, nlat, 1, & + UNDEF_SSI_) + endif call GSI_GridCompFlipLons_(vtybuf) vtybuf=real(nint(vtybuf),r_single) - where(vtybuf < 0.0_r_single) vtybuf = 0.0_r_single - where(vtybuf > 13.0_r_single) vtybuf = 13.0_r_single + rvege = float(nvege_type-1) + where(vtybuf < 0.0_r_single) vtybuf = 0.0_r_single + where(vtybuf > rvege) vtybuf = rvege if (wrt_ncep_sfc_grd) then call wryte(36,4*nlat*nlon,vtybuf) endif if (.not.use_sfcio ) then - call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & - STATUS, jrec=16, fld=sfcbuf) - VERIFY_(STATUS) + if (use_sfcnc) then + stybuf=sfcvars%stype + else + call ncep_rwsurf_ ( .false., 'ncepsfc', glat2, glon, & + STATUS, jrec=16, fld=sfcbuf) + VERIFY_(STATUS) + endif else if (associated(data%stype)) then sfcbuf=data%stype @@ -2741,11 +2790,13 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) call die(Iam,' stype not found in NCEP input file', 99) endif endif - call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) - call GSI_GridCompFlipLons_(sfcbufpp) - call hinterp2 ( sfcbufpp ,glon, glat , & - stybuf, nlon, nlat, 1, & - UNDEF_SSI_) + if (.not.bypass_interp) then + call GSI_GridCompSP2NP_(sfcbufpp,sfcbuf) + call GSI_GridCompFlipLons_(sfcbufpp) + call hinterp2 ( sfcbufpp ,glon, glat , & + stybuf, nlon, nlat, 1, & + UNDEF_SSI_) + endif call GSI_GridCompFlipLons_(stybuf) stybuf=real(nint(stybuf),r_single) where(stybuf < 0.0_r_single) stybuf = 0.0_r_single @@ -2754,6 +2805,9 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) call wryte(36,4*nlat*nlon,stybuf) call baclose(36,iret) endif + if (use_sfcnc) then + call nc_ncepsfc_vars_final(sfcvars) + endif deallocate(sfcbuf, sfcbufpp, stat=STATUS) VERIFY_(STATUS) @@ -2810,7 +2864,7 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit) #endif end subroutine GSI_GridCompGetNCEPsfcFromFile_ - + !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- @@ -3421,10 +3475,14 @@ subroutine getnexttoken1_(cf,token,myname,rc) character(len=*),intent(in ) :: myname integer ,intent(out) :: rc character(len=*), parameter :: IAm='GSI_GridComp.getnexttoken1_' + logical :: tableEnd token="" - call ESMF_ConfigNextLine(cf, rc=rc) - if(rc/=0) return ! end-of-file is expected. No error message is produced. + call ESMF_ConfigNextLine(cf, tableEnd=tableEnd, rc=rc) + if(tableEnd) then ! end-of-file is expected. No error message is produced. + rc=-1 ! end-of-table is expected. No error message is produced + return + endif call ESMF_ConfigGetAttribute(cf, token, rc=rc) if(rc/=0) then @@ -3432,7 +3490,6 @@ subroutine getnexttoken1_(cf,token,myname,rc) call MAPL_Abort() endif - if(token=='::') rc=-1 ! end-of-table is expected. No error message is produced end subroutine getnexttoken1_ subroutine getnexttoken2_(cf,token,myname) diff --git a/GEOSaana_GridComp/GSI_GridComp/analyzer b/GEOSaana_GridComp/GSI_GridComp/analyzer index 9d33f84..e5869d0 100755 --- a/GEOSaana_GridComp/GSI_GridComp/analyzer +++ b/GEOSaana_GridComp/GSI_GridComp/analyzer @@ -658,11 +658,25 @@ sub init { # NOTE: for a while, when using GSI, will need both because of spectral transforms # -------------------------------------------------------------------------------- if ( $doTRANSF ) { - if (-e "$fvInput/$myetc/newncepsfc.${jcap}" ) { - Assignfn( "$fvInput/$myetc/newncepsfc.${jcap}", "ncepsfc"); + if ( ! -e $rcname ) {die ">>>> ERROR <<< cannot find $rcname" }; + my $nveg = `nmlread.py $rcname GRIDOPTS nvege_type`; chomp($nveg); + print " this value: $nveg \n"; + if ( $nveg == 13 ) { + if (-e "$fvInput/$myetc/newncepsfc.${jcap}" ) { + Assignfn( "$fvInput/$myetc/newncepsfc.${jcap}", "ncepsfc"); + } else { + print "WARNING, analyzer: jcap inconsistent with available NCEP SFC file, reset to 254\n"; + Assignfn( "$fvInput/$myetc/newncepsfc.254", "ncepsfc"); + } } else { - print "WARNING, analyzer: jcap inconsistent with available NCEP SFC file, reset to 254\n"; - Assignfn( "$fvInput/$myetc/newncepsfc.254", "ncepsfc"); + my $mylat = `nmlread.py $rcname GRIDOPTS NLAT`; chomp($mylat); + my $mylon = `nmlread.py $rcname GRIDOPTS NLON`; chomp($mylon); + print "analyzer: NCEPSFC /discover/nobackup/projects/gmao/dadev/rtodling/TEST/GSIFIXSFC/JEDI/jedi.crtmsrf.${mylon}x${mylat}.nc4 \n"; + if (-e "/discover/nobackup/projects/gmao/dadev/rtodling/TEST/GSIFIXSFC/JEDI/jedi.crtmsrf.${mylon}x${mylat}.nc4" ) { + Assignfn( "/discover/nobackup/projects/gmao/dadev/rtodling/TEST/GSIFIXSFC/JEDI/jedi.crtmsrf.${mylon}x${mylat}.nc4", "ncepsfc"); + } else { + die ">>>> ERROR <<< cannot find jedi.crtmsrf.${mylon}x${mylat}.nc4"; + } } } @@ -1329,10 +1343,14 @@ sub ana { if ( $stg4hyb && ($stg4hyb ne "/dev/null") ) { print "satbias_out to $stg4hyb \n"; cp("satbias_out","$stg4hyb/$expid.ana.satbias.${nymd}_${hh}z.txt"); -# if ( $NEWRADBC && -e satbias_pc.out ) { -# print "satbias_pc.out to $stg4hyb \n"; -# cp ("satbias_pc.out", "$stg4hyb/$expid.ana.satbias_pc.${nymd}_${hh}z.txt"); -# } + if ( $NEWRADBC && -e satbias_pc.out ) { + print "satbias_pc.out to $stg4hyb \n"; + cp ("satbias_pc.out", "$stg4hyb/$expid.ana.satbias_pc.${nymd}_${hh}z.txt"); + } + if ( ($doABC > 0) && -e aircftbias_out ) { + print "aircftbias_out to $stg4hyb \n"; + cp("aircftbias_out","$stg4hyb/$expid.ana.acftbias.${nymd}_${hh}z.txt"); + } } } diff --git a/GEOSaana_GridComp/GSI_GridComp/m_nc_ncepsfc.f90 b/GEOSaana_GridComp/GSI_GridComp/m_nc_ncepsfc.f90 new file mode 100644 index 0000000..f5e82cf --- /dev/null +++ b/GEOSaana_GridComp/GSI_GridComp/m_nc_ncepsfc.f90 @@ -0,0 +1,349 @@ +module m_nc_ncepsfc +use netcdf +implicit none +private + +public :: nc_ncepsfc_vars_init +public :: nc_ncepsfc_vars_final +public :: nc_ncepsfc_vars +public :: nc_ncepsfc_dims +public :: nc_ncepsfc_read +public :: nc_ncepsfc_write +public :: nc_ncepsfc_getpointer + +type nc_ncepsfc_vars + logical :: initialized=.false. + integer :: nlon,nlat + integer :: nveg + real(4),pointer,dimension(:,:) :: vtype,stype,vfrac + real(4),pointer,dimension(:,:) :: v2d +end type nc_ncepsfc_vars + +character(len=*), parameter :: myname = 'm_nc_ncepsfc' + +integer, parameter :: nv2dx = 3 +character(len=5),parameter :: cvars2dx(nv2dx) = (/ 'stype', 'vtype', 'vfrac' /) + +interface nc_ncepsfc_dims; module procedure & + read_dims_ ; end interface +interface nc_ncepsfc_read; module procedure & + read_ncepsfc_ ; end interface +interface nc_ncepsfc_write; module procedure & + write_ncepsfc_ ; end interface +interface nc_ncepsfc_vars_init; module procedure & + init_ncepsfc_vars_ ; end interface +interface nc_ncepsfc_vars_final; module procedure & + final_ncepsfc_vars_ ; end interface +interface nc_ncepsfc_getpointer + module procedure get_pointer_2d_ +end interface + +contains + +subroutine read_dims_ (fname,nlat,nlon,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + integer, intent(out) :: rc + integer, intent(out) :: nlat,nlon + integer, intent(in), optional :: myid, root + +! This will be the netCDF ID for the file and data variable. + integer :: ncid, varid, ier + integer :: mype_,root_ + +! Local variables + character(len=*), parameter :: myname_ = myname//"::dims_" + +! Return code (status) + rc=0; mype_=0; root_=0 + if(present(myid) .and. present(root) ) then + mype_ = myid + root_ = root + endif + +! Open the file. NF90_NOWRITE tells netCDF we want read-only access to +! the file. + + call check_( nf90_open(fname, NF90_NOWRITE, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Read global attributes + call check_( nf90_inq_dimid(ncid, "lon", varid), rc, mype_, root_) + call check_( nf90_inquire_dimension(ncid, varid, len=nlon), rc, mype_, root_ ) + call check_( nf90_inq_dimid(ncid, "lat", varid), rc, mype_, root_ ) + call check_( nf90_inquire_dimension(ncid, varid, len=nlat), rc, mype_, root_ ) + +! Close the file, freeing all resources. + call check_( nf90_close(ncid), rc, mype_, root_ ) + + return + +end subroutine read_dims_ + +subroutine read_ncepsfc_ (fname,sfcvar,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + type(nc_ncepsfc_vars),intent(inout) :: sfcvar ! background error variables + integer, intent(out) :: rc + integer, intent(in), optional :: myid,root ! accommodate MPI calling programs + +! This will be the netCDF ID for the file and data variable. + integer :: ncid, varid + +! Local variables + character(len=*), parameter :: myname_ = myname//"::read_" + character(len=4) :: cindx + integer :: nv,nlat,nlon + integer :: ndims_, nvars_, ngatts_, unlimdimid_ + integer :: nlat_,nlon_ + integer :: mype_,root_ + real(4), allocatable :: data_in(:,:,:) + logical :: verbose + logical :: init_ + + +! Return code (status) + rc=0; mype_=0; root_=0 + verbose=.true. + init_=.false. + if(present(myid).and.present(root) )then + if(myid/=root) verbose=.false. + mype_ = myid + root_ = root + endif + +! Get dimensions + call read_dims_ (fname,nlat_,nlon_,rc, mype_,root_) + + init_ = sfcvar%initialized + if ( init_ ) then +! Set dims + nlat=sfcvar%nlat + nlon=sfcvar%nlon + +! Consistency check + if (nlon_ /= nlon .or. nlat_ /=nlat) then + rc=1 + if(myid==root) then + print *, 'nlat(file) = ', nlat_, 'nlat(required) = ', nlat + print *, 'nlon(file) = ', nlon_, 'nlon(required) = ', nlon + print *, myname_, 'Inconsistent dimensions, aborting ... ' + endif + return + endif + else +! Set dims + nlat=nlat_ + nlon=nlon_ + call init_ncepsfc_vars_(sfcvar,nlon,nlat) + endif + +! Open the file. NF90_NOWRITE tells netCDF we want read-only access to +! the file. + + call check_( nf90_open(fname, NF90_NOWRITE, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Read global attributes +! call check_( nf90_inquire(ncid, ndims_, nvars_, ngatts_, unlimdimid_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lon", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlon_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lat", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlat_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lev", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlev_), rc, mype_, root_ ) + +! Read in lat/lon fields + allocate(data_in(nlon,nlat,1)) + do nv = 1, nv2dx + call check_( nf90_inq_varid(ncid, trim(cvars2dx(nv)), varid), rc, mype_, root_ ) + call check_( nf90_get_var(ncid, varid, data_in(:,:,1)), rc, mype_, root_ ) + if(trim(cvars2dx(nv))=="vfrac" ) then + sfcvar%vfrac = transpose(data_in(:,:,1)) + endif + if(trim(cvars2dx(nv))=="vtype" ) then + sfcvar%vtype = transpose(data_in(:,:,1)) + endif + if(trim(cvars2dx(nv))=="stype" ) then + sfcvar%stype = transpose(data_in(:,:,1)) + endif + enddo + deallocate(data_in) + + sfcvar%nveg = nint(maxval(sfcvar%vtype))+1 + +! Close the file, freeing all resources. + call check_( nf90_close(ncid), rc, mype_, root_ ) + + if(verbose) print *,"*** Finish reading file: ", trim(fname) + + return + +end subroutine read_ncepsfc_ + +subroutine write_ncepsfc_ (fname,sfcvar,plevs,lats,lons,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + type(nc_ncepsfc_vars),intent(in) :: sfcvar ! background error variables + real(4), intent(in) :: lats(:) ! latitudes per GSI: increase index from South to North Pole + real(4), intent(in) :: lons(:) ! longitudea per GSI: increase index from East to West + real(4), intent(in) :: plevs(:) + integer, intent(out) :: rc + integer, intent(in), optional :: myid,root ! accommodate MPI calling programs + + character(len=*), parameter :: myname_ = myname//"::write_" + integer, parameter :: NDIMS = 3 + +! When we create netCDF files, variables and dimensions, we get back +! an ID for each one. + character(len=4) :: cindx + integer :: ncid, dimids(NDIMS) + integer :: x_dimid, y_dimid, z_dimid + integer :: lon_varid, lat_varid, lev_varid + integer :: ii,jj,nl,nv,nn,nlat,nlon + integer :: mype_,root_ + integer, allocatable :: varid1d(:), varid2d(:), varid2dx(:), varidMLL(:) + logical :: verbose + +! This is the data array we will write. It will just be filled with +! a progression of integers for this example. + real(4), allocatable :: data_out(:,:,:) + +! Return code (status) + rc=0; mype_=0; root_=0 + verbose=.true. + if(present(myid).and.present(root) )then + if(myid/=root) verbose=.false. + mype_ = myid + root_ = root + endif + +! Set dims + nlat=sfcvar%nlat + nlon=sfcvar%nlon + +! Always check the return code of every netCDF function call. In +! this example program, wrapping netCDF calls with "call check()" +! makes sure that any return which is not equal to nf90_noerr (0) +! will print a netCDF error message and exit. + +! Create the netCDF file. The nf90_clobber parameter tells netCDF to +! overwrite this file, if it already exists. + call check_( nf90_create(fname, NF90_CLOBBER, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Define the dimensions. NetCDF will hand back an ID for each. + call check_( nf90_def_dim(ncid, "lon", nlon, x_dimid), rc, mype_, root_ ) + call check_( nf90_def_dim(ncid, "lat", nlat, y_dimid), rc, mype_, root_ ) + + call check_( nf90_def_var(ncid, "lon", NF90_REAL, x_dimid, lon_varid), rc, mype_, root_ ) + call check_( nf90_def_var(ncid, "lat", NF90_REAL, y_dimid, lat_varid), rc, mype_, root_ ) + call check_( nf90_def_var(ncid, "lev", NF90_REAL, z_dimid, lev_varid), rc, mype_, root_ ) + + call check_( nf90_put_att(ncid, lon_varid, "units", "degress"), rc, mype_, root_ ) + call check_( nf90_put_att(ncid, lat_varid, "units", "degress"), rc, mype_, root_ ) + call check_( nf90_put_att(ncid, lev_varid, "units", "hPa"), rc, mype_, root_ ) + +! The dimids array is used to pass the IDs of the dimensions of +! the variables. Note that in fortran arrays are stored in +! column-major format. + dimids = (/ x_dimid, y_dimid, z_dimid /) + +! Define variables. + allocate(varid2dx(nv2dx)) + do nv = 1, nv2dx + call check_( nf90_def_var(ncid, trim(cvars2dx(nv)), NF90_REAL, (/ x_dimid, y_dimid /), varid2dx(nv)), rc, mype_, root_ ) + enddo + +! End define mode. This tells netCDF we are done defining metadata. + call check_( nf90_enddef(ncid), rc, mype_, root_ ) + +! Write coordinate variables data + call check_( nf90_put_var(ncid, lon_varid, lons ), rc, mype_, root_ ) + call check_( nf90_put_var(ncid, lat_varid, lats ), rc, mype_, root_ ) + +! Write out lat/lon fields + allocate(data_out(nlon,nlat,1)) + do nv = 1, nv2dx + if(trim(cvars2dx(nv))=="vfrac" ) then + data_out(:,:,1) = transpose(sfcvar%vfrac) + endif + if(trim(cvars2dx(nv))=="vtype" ) then + data_out(:,:,1) = transpose(sfcvar%vtype) + endif + if(trim(cvars2dx(nv))=="stype" ) then + data_out(:,:,1) = transpose(sfcvar%stype) + endif + call check_( nf90_put_var(ncid, varid2dx(nv), data_out(:,:,1)), rc, mype_, root_ ) + enddo + deallocate(data_out) + +! Close file + call check_( nf90_close(ncid), rc, mype_, root_ ) + + deallocate(varidMLL) + deallocate(varid2d) + deallocate(varid1d) + + if(verbose) print *,"*** Finish writing file: ", trim(fname) + + return + +end subroutine write_ncepsfc_ + +subroutine init_ncepsfc_vars_(vr,nlon,nlat) + + integer,intent(in) :: nlon,nlat + type(nc_ncepsfc_vars) vr + + if(vr%initialized) return + + vr%nlon=nlon + vr%nlat=nlat + +! allocate arrays + allocate(vr%vfrac(nlat,nlon),vr%vtype(nlat,nlon),vr%stype(nlat,nlon)) + vr%initialized=.true. + end subroutine init_ncepsfc_vars_ + + subroutine final_ncepsfc_vars_(vr) + type(nc_ncepsfc_vars) vr +! deallocate arrays + if(.not. vr%initialized) return + deallocate(vr%vtype,vr%vfrac,vr%stype) + vr%initialized=.false. +end subroutine final_ncepsfc_vars_ + +subroutine get_pointer_2d_ (vname, sfcvar, ptr, rc ) +implicit none +character(len=*), intent(in) :: vname +type(nc_ncepsfc_vars) sfcvar +real(4),pointer,intent(inout) :: ptr(:,:) +integer,intent(out) :: rc +rc=-1 +if(trim(vname)=='vtype') then + ptr => sfcvar%vtype + rc=0 +endif +if(trim(vname)=='vfrac') then + ptr => sfcvar%vfrac + rc=0 +endif +if(trim(vname)=='stype') then + ptr => sfcvar%stype + rc=0 +endif +end subroutine get_pointer_2d_ + +subroutine check_(status,rc, myid, root) + integer, intent ( in) :: status + integer, intent (out) :: rc + integer, intent ( in) :: myid, root + rc=0 + if(status /= nf90_noerr) then + if(myid==root) print *, trim(nf90_strerror(status)) + rc=999 + end if +end subroutine check_ + +end module m_nc_ncepsfc From d333c06507cf5ccd0f3af170de8a34b1ac09a60c Mon Sep 17 00:00:00 2001 From: Ricardo Todling Date: Thu, 24 Oct 2024 16:15:57 -0400 Subject: [PATCH 2/2] info on change --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index cdd9eef..d667087 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - add fix to read_bufrtovs to handle ta2tb=.true. when there are multiple versions of SpcCoeff.bin file for single instrument/platform +- add ability to handle nvege_type=20 - no longer doing + fishy interpolation. ### Fixed