diff --git a/CHANGELOG.md b/CHANGELOG.md index 87275fd9..aec5efaf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added statsNx.rc for screen level variable fstats +### Changed + +- Added options for land-only and screen-level variables in fstats.x and g5fcst_stats.pl + +## [Unreleased] + +### Added + ### Changed - Added a few new tags to `regrid.pl` diff --git a/GEOS_Util/post/CMakeLists.txt b/GEOS_Util/post/CMakeLists.txt index 49656e96..3d5d901e 100644 --- a/GEOS_Util/post/CMakeLists.txt +++ b/GEOS_Util/post/CMakeLists.txt @@ -141,6 +141,7 @@ set (post_rc plot.rc post.rc stats.rc + statsNx.rc time_ave.rc 3CH.rc ) @@ -151,7 +152,7 @@ install ( ) install ( - FILES stats.rc time_ave.rc + FILES stats.rc time_ave.rc statsNx.rc DESTINATION etc ) diff --git a/GEOS_Util/post/g5fcst_stats.pl b/GEOS_Util/post/g5fcst_stats.pl index 436601c4..e764b141 100755 --- a/GEOS_Util/post/g5fcst_stats.pl +++ b/GEOS_Util/post/g5fcst_stats.pl @@ -1,6 +1,6 @@ #!/usr/bin/env perl #======================================================================= -# name - g5fcst_stats.pl +# name - g5fcst_stats.pl # purpose - script to submit jobs to calculate forecast statistics # # key global variables - @@ -29,7 +29,11 @@ my ($offset_sec, $pesto, $progdir, $progtype, $secs_per_day, $secs_per_hour); my ($statsX, $storedir, $tau_freq, $tau_fsec, $vanadir, $vanatype, $vexpid); my (%opts, @acqIDs, @statsIDs); +my ($landonly, $landonly_dflt); +my ($nxonly, $nxonly_dflt); +$landonly_dflt = "no"; +$nxonly_dflt = "no" ; $fhours_dflt = 123; $localID = $$; $ncsuffix = ".nc4"; @@ -135,7 +139,11 @@ if ($ihh eq "00" and $nv == 1) { $fdir = "$vanadir/Y$nyyyy/M$nmm"; - $fname = "$vexpid.$vanatype.inst3d_met_p.${vdate}_${vhh}${vmn}z"; + if ($nxonly eq "yes") { + $fname = "$vexpid.$vanatype.${vdate}_${vhh}${vmn}z"; + } else { + $fname = "$vexpid.$vanatype.inst3d_met_p.${vdate}_${vhh}${vmn}z"; + } $ffile0 = "$fdir/$fname$ncsuffix"; } else { @@ -167,7 +175,11 @@ push @alist, "$vexpid.$vanatype.${vdate}_${vhh}z+${vdate}_${vhh}${vmn}z"; } else { - push @alist, "$vexpid.$vanatype.inst3d_met_p.${vdate}_${vhh}${vmn}z"; + if ($nxonly eq "yes" ) { + push @alist, "$vexpid.$vanatype.${vdate}_${vhh}${vmn}z"; + } else { + push @alist, "$vexpid.$vanatype.inst3d_met_p.${vdate}_${vhh}${vmn}z"; + } } # loop through naming convention options @@ -192,9 +204,12 @@ #---------------------------------------------------- ($vdate, $vtime) = tick($vdate, $vtime, $tau_fsec); } - + if ($nxonly eq "yes") { + @climfiles = (<$ENV{SHARE}/dao_ops/verification/stats/MERRA-2.inst1_2d_asm_Nx.198501_201412.clim_??z.576x361.data.nc4>); + } else { @climfiles = (<$ENV{SHARE}/dao_ops/verification/stats/MERRA-2.inst3_3d_asm_Np.198501_201412.clim_??z.576x361.data.nc4>); #--@climfiles = (<$ENV{SHARE}/dao_ops/verification/stats/merrasc.197901-200812.clim_??z.144x91.data.nc>); + } $climfilecnt = scalar(@climfiles); if ($climfilecnt < 4) { rmtree($fstatswork) if -d $fstatswork; @@ -221,15 +236,23 @@ } else { push @fcst_fnames, $ffile0 if $ffile0; - push @fcst_fnames, "$expid.prog.$progtype.${ndate}_${ihh}z+*$ncsuffix"; - push @ana_fnames, "$vexpid.asm.inst3d_met_p.*$ncsuffix"; + push @fcst_fnames, "$expid.prog.$progtype.${ndate}_${ihh}z+*$ncsuffix"; + if ($nxonly eq "yes" ) { + push @ana_fnames, "$vexpid.$vanatype.*$ncsuffix"; + } else { + push @ana_fnames, "$vexpid.asm.inst3d_met_p.*$ncsuffix"; + } } } else { push @fcst_fnames, $ffile0 if $ffile0; push @fcst_fnames, "$expid.prog.$progtype.${ndate}_${ihh}z+*$ncsuffix"; - push @ana_fnames, "$vexpid.asm.inst3d_met_p.*$ncsuffix"; - push @ana_fnames, "$vexpid.ana.inst3d_met_p.*$ncsuffix"; + if ($nxonly eq "yes" ) { + push @ana_fnames, "$vexpid.$vanatype.*$ncsuffix" + } else { + push @ana_fnames, "$vexpid.asm.inst3d_met_p.*$ncsuffix"; + push @ana_fnames, "$vexpid.ana.inst3d_met_p.*$ncsuffix"; + } } @fcstlist = sort keys %fcsthash; @vanalist = sort keys %vanahash; @@ -282,6 +305,8 @@ $args{"vhh0"} = $vhh0; $args{"vhh1"} = $vhh1; $args{"nodes"} = $nodes; + $args{"landonly"} = $landonly; + $args{"nxonly"} = $nxonly; # submit job to calculate stats #------------------------------ @@ -344,6 +369,10 @@ sub init { "nodes=s" => \$nodes, + "landonly=s" => \$landonly, + + "nxonly=s" => \$nxonly, + "das" => \$dasFLG, "np" => \$noprompt, @@ -359,7 +388,15 @@ sub init { $expid = shift @ARGV; $idate = shift @ARGV; $ndays = shift @ARGV; - + + # landonly option + #---------------- + $landonly = $landonly_dflt unless $landonly; + + #nxonly option + #-------------- + $nxonly = $nxonly_dflt unless $nxonly; + # initial fcst hour and offset #----------------------------- unless ($ihh) { @@ -536,9 +573,12 @@ sub init { # find stats.rc #-------------- $fv_etcdir = dirname($Bin) ."/etc"; - die "Error. Cannot find $fv_etcdir/stats.rc" unless -e "$fv_etcdir/stats.rc"; + if ($nxonly eq "yes" ) { + die "Error. Cannot find $fv_etcdir/statsNx.rc" unless -e "$fv_etcdir/statsNx.rc"; + } else { + die "Error. Cannot find $fv_etcdir/stats.rc" unless -e "$fv_etcdir/stats.rc"; + } } - #======================================================================= # name - das_check # purpose - check for hidden files in $FVHOME/fcst and $FVHOME directories @@ -636,11 +676,14 @@ sub submit_calcjob { my (@rmfilelist, $yyyy, $mm, $dd); my ($logdir, $logfile1, $logfile2, $jobname, $jobdate, $jobfile, $jobtype); my ($cmd, $jobID, $jobIDline); - my (@levs, @levels_19, @levels_11); + my (@levs, @levels_19, @levels_11, @levels_1); my ($mynodes,$usrnodes); my ($qos, $partition); my ($ntspn, $npn); - + my ($landonly, $landmaskdirfile); + my ($nxonly); + my ($whichrc); + @levels_19 = ( 1000.0, 975.0, 950.0, 925.0, 900.0, 850.0, 800.0, 750.0, 700.0, 600.0, 500.0, 400.0, @@ -651,6 +694,8 @@ sub submit_calcjob { 500.0, 400.0, 300.0, 250.0, 200.0, 150.0, 100.0 ); + @levels_1 = ( 1000.0 ) ; + if ($vexpid eq "ecmwf") { @levs = @levels_11 } else { @levs = @levels_19 } @@ -668,6 +713,22 @@ sub submit_calcjob { $vdate1 = $args{"vdate1"}; $vhh1 = $args{"vhh1"}; $usrnodes = $args{"nodes"}; + $landonly = $args{"landonly"}; + $nxonly = $args{"nxonly"}; + print "nxonly: $nxonly\n"; + print "landonly: $landonly\n"; + + if ($landonly eq "yes") { + $landmaskdirfile = "$jobdir/landmaskfile.nc4"; + die "Error. Cannot find $landmaskdirfile" unless -e "$landmaskdirfile"; + } + + $whichrc = "stats.rc"; + if ($nxonly eq "yes") { + @levs = @levels_1; + $whichrc = "statsNx.rc"; + print "whichrc: $whichrc\n"; + } foreach (@fcstlist, @vanalist) { push @rmfilelist, basename($_) }; @@ -736,12 +797,24 @@ sub submit_calcjob { -cli @climfiles \\ -tag $expid.${ihh}z \\ -nfreq ${tau_freq}0000 \\ +EOF +if ( $landonly eq "yes" ) { +print FH <<"EOF"; + -land ${landmaskdirfile} \\ +EOF +} + +print FH <<"EOF"; -levs @levs \\ -o $expid.fstats.log.$jobdate.txt \\ -verif gmao \\ -fcsrc gmao \\ -fhour $fhours \\ - -rc $fv_etcdir/stats.rc + -rc $fv_etcdir/$whichrc +EOF + +print FH <<"EOF"; + @ calc_status = \$status $pesto -arc $arcfile \\ @@ -940,7 +1013,8 @@ sub verify_values { print "progdir: $progdir\n"; print "progtype: $progtype\n"; print "fs_tag: $fs_tag\n\n"; - + print "landonly: $landonly\n\n"; + print "nxonly: $nxonly\n\n"; print "dryrun: $dryrun\n\n" if $dryrun; $ans = query("Continue (y/n):", "y"); @@ -986,6 +1060,12 @@ sub usage { [dirname(\$FVHOME) or \$NOBACKUP] -noarchive do not archive outputs [archives by default] + -landonly calculate stats only over land surface (yes/no) + [no by default] + + -nxonly calculate stats of 2d Nx collection only (yes/no) + [no by default] + -nodes nodesname specify nodes (e.g., sky, hasw, or cas) -das check for DAS hidden files before attempting to fetch files and set no prompt; requires \$FVHOME environment variable; diff --git a/GEOS_Util/post/stats.F90 b/GEOS_Util/post/stats.F90 index 6b616b5e..03dfd2f8 100644 --- a/GEOS_Util/post/stats.F90 +++ b/GEOS_Util/post/stats.F90 @@ -56,6 +56,13 @@ program stats real, allocatable :: rms(:,:,:,:,:) ! Note: Hardwired for 100 time periods (Max) real*4 dum(nr) +!for land-only option +! ------------------- + real, allocatable :: landmask(:,:) + logical landonly + data landonly /.false./ + character*256 lndfname + ! Original Levels ! --------------- ! real zlev0(10) @@ -277,6 +284,12 @@ end subroutine init_levs if( trim(arg(n)).eq.'-o' ) read(arg(n+1),fmt='(a)') outfile if( trim(arg(n)).eq.'-verif' ) read(arg(n+1),*) averify if( trim(arg(n)).eq.'-fcsrc' ) read(arg(n+1),*) fcsource +!for land-only option +!-------------------- + if( trim(arg(n)).eq.'-land' ) then + read(arg(n+1),fmt='(a)') lndfname + landonly = .true. + endif enddo endif @@ -335,7 +348,7 @@ end subroutine init_levs if(failed) then print* ,' ERROR: when GMAOpy output requested the following must be specified:' print* ,' -fcsrc FORECAST (e.g., gmao)' - print* ,' -verif VERIFICATION (e.g., ncep)' + print* ,' -verif VERIFICATION (e.g., ecmwf or ncep)' print* ,' Aborting ...' call exit(1) endif @@ -702,6 +715,14 @@ end subroutine init_levs dl = 2*pi/im dp = pi/(jm-1) +!for land-only option +!-------------------- + if ( landonly ) then + allocate ( landmask (im,jm ) ) + print *, ' landmask filename: ', trim(lndfname) + call load_landmask ( lndfname,landmask,im,jm,undef ) + endif + ! Loop over Forecast Times ! ------------------------ nt = 0 @@ -792,6 +813,16 @@ end subroutine init_levs call writit ( fields_3d(n)%clim,im,jm,nl,cundef,undef ) enddo +!for land-only option +!-------------------- + if (landonly) then + do n=1,n2d + call landmk ( fields_2d(n)%fcst,im,jm,1 ,landmask,undef ) + enddo + do n=1,n3d + call landmk ( fields_3d(n)%fcst,im,jm,nl,landmask,undef ) + enddo + endif ! Loop over Geographical Regions ! ------------------------------ @@ -2360,6 +2391,7 @@ subroutine bin_10x10 ( z,im,jm,z10x10 ) subroutine read_anal( nymd,nhms,fields_2d,fields_3d,n2d,n3d,idim,jdim,nl,zlev,ana_files,num_ana_files,undef ) use stats_mod + use MAPL_ConstantsMod implicit none type(fields) :: fields_2d(n2d) type(fields) :: fields_3d(n3d) @@ -2383,6 +2415,7 @@ subroutine read_anal( nymd,nhms,fields_2d,fields_3d,n2d,n3d,idim,jdim,nl,zlev,an character*256, allocatable :: vunits(:) real, allocatable :: q(:,:) + real, allocatable :: sp(:,:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) @@ -2400,6 +2433,8 @@ subroutine read_anal( nymd,nhms,fields_2d,fields_3d,n2d,n3d,idim,jdim,nl,zlev,an logical check_names logical shift, defined, first integer num, LL, i,j,k, loc, ndates, len + real Td, pp, ee + real tice, epsln, ec0, ec1, ec2 data id /0/ data num /0/ data shift /.false./ @@ -2522,6 +2557,34 @@ subroutine read_anal( nymd,nhms,fields_2d,fields_3d,n2d,n3d,idim,jdim,nl,zlev,an if( check_names( vname(n),fields_2d(m)%alias(k) ) ) then found_2d(m) = .true. call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,0,1,q,rc ) + + if( trim(vname(n)) == 'N2_metre_dewpoint_temperature' ) then + tice= MAPL_TICE + epsln = MAPL_EPSILON + !! empirical coeff. (Bolton 1980) + ec0 = 6.112 + ec1 = 17.67 + ec2 = 243.5 + + allocate( sp(im, jm) ) + call gfio_getvar ( id,'Surface_pressure',nymd,nhms,im,jm,0,1,sp,rc ) + do j=1,jm + do i=1,im + if( defined(q(i,j),undef) .and. defined(sp(i,j),undef) )then + Td = q(i,j)-tice ! to C + pp = sp(i,j)/100. ! to mb + ee= ec0*exp((ec1*Td)/(Td + ec2)) + q(i,j) = (epsln * ee)/(pp - (1.0-epsln) * ee) + else + print*, 'Td conversion fails,set q2m to UNDEF' + q(i,j) = undef + endif + enddo + enddo + deallocate ( sp ) + endif + !------ + if( shift ) call hshift ( q,im,jm ) if( check_names( fields_2d(m)%name,'p' ) ) then do j=1,jm @@ -3181,6 +3244,106 @@ subroutine writit ( q,im,jm,lm,qundef,undef ) return end +!for land-only option +!-------------------- + subroutine landmk ( q,im,jm,lm,landmask,undef ) + implicit none + integer im,jm,lm + real q(im,jm,lm) + real*4 undef + real*4 landmask(im,jm) + integer i,j,L + + do L=1,lm + do j=1,jm + do i=1,im + if(landmask(i,j) == 0 ) then + q(i,j,L) = undef + endif + enddo + enddo + enddo + return + end + +! for land-only option +!-------------------- + subroutine load_landmask (filename,frland,imbin,jmbin,undef ) + use stats_mod + implicit none + + character*256 filename + integer id,im,jm,lm,nvars,rc + integer ntime,ngatts,timinc + character*256 title + character*256 source + character*256 contact + character*256 levunits + character*256, allocatable :: vname(:) + character*256, allocatable :: vtitle(:) + character*256, allocatable :: vunits(:) + + real, allocatable :: q(:,:) + real, allocatable :: lat(:) + real, allocatable :: lon(:) + real, allocatable :: lev(:) + real, allocatable :: vrange(:,:) + real, allocatable :: prange(:,:) + integer, allocatable :: yymmdd(:) + integer, allocatable :: hhmmss(:) + integer, allocatable :: kmvar(:) + + integer imbin,jmbin, n + character*256 landname + real frland (imbin,jmbin ) + logical shift + real undef + + call gfio_open ( trim(filename) ,1,id,rc ) + call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) + print *, ' landmask dimension: ',im, jm + allocate ( lon(im) ) + allocate ( lat(jm) ) + allocate ( lev(lm) ) + allocate ( yymmdd(ntime) ) + allocate ( hhmmss(ntime) ) + allocate ( vname(nvars) ) + allocate ( vtitle(nvars) ) + allocate ( vunits(nvars) ) + allocate ( kmvar(nvars) ) + allocate ( vrange(2,nvars) ) + allocate ( prange(2,nvars) ) + call gfio_inquire ( id,im,jm,lm,ntime,nvars, & + title,source,contact,undef, & + lon,lat,lev,levunits, & + yymmdd,hhmmss,timinc, & + vname,vtitle,vunits,kmvar, & + vrange,prange,rc ) + shift = .false. + if( lon(1).eq.0.0 ) then + print *, 'landmask begins at lon: ',lon(1) + print *, 'Horizontal Shift will be performed' + shift = .true. + endif + + do n=1, nvars + if ( trim(vname(n)) == 'FRLAND') then + landname = trim(vname(n)) + print *, ' landmask varname :', trim(landname) + endif + enddo + + allocate (q(im,jm)) + call gfio_getvar ( id,trim(landname),yymmdd,hhmmss, & + im,jm,0,1,q,rc ) + + if( shift ) call hshift ( q,im,jm ) + call bin ( q,im,jm,frland,imbin,jmbin,undef,0 ) + call gfio_close(id,rc) + return + end + + subroutine usage() print *, "Usage: " print * @@ -3190,6 +3353,7 @@ subroutine usage() print *, " <-tag tag>" print *, " <-pref PREF>" print *, " <-nfreq HHMMSS>" + print *, " <-landonly landmaskfile>" print *, " <-syserror SYSERR>" print * print *, "where:" @@ -3205,6 +3369,7 @@ subroutine usage() print *, " -pref PREF : Optional Reference Pressure for STD Output Printing" print *, " (Default: 500-mb)" print *, " -nfreq HHMMSS : Optional Frequency for Forecast Time Periods" + print *, " -landonly landmaskfile : Optional calculate stats overland" print *, " -syserr SYSERR File : Optional Systematic Error File" print *, " to be Subtracted from Forecasts" print * diff --git a/GEOS_Util/post/statsNx.rc b/GEOS_Util/post/statsNx.rc new file mode 100644 index 00000000..46bfc4af --- /dev/null +++ b/GEOS_Util/post/statsNx.rc @@ -0,0 +1,86 @@ +# ------------------- +# Stats RC Parameters +# ------------------- +# +# COLLECTIONS: List of Output Collections +# NOTE: A collection name of: default +# will produce stats files with names consistent with older versions. +# +# fields_2d: List of 2D Fields for each Collection +# fields_3d: List of 3D Fields for each Collection +# +# TYPE_x: Type Attribute for Field x (can be from ANY collection) +# Options: scalar ---- (implies wave-0 pole value) +# vector ---- (implies wave-1 pole value) +# aerosol --- (implies wave-0 pole value and math: log(x+0.01)) +# +# ALIAS_x: All Possible Aliases for Field x +# +# --------------------------------------------------------------------------------------------------- + +## COLLECTIONS: default + COLLECTIONS: inst1_2d_asm_Nx + +# Collection Descriptions: +# ------------------------ + + default.fields_2d: p + default.fields_3d: u v t q h + + inst1_2d_asm_Nx.fields_2d: u10m v10m t2m q2m + + aerosol.fields_2d: tau du ss su bc oc + +# --------------------------------------------------------------------------------------------------- + + TYPE_p: scalar + ALIAS_p: slp slprs Surface_pressure Mean_sea_level_pressure msl + + TYPE_u10m: vector + ALIAS_u10m: u10m N10_metre_U_wind + + TYPE_v10m: vector + ALIAS_v10m: v10m N10_metre_V_wind + + TYPE_t2m: scalar + ALIAS_t2m: t2m N2_metre_temperature + + TYPE_q2m: scalar + ALIAS_q2m: qv2m N2_metre_dewpoint_temperature + +# --------------------------------------------------------------------------------------------------- + + TYPE_u: vector + ALIAS_u: u uwnd ugrd ugrdprs u_velocity + + TYPE_v: vector + ALIAS_v: v vwnd vgrd vgrdprs v_velocity + + TYPE_t: scalar + ALIAS_t: t tmpu tmp tmpprs temperature + + TYPE_q: scalar + ALIAS_q: q sphu qv + + TYPE_h: scalar + ALIAS_h: h hght hgt hgtprs height + +# --------------------------------------------------------------------------------------------------- + + TYPE_tau: aerosol + ALIAS_tau: totexttau + + TYPE_du: aerosol + ALIAS_du: duexttau + + TYPE_ss: aerosol + ALIAS_ss: ssexttau + + TYPE_su: aerosol + ALIAS_su: suexttau + + TYPE_bc: aerosol + ALIAS_bc: bcexttau + + TYPE_oc: aerosol + ALIAS_oc: ocexttau