Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow for consistency of NCEP-sfc fields with those used in JEDI #180

Merged
merged 2 commits into from
Nov 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
1 change: 1 addition & 0 deletions GEOSaana_GridComp/GSI_GridComp/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
139 changes: 98 additions & 41 deletions GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand All @@ -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
!-------------------------------------------------------------------------
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -2680,22 +2711,29 @@ 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
else
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
Expand All @@ -2707,45 +2745,58 @@ 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
else
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
else
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
Expand All @@ -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)
Expand Down Expand Up @@ -2810,7 +2864,7 @@ subroutine GSI_GridCompGetNCEPsfcFromFile_(lit)
#endif

end subroutine GSI_GridCompGetNCEPsfcFromFile_

!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
Expand Down Expand Up @@ -3421,18 +3475,21 @@ 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
write(6,'(2a,i5)') myname, ': ESMF_ConfigGetAttribute(field #1) error, rc=', 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)
Expand Down
34 changes: 26 additions & 8 deletions GEOSaana_GridComp/GSI_GridComp/analyzer
Original file line number Diff line number Diff line change
Expand Up @@ -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";
}
}
}

Expand Down Expand Up @@ -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");
}
}
}

Expand Down
Loading