gblevents.f

Go to the documentation of this file.

C> @file
C> @brief Pre/Post processing of of prepbufr events
      MODULE gblevn_module
 
      IMPLICIT NONE
 
      SAVE
 
      INTEGER                    :: imax
      INTEGER                    :: jmax
      INTEGER                    :: kmax
      INTEGER                    :: kmaxs
      INTEGER*4                  :: idvc
      INTEGER*4                  :: idsl
      INTEGER*4                  :: nvcoord
      INTEGER*4                  :: sfcpress_id
      INTEGER*4                  :: thermodyn_id
      REAL (kind=8), allocatable :: iar14t(:,:,:) 
      REAL (kind=8), allocatable :: iar15u(:,:,:) 
      REAL (kind=8), allocatable :: iar16v(:,:,:) 
      REAL (kind=8), allocatable :: iar17q(:,:,:) 
      REAL (kind=8), allocatable :: iar12z(:,:) 
      REAL (kind=8), allocatable :: iar13p(:,:) 
      REAL (kind=8), allocatable :: iarpsi(:,:,:) 
      REAL (kind=8), allocatable :: iarpsl(:,:,:) 
      REAL (kind=8), allocatable :: iarpsd(:,:,:) 
      REAL (kind=4), allocatable :: vcoord(:,:)  
      REAL                       :: dlat
      REAL                       :: dlon
 
      END MODULE gblevn_module
 
C> SUBPROGRAM:    GBLEVENTS   PRE/POST PROCESSING OF PREPBUFR EVENTS 
C>   PRGMMR: DENNIS KEYSER    ORG: EMC        DATE: 2017-02-22
C>
C> RUNS IN TWO MODES: "PREVENTS" AND "POSTEVENTS".  IN THE
C> PREVENTS MODE, PREPARES OBSERVATIONAL PREPBUFR REPORTS FOR
C> SUBSEQUENT QUALITY CONTROL AND ANALYSIS PROGRAMS.  THIS IS DONE
C> THROUGH THE FOLLOWING: INTERPOLATION OF GLOBAL FIRST GUESS {FROM
C> EITHER SIGIO (SIGMA OR HYBRID) OR NEMSIO INPUT} TO PREPBUFR
C> OBSERVATION LOCATIONS WITH ENCODING OF FIRST GUESS VALUES INTO
C> PREPBUFR REPORTS; ENCODING OF "PREVENT" AND/OR "VIRTMP" EVENTS INTO
C> PREPBUFR REPORTS; IN CERTAIN CASES, ENCODING OF A DERIVED PMSL INTO
C> SURFACE PREPBUFR REPORTS; AND ENCODING OF OBSERVATION ERRORS FROM
C> THE ERROR SPECIFICATION FILE INTO PREPBUFR REPORTS.  IN THE
C> POSTEVENTS MODE, AFTER ALL QUALITY CONTROL AND ANALYSIS PROGRAMS
C> HAVE RUN, INTERPOLATES THE GLOBAL ANALYSIS {FROM EITHER SIGIO
C> (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT} TO PREPBUFR OBSERVATION
C> LOCATIONS AND ENCODES THESE ANALYZED VALUES INTO PREPBUFR REPORTS.
C> THE REMAINDER OF THIS ABSTRACT APPLIES ONLY TO THE PREVENTS MODE.
C> THE "PREVENT" EVENT CAN CHANGE A QUALITY MARKER TO FLAG AN
C> OBSERVATION DATUM FOR NON-USE BY SUBSEQUENT QC AND ANALYSIS
C> PROGRAMS (FILTERING).  EXAMPLES WHERE THIS SUBROUTINE WILL WRITE
C> AN EVENT TO FLAG A DATUM INCLUDE: THE OBSERVATION ERROR FOR THAT
C> DATUM IS READ IN AS MISSING IN THE INPUT ERROR FILE, THE DATUM
C> ITSELF VIOLATES A GROSS OR "SANITY" CHECK, OR THE OBSERVED
C> PRESSURE DATUM IS MORE THAN 100 MB BELOW THE GUESS SURFACE
C> PRESSURE.  THE "VIRTMP" EVENT CAN CHANGE THE SPECIFIC HUMIDITY
C> OBSERVATION (RE-CALCULATED) AS WELL AS THE TEMPERATURE
C> OBSERVATION (FROM SENSIBLE TO VIRTUAL TEMPERATURE, BASED ON
C> JUST-CALCULATED SPECIFIC HUMIDITY).  CURRENTLY THIS APPLIES ONLY
C> TO SURFACE (LAND, MARINE AND MESONET) DATA TYPES, POSSIBLY TO
C> RAOB, DROP AND MULTI-LEVEL RECCO DATA TYPES IF THE SWITCH
C> "ADPUPA_VIRT" IS TRUE (NORMALLY, HOWEVER IT IS FALSE) [OTHER DATA
C> TYPES WITH REPORTED SENSIBLE TEMPERATURE EITHER HAVE MISSING
C> MOISTURE (E.G., ALL AIRCRAFT TYPES EXCEPT FOR SOME ACARS, SATELLITE
C> WIND TYPES), FLAGGED MOISTURE (E.G., SOME ACARS) OR CALCULATE
C> SPECIFIC HUMIDITY/VIRTUAL TEMPERATURE IN SUBSEQUENT PROGRAMS (E.G.,
C> RAOBS, DROPS AND MULTI-LEVEL RECCOS WHICH CALCULATE THESE IN
C> PROGRAM "CQCBUFR", IN WHICH CASE THE SWITCH "ADPUPA_VIRT" HERE MUST
C> BE FALSE!)]. FOR CASES WHERE THE SWITCH "DOBERR" IS FALSE, THE
C> OBSERVATION ERROR FOR ALL DATA REMAINS MISSING IN THE PREPBUFR
C> FILE.  IN THIS CASE, THE INPUT ERROR FILE IS USUALLY A NULL FILE
C> AND THE "PREVENT" EVENT TO FLAG THE DATUM IS NOT INVOKED.  FOR
C> CASES WHERE THE SWITCH "DOFCST" IS FALSE, IF THE SWITCH "SOME_FCST"
C> IS ALSO FALSE, THEN FORECAST VALUES ARE NOT ENCODED FOR ANY MESSAGE
C> TYPE; IF "SOME_FCST" IS TRUE THEN FORECAST VALUES ARE ENCODED, BUT
C> ONLY FOR REPORTS IN THOSE MESSAGE TYPES FOR WHICH A GUESS VALUE IS
C> NEEDED BY SUBSEQUENT QC PROGRAMS.  IT SHOULD BE NOTED THAT THE
C> FILTERING OF DATA ASSOCIATED WITH THE "PREVENT" EVENT PROCESSING IS
C> NOT INVOKED IF ALL THREE ARE TRUE: DOBERR= FALSE, THE FORECAST
C> VALUES ARE MISSING (DOFCST=FALSE & SOME_FCST=TRUE & MESSAGE TYPE IS
C> NOT "ADPUPA", "AIRCFT", "AIRCAR", "PROFLR", OR "VADWND" -- OR  --
C> DOFCST=FALSE & SOME_FCST=FALSE), AND "VIRTMP" EVENT PROCESSING IS
C> NOT INVOKED (EITHER MESSAGE TYPE IS NOT "ADPSFC", "SFCSHP" OR
C> "MSONET" WHEN "ADPUPA_VIRT" IS FALSE, OR MESSAGE TYPE IS NOT
C> "ADPSFC", "SFCSHP", "MSONET" OR "ADPUPA" WHEN "ADPUPA_VIRT" IS
C> TRUE). ALSO, IF VIRTUAL TEMPERATURE PROCESSING IS PERFORMED, ALL
C> SURFACE REPORTS WITH MISSING PMSL WILL ENCODE A DERIVED PMSL INTO
C> PREPBUFR IF THE SWITCH DOPMSL IS TRUE AND A VIRTUAL TEMPERATURE WAS
C> SUCCESSFULLY CALCULATED.
C>
C> PROGRAM HISTORY LOG:
C> -1999-07-01 D. A. KEYSER -- ORIGINAL AUTHOR (ADAPTED FROM PREVENTS
C>     SUBROUTINE IN PREPDATA PROGRAM, BUT NOW GENERALIZED FOR
C>     POSTEVENTS MODE)
C> -1999-07-12 D. A. KEYSER -- MODIFIED TO INTERPOLATE MODEL SPECIFIC
C>     HUMIDITY TO OBSERVATION LOCATION WHEN OBS. SPECIFIC HUMIDITY IS
C>     MISSING AS LONG AS OBS. TEMPERATURE IS NON-MISSING
C> -1999-09-09 D. A. KEYSER -- ADDED "VADWND" TO THE LIST OF MESSAGE
C>     TYPES FOR WHICH FORECAST VALUES MUST BE ENCODED, EVEN WHEN
C>     DOFCST=FALSE (NECESSARY BECAUSE THE NEW PROGRAM CQCVAD NEEDS THE
C>     BACKGROUND DATA)
C> -1999-09-09 D. A. KEYSER -- CHANGES TO MAKE CODE MORE PORTABLE; 
C>     'TFC' NOW GENERATED FOR VADWND MESSAGE TYPES EVEN THOUGH TOB IS
C>     MISSING (NEEDED BY CQCVAD PROGRAM)
C> -1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE
C>     NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT
C>     TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD
C>     QUAL. MARKER (8)
C> -2000-09-21 D. A. KEYSER -- THE PRESSURE LEVEL ABOVE WHICH ALL SPEC.
C>     HUMIDITY QUAL. MARKERS ARE "REJECTED" (Q.M. SET TO 9) IS NOW READ
C>     IN AS A N-LIST SWITCH (QTOP_REJ), BEFORE IT WAS HARDWIRED TO 300
C>     MB
C> -2000-12-13 D. A. KEYSER -- WILL NO LONGER PERFORM VIRTUAL TEMPERATURE
C>     PROCESSING FOR ACARS DATA SINCE MOISTURE IS FLAGGED RIGHT NOW
C>     (ACARS MOISTURE ONLY WRITTEN INTO PREPBUFR FILE FOR STATISTICAL
C>     REASONS)
C> -2001-02-02 D. A. KEYSER -- RESTORED LEGACY LOGIC TO FLAG CERTAIN
C>     SATELLITE TEMPERATURE SOUNDINGS EITHER BELOW 100 MB (TEMP. OBS)
C>     OR ON ALL LEVELS (SPEC. HUM. OBS), CONTROLLED BY NEW NAMELIST
C>     SWITCH "SATMQC"
C> -2001-09-27 D. A. KEYSER -- 'TFC' AND 'QFC' NOW GENERATED FOR REPORT
C>     TYPE 111 (SYNDAT REPORTS AT STORM CENTER) EVEN THOUGH "TOB" AND
C>     "QOB" ARE MISSING (NEEDED BY SYNDATA PROGRAM); IN PREPARATION FOR
C>     CHANGE FROM T170L42 TO T254L64 SGES, NOW MAKES COEFFICIENT ARRAYS
C>     ALLOCATABLE TO ALLOW THEM TO OBTAIN MEMORY FROM "HEAP" RATHER
C>     THAN FROM "STACK", ALSO HAVE INCREASED THE MAX NUMBER OF LEVELS
C>     IN ARRAYS FROM 42 TO 64, FINALLY ALSO NO LONGER STOPS WITH C.
C>     CODE 70 IF EVEN NUMBER OF LONGITUDES IN SIGMA GUESS (IMAX,
C>     HARDWIRED TO 384) IS .LT. SPECTRAL RESOLUTION (JCAP) * 2
C> -2001-10-10 D. A. KEYSER -- AT PREPBUFR CENTER DATES WITH AN HOUR THAT
C>     IS NOT A MULTIPLE OF 3 (WHEN A GLOBAL SIGMA GUESS/ANAL FILE IS
C>     NOT AVAILABLE; E.G., IN RUC2A RUNS) NOW PERFORMS A LINEAR
C>     INTERPOLATION BETWEEN SPECTRAL COEFFICIENTS IN 2 SPANNING SIGMA
C>     GUESS/ANAL FILES 3-HRS APART TO CENERATE A GUESS/ANAL FILE VALID
C>     AT THE PREPBUFR CENTER TIME
C> -2002-05-10 D. A. KEYSER -- ADDED "AIRCAR" TO THE LIST OF TABLE A
C>     MESSAGE TYPES THAT WILL STILL HAVE THE BACKGROUND ENCODED WHEN
C>     DOFCST IS FALSE (BECAUSE ACARS ARE NOW Q.C.'d IN PREPOBS_ACARSQC
C>     PROGRAM)
C> -2003-09-02 D. A. KEYSER -- ADDED "MSONET" TO THE LIST OF TABLE A
C>     MESSAGE TYPES THAT WILL HAVE THE VIRTUAL TEMPERATURE CALCULATED;
C>     DOES NOT CALL UFBINT FOR OUTPUTTING DATA IF "NLEV" (4'TH
C>     ARGUMENT) IS ZERO (NOW CAN ONLY HAPPEN FOR GOESND FORECAST DATA
C>     WHEN ONLY RADIANCES ARE PRESENT)
C> -2004-08-30 D. A. KEYSER -- NOW INCLUDES THE 4 LAYER PWATERS, THESE
C>     GET AN OBS. ERROR (EACH THE SAME AS TOTAL PWATER) AND AN EVENT
C>     IS GENERATED WITH A REJECTED Q.M. FOR THE 4 LAYER PWATERS IF THE
C>     PWATER OBS. ERROR READ IN IS MISSING (THIS CHANGE ALLOWS THE ETA/
C>     GSI TO PROCESS OBS. ERRORS IN THE PREPBUFR FILE THE SAME AS THE
C>     ETA/3DVAR DID WHEN READING THE OBS. ERRORS FROM AN EXTERNAL
C>     FILE); FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY OF THE REPORTED
C>     (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS
C>     NEW REASON CODE 3
C> -2004-09-10 D. T. KLEIST -- ADDED CAPABILITY TO READ GUESS FIELDS FROM
C>     EITHER HYBRID OR, AS BEFORE, SIGMA GLOBAL FORECAST FILES
C> -2005-01-03 D. A. KEYSER -- FIXED ERROR READING CDAS SGES FILE WHICH
C>     STILL HAS A 207-WORD HEADER (T62) {2004-09-10 CHANGE ASSUMED ALL
C>     SGES FILES HAD A 226-WORD HEADER (T254), BUT THIS IS VALID ONLY
C>     FOR GFS SGES)
C> -2006-05-05 R. E. TREADON -- CHANGE VERTICAL INTERPOLATION TO DIRECTLY
C>     USE PRESSURE PROFILE, NOT PRESSURE PROFILE CONVERTED TO SIGMA.  
C>     THIS CHANGE IS IN SUBROUTINE GBLEVN03.  AS A RESULT OF THIS 
C>     CHANGE, SUBROUTINE GBLEVN07 WAS REMOVED.
C> -2006-07-14 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "SOME_FCST"
C>     WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOFCST" IS FALSE: IF
C>     DOFCST=F AND SOME_FCST=T THEN, JUST AS BEFORE WHEN DOFCST=F, A
C>     FORECAST WILL STILL BE ENCODED FOR REPORTS IN CERTAIN MESSAGE
C>     TYPES USED IN SUBSEQUENT Q.C. PROGRAMS (I.E, "ADPUPA", "AIRCFT",
C>     "AIRCAR", "PROFLR" OR "VADWND") (THE DEFAULT FOR SOME_FCST IS
C>     TRUE); HOWEVER IF DOFCST=F AND SOME_FCST=F THEN A FORECAST WILL
C>     NOT BE ENCODED INTO REPORTS IN ANY MESSAGE TYPE (THIS ALLOWS
C>     THIS PROGRAM TO ENCODE OBS ERRORS AND/OR VIRTUAL TEMPERATURE
C>     EVENTS INTO A PREPBUFR FILE WITHOUT ENCODING A FORECAST); ADDED
C>     NEW NAMELIST SWITCH "ADPUPA_VIRT" WHICH, WHEN TRUE, INCLUDES
C>     REPORTS IN MESSAGE TYPE ADPUPA (I.E., RAOBS, DROPS, MULTI-LEVEL
C>     RECCOS) IN THE "VIRTMP" PROCESSING (PROCESSING THEM WITH SAME
C>     LOGIC AS IN SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR)
C>     {NORMALLY "ADPUPA_VIRT" IS FALSE (DEFAULT) BECAUSE SUBSEQUENT
C>     PROGRAM PREPOBS_CQCBUFR PERFORMS THIS FUNCTION}
C> -2007-09-14 S. MOORTHI -- ADDED CAPABILITY TO READ GENERALIZED SIGMA/
C>     HYBRID FILES FROM THE GFS USING "SIGIO" UTILITY; ALSO, CLEANED UP
C>     SOME CODE; NEW ERROR CONDITION CODES 70 AND 71 ADDED
C> -2007-09-14 D. A. KEYSER -- FUNCTION OEFG01, WHICH RETURNS THE OBS
C>     ERROR FOR A REQUESTED VARIABLE INTERP. TO A DEFINED PRESSURE
C>     LEVEL FOR A DEFINED REPORT TYPE, MODIFIED TO USE EXACT LOGIC AS
C>     IN GSI (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE
C>     TYPE, LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO
C>     MB, A FEW OTHER MINOR CHANGES) - THIS WILL ALLOW GSI TO READ OBS
C>     ERROR DIRECTLY OUT OF PREPBUFR FILE RATHER THAN OUT OF AN
C>     EXTERNAL FILE; FOR PW TYPES, NOW PASSES REPORTED SURFACE PRESSURE
C>     (PRSS * 0.01) INTO FUNCTION OEFG01 RATHER THAN VERTICAL
C>     COORDINATE PRESSURE (POB), SINCE LATTER IS ALWAYS MISSING FOR
C>     THESE TYPES (DOESN'T CHANGE VALUE COMING OUT OF OEFG01 SINCE IT
C>     IS CONSTANT ON ALL LEVELS ANYWAY FOR PW); IN SUBR. GBLEVN02, Q.M.
C>     9 IS NOW ASSIGNED TO A VARIABLE ONLY IF ITS OBS ERROR IS MISSING,
C>     OR IN THE CASE OF MOISTURE IF THE LEVEL IS ABOVE PRESSURE LEVEL
C>     "QTOP_REJ" OR IF ITS TEMPERATURE OBS ERROR IS MISSING, ALL OTHER
C>     EVENT (E.G., GROSS CHECK ERRORS) ASSIGN Q.M. 8 (EVEN IF OBS ERROR
C>     IS MISSING), PRIOR TO THIS ONLY REJECTION OF PRESSURE ON LEVEL
C>     RESULTED IN Q.M. 8, ALL OTHER REJECTIONS GOT Q.M. 9 - THIS MEANS
C>     TRULY "BAD" OBS WILL NOW ALWAYS GET Q.M. 8 AND ONLY OBS FLAGGED
C>     FOR NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M.
C>     9 (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES
C>     NOT EVEN CONSIDER OBS WITH Q.M. 8); CORRECTED ERROR WHICH
C>     MISTAKENLY ASSIGNED REASON CODE OF 9 INSTEAD OF 3 TO MOISTURE
C>     WITH MISSING OBS ERROR; IN SUBR. GBLEVN02, Q.M. 9 WILL NOT BE
C>     ASSIGNED TO A VARIABLE IF THAT VARIABLE ALREADY HAS A "BAD" Q.M.
C>     (I.E., > 3 BUT < 15), IN FACT THE "PREVENT" EVENT WHICH WOULD
C>     ASSIGN Q.M. 9 IS SKIPPED ENTIRELY (DO NOT WANT THE GSI TO MONITOR
C>     THE OBS WHICH REALLY ARE ARE "BAD"); IN SUBR. GBLEVN08, FOR NON-
C>     "ADPUPA" TYPES, Q.M. 9 IS NOW ASSIGNED TO CALCULATED VIRT. TEMPS
C>     IF THE MOISTURE Q.M.  IS 9 OR 15 AND ORIG. TEMP NOT "BAD", THESE
C>     "VIRTMP" EVENTS RECEIVE NEW REASON CODE 4, HAD RECEIVED Q.M. 8
C>     WITH REASON CODE 2 LIKE VIRT. TEMPS CALCULATED FROM "BAD"
C>     MOISTURE - THIS MEANS ONLY TRULY "BAD" VIRT. TEMPS WILL NOW GET
C>     Q.M. 8 AND VIRT. TEMPS FLAGGED FOR NON-USE BY ASSIMILATION (BUT
C>     STILL "GOOD") WILL NOW GET Q.M. 9 (GSI MONITORS, BUT DOES NOT
C>     USE, OBS WITH Q.M. 9, BUT IT DOES NOT EVEN CONSIDER OBS WITH Q.M.
C>     8); IN SUBR. GBLEVN08, FOR "ADPUPA" TYPES, Q.M. 3 IS NOW ASSIGNED
C>     TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS TRULY BAD
C>     (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM IS 1 OR 2
C>     AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN QQM WAS 9 OR
C>     15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES, ENCODES A
C>     SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP"
C>     EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO WHAT IS
C>     ALREADY DONE FOR "RASSDA" TYPES)
C> -2010-01-29 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "RECALC_Q"
C>     WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOVTMP" IS FALSE: IF
C>     DOVTMP=F AND RECALC_Q=T THEN, JUST AS BEFORE WHEN DOVTMP=F, SPEC.
C>     HUMIDITY IS STILL RE-CALCULATED AND THE EVENT IS ENCODED INTO THE
C>     PREPBUFR FILE (BUT VIRTUAL TEMP. IS NOT ENCODED) (THE DEFAULT FOR
C>     RECALC_Q IS TRUE), HOWEVER IF DOVTMP=F AND RECALC_Q=F THEN SPEC.
C>     HUMIDITY IS NOT RE-CALCULATED (AND NEITHER IS VIRTUAL
C>     TEMPERATURE) (THIS ALLOWS THIS PROGRAM TO BYPASS ALL "VIRTMP"
C>     EVENT PROCESSING); ADDED NEW NAMELIST SWITCH "DOPREV" WHICH,
C>     WHEN TRUE, WRITES "PREVENT" EVENTS INTO THE PREPBUFR FILE (IT
C>     ALWAYS DID THIS BEFORE) (DEFAULT), BUT NOW ALLOWS THE PROGRAM
C>     TO BYPASS "PREVENT" EVENT PROCESSING WHEN DOPREV=F; INITIALIZED
C>     ARRAY IDATE AS ZERO IN SUBR. GBLEVN10, CORRECTED BUG WHICH
C>     EXPOSED PREVIOUSLY HIDDEN MEMORY CLOBBERING WHEN CALLING PROGRAMS
C>     WERE LINKED TO NEW BUFRLIB; RULES IN SUBROUTINE GBLEVN02 REFINED
C>     TO INCLUDE FULL SFC PRESSURE SANITY CHECK FOR ALL SFC REPORTS
C>     (MASS, 18x, & WIND, 28x), BEFORE ONLY DONE FOR SFC MASS REPORTS
C>     (18x) AND STILL NOT DONE FOR NON-SFC WIND REPORTS SINCE LOWEST
C>     LEVEL PRESSURE NOT NECESSARILY AT THE SFC), AS A RESULT 28x WINDS
C>     WILL NOW GET QM=8 IF PRESSURE FAILS SANITY CHECK (OFTEN HAPPENS
C>     IN MESONET REPORTS) (GSI WAS ALREADY NOT USING THESE WINDS SINCE
C>     PRESSURE QM SET TO 8 ALL ALONG)
C> -2012-11-20  J. WOOLLEN  INITIAL PORT TO WCOSS. ADDED CALL TO BUFRLIB
C>     ROUTINE GETBMISS TO ADAPT BMISS TO LINUX ENVIRONMENT IF NEED BE
C>     {I.E., OBTAINS BUFRLIB MISSING (BMISS) VIA CALL TO GETBMISS
C>     RATHER THAN HARDWIRING IT TO 10E10 (10E10 CAN CAUSE INTEGER
C>     OVERFLOW ON WCOSS - SEE CALLING PROGRAM FOR MORE INFO)}
C> -2013-02-13  D. A. KEYSER -- FINAL CHANGES TO RUN ON WCOSS: USE
C>     FORMATTED PRINT STATEMENTS WHERE PREVIOUSLY UNFORMATTED PRINT WAS
C>     > 80 CHARACTERS; RENAME ALL REAL(8) VARIABLES AS *_8
C> -2013-04-12  D. A. KEYSER -- IN SUBROUTINE GBLEVN08, DON'T ALLOW
C>     CALCULATED Q TO BE < 0 WHICH CAN OCCUR ON WCOSS FOR CASES OF
C>     HORRIBLY BAD MESONET DATA
C> -2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH,
C>     WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION
C>     PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST-
C>     COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN
C>     REPORTED PMSL IS MISSING (DONE IN SUBROUTINE GBLEVN08).  DOVTMP
C>     MUST BE TRUE AND DOANLS MUST BE FALSE ("PREVENTS" MODE). THE
C>     DERIVED PMSL EITHER GETS A QUALITY MARK ("PMQ") OF 3 OR INHERITS
C>     THE STATION PRESSURE QUALITY MARK ("PQM"), WHICHEVER IS GREATER.
C>     THE VALUE OF THE PMSL INDICATOR (NEW MNEMONIC "PMIN") IS SET TO 1
C>     TO DENOTE PMSL WAS DERIVED RATHER THAN OBSERVED.  THE DEFAULT FOR
C>     "DOPMSL" IS FALSE (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS).  IT
C>     IS FORCED TO BE FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE). IN
C>     SUBROUTINE GBLEVN02, SFCSHP REPORTS WITH CALM WINDS AND NON-
C>     MISSING BACKGROUND U- OR V-COMPONENT WIND .GE. 5 M/SEC ARE
C>     FLAGGED WITH Q.M. 8 (EVENT PGM "PREVENT", REASON CODE 8).
C> -2014-05-08 JWhiting -- altered print statement (2 format) in GBLEVN10 
C>     subroutine; increased field width for spectral resolution to 
C>     accommodate models w/ up to 5-digit resolution (I3 to I5).
C> -2016-06-13 FANGLIN YANG AND RUSS TREADON -- HANG LEI ADDED NEMSIO 
C>     TO SUBROUTINE GBLEVN10 AND REMOVED ALL SIGIO CAPABILITY.
C>     THIS UPDATE RESTORES GBLEVN10 FOR PROCESSING SIGIO INPUT, AND 
C>     ADDS A NEW GBLEVN12 FOR PROCESSING NEMSIO INPUT. THE INPUT GFS 
C>     FILE TYPE SIGIO VS NEMSIO IS NOW DETERMINED IN THE MAIN PROGRAM. 
C>     THE CODE IS ALSO UPDATED TO REMOVE BUGS. SUBROUTINE SIGIO_MODPR
C>     IS USED TO COMPUTE LAYER AND INTERFACE PRESSURE FOR NEMSIO INPUT.
C> -2017-02-17 D Keyser & J Whiting -- In subroutine GBLEVN12, removed
C>     references to multiple input files, since only 1 nemsio formatted
C>     input is needed (no interpolation is attempted; c.f.; what is
C>     done with sigio formatted inputs); GBLEVN12 routine still retains
c     some structures (esp array references) left over from multi file
C>     input. Updated comments and docblock to account for new NEMSIO
C>     input.
C> -2017-02-22 D. Keyser -- Further changes to subr. GBLEVN12 to remove
C>     array and logic references to multiple input files.
C> -2019-10-31 Hang Lei -- Add GBLEVN13 to process netcdf input.
C>     
C>
C> USAGE:    CALL GBLEVENTS(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET,
C>          $               NEWTYP)
C>
C>   INPUT ARGUMENT LIST:
C> @param IDATEP   - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
C> @param IUNITF   - 2-WORD ARRAY:
C>        For SIGIO input:
C>              - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
C>                IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
C>                AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
C>                IS A MULTIPLE OF 3
C>              - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
C>                OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
C>                THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
C>        For NEMSIO input:
C>              - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
C>                (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
C>                THE DATE IN IDATEP
C>              - WORD 2 - NOT USED, SHOULD BE A NULL FILE
C> @param IUNITE   - UNIT NUMBER OF INPUT OBSERVATION ERROR FILE
C>              - (USED ONLY IN PREVENTS MODE)
C> @param IUNITP   - UNIT NUMBER OF OUTPUT PREPBUFR DATA SET
C> @param IUNITS   - UNIT NUMBER OF "PREVENT" EVENTS DATA FILTERING
C>              - SUMMARY PRINT FILE
C>              - (USED ONLY IN PREVENTS MODE)
C> @param SUBSET   - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
C>              - REPORT BEING PROCESSED
C> @param NEWTYP   - INDICATOR IF THE BUFR MESSAGE TABLE A ENTRY HAS
C>              - CHANGED FROM THAT OF THE PREVIOUS REPORT (=0 - NO,
C>              - =1 - YES)
C>                
C>
C>   INPUT FILES:
C>   -  UNIT 05  - STANDARD INPUT (DATA CARDS - SEE NAMELIST
C>                DOCUMENTATION BELOW)
C>                (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS
C>                       SUBROUTINE RUNS IN POSTEVENTS MODE)
C>   -  UNIT AA  - PREPBUFR DATA SET
C>                (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN
C>                INPUT ARGUMENT LIST)
C>   -  UNIT BB  - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE)
C>                FILE
C>                (WHERE BB IS UNIT NUMBER DEFINED AS IUNITF(1) IN
C>                INPUT ARGUMENT LIST)
C>   -  UNIT CC  - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE)
C>                FILE
C>                (WHERE CC IS UNIT NUMBER DEFINED AS IUNITF(2) IN
C>                INPUT ARGUMENT LIST)
C>                NOTE: only valid for SIGIO input
C>   -  UNIT DD  - OBSERVATION ERROR FILE (WHERE DD IS UNIT NUMBER
C>                DEFINED AS IUNITE IN INPUT ARGUMENT LIST)
C>                (USED ONLY IN PREVENTS MODE)
C>
C>   OUTPUT FILES:
C>   -  UNIT 06  - STANDARD OUTPUT PRINT
C>   -  UNIT AA  - PREPBUFR DATA SET
C>                (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN
C>                INPUT ARGUMENT LIST)
C>   -  UNIT DD  - "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE
C>                (WHERE DD IS UNIT NUMBER DEFINED AS IUNITS IN
C>                INPUT ARGUMENT LIST)
C>                (USED ONLY IN PREVENTS MODE)
C>
C>   SUBPROGRAMS CALLED:
C>     UNIQUE:      GBLEVN02  GBLEVN03  GBLEVN04  GBLEVN06  OEFG01
C>                  GBLEVN08  GBLEVN10  GBLEVN11  GBLEVN11D GBLEVN12
C>                  GBLEVN13
C>                  GETLATS
C>                  MODULES: GBLEVN_MODULE SIGIO_MODULE     SIGIO_R_MODULE
C>                           NEMSIO_MODULE NEMSIO_OPENCLOSE NEMSIO_READ
C>                           NEMSIO_WRITE
C>     LIBRARY:
C>       SIGIO    - SIGIO_RROPEN       SIGIO_RRHEAD    SIGIO_SCLOSE  SIGIO_ALDATS
C>                  SIGIO_ALDATM       SIGIO_RRDATS    SIGIO_RRDATM  SIGIO_AXDATS
C>                  SIGIO_AXDATM       SIGIO_MODPR     SIGIO_CNVTDV
C>       SPLIB    - SPTEZM             SPTEZMV         SPLAT
C>       W3NCO    - W3MOVDAT           ERREXIT
C>       BUFRLIB  - UFBINT             UFBQCD          GETBMISS      IBFMS
C>       NEMSIO   - NEMSIO_OPEN        NEMSIO_CLOSE    NEMSIO_INIT
C>                  NEMSIO_GETFILEHEAD NEMSIO_READRECV NEMSIO_FINALIZE
C>                  NEMSIO_GETHEADVAR  NEMSIO_READRECVw34
C>
C>   EXIT STATES:
C>   -  COND =   0 - SUCCESSFUL RUN
C>   -  COND =  60 - OBSERVATION ERROR TABLE EMPTY OR DOES NOT EXIST
C>   -  COND =  61 - VARIABLE NLTD .NE. VARIABLE NLEV
C>   -  COND =  62 - VARIABLE NLTQ .NE. VARIABLE NLEV
C>   -  COND =  63 - VARIABLE NLQQ .NE. VARIABLE NLEV
C>   -  COND =  68 - DATE OF FIRST GUESS/ANALYSIS FILE(S) DOES NOT MATCH,
C>                  OR AT LEAST SPAN, THE CENTER DATE FOR THE PREPBUFR
C>                  FILE
C>   -  COND =  69 - FOR SIGIO INPUT GLOBAL FILES, VARIABLE KMAX TOO BIG
C>                 - UNABLE TO TRANSFORM FIRST GUESS OR ANALYSIS FILE(S)
C>   -  COND =  70 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RROPEN
C>                  RETURNED WITH NON-ZERO R.C.
C>   -  COND =  71 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RRHEAD
C>                  RETURNED WITH NON-ZERO R.C.
C>
C>
C> REMARKS: THIS SUBROUTINE MAY NOT WORK CORRECTLY IN THE EIGHT BYTE
C>     INTEGER W3NCO (_8) LIBRARY.  PLEASE COMPILE APPLICATION CODE
C>     USING A FOUR BYTE REAL W3NCO LIBRARY (_4 OR _d).
C>
C>     THIS ROUTINE PROCESSES ONE REPORT AT A TIME.  IT EXPECTS THAT THE
C>     CALLING PROGRAM HAS ALREADY ENCODED THE REPORT INTO THE PREPBUFR
C>     FILE VIA THE UFBINT OR UFBCPY ROUTINES.  THE CALLING PROGRAM
C>     SHOULD THEN CALL THIS ROUTINE AND, UPON ITS RETURN, THE CALLING
C>     PROGRAM SHOULD CALL WRITSB TO ACTUALLY WRITE THE UPDATED SUBSET
C>     (REPORT) INTO THE BUFR MESSAGE.
C>C
C>  ***** VARIABLES IN NAMELIST PREVDATA READ IN BY THIS SUBROUTINE *****
C>               (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS
C>                      SUBROUTINE RUNS IN POSTEVENTS MODE - DOANLS=TRUE
C>                      AND ALL OTHER VARIABLES ARE SET TO FALSE)
C>C
C>C
C>  -  DOPREV  - WRITE "PREVENT" EVENT INTO THE PREPBUFR FILE?
C>              DOPREV = .TRUE.  ---> YES                       (DEFAULT)
C>              DOPREV = .FALSE. ---> NO
C>  -  DOVTMP, ADPUPA_VIRT & RECALC_Q:
C>      DOVTMP  - WRITE VIRTUAL TEMPERATURE EVENT ("VIRTMP") INTO THE
C>                PREPBUFR FILE (I.E., RE-CALCULATE SPECIFIC HUMIDITY
C>                THEN CALCULATE VIRTUAL TEMPERATURE) FOR THE FOLLOWING
C>                TYPES OF REPORTS:
C>                 ADPUPA_VIRT = .FALSE.  ---> SURFACE LAND, MARINE,
C>                                             MESONET AND RASS REPORTS?
C>                 ADPUPA_VIRT = .TRUE.   ---> SURFACE LAND, MARINE,
C>                                             MESONET RASS, RAOB, DROP
C>                                             AND MULTI-LEVEL RECCO
C>                                             REPORTS?
C>              FOR ALL TYPES EXCEPT RASS, THIS WILL ATTEMPT TO
C>              CALCULATE VIRTUAL TEMPERATURE FROM SENSIBLE TEMPERATURE
C>              AND THE JUST RE-CALCULATED SPECIFIC HUMIDITY AND ENCODE
C>              IT AS A STACKED EVENT IN THE PREPBUFR FILE.  FOR RASS
C>              REPORTS THIS WILL JUST ENCODE THE REPORTED TEMPERATURE
C>              AS A STACKED EVENT IN THE PREPBUFR FILE SINCE THE
C>              REPORTED TEMPERATURE IS ALREADY VIRTUAL (NO MOISTURE IS
C>              PRESENT SO Q IS NOT RE-CALCULATED FOR RASS REPORTS).
C>                DOVTMP = .TRUE.  ---> YES                     (DEFAULT)
C>                DOVTMP = .FALSE.
C>                  RECALC_Q = .TRUE.  ---> RE-CALCULATE SPECIFIC
C>                                           HUMIDITY BUT DO NOT THEN
C>                                           CALCULATE VIRTUAL
C>                                           TEMPERATURE        (DEFAULT)
C>                  RECALC_Q = .FALSE. ---> NO, DO NOT RE-CALCULATE
C>                                           SPECIFIC HUMIDITY AND DO
C>                                           NOT CALCULATE VIRTUAL
C>                                           TEMPERATURE
C>      {NOTE1: FOR SURFACE LAND, MARINE AND MESONET REPORTS, (AND
C>              RAOB, DROP AND MULTI-LEVEL RECCO REPORTS IF
C>              "ADPUPA_VIRT"=TRUE) DOVTMP=FALSE WILL STILL RE-CALCULATE
C>              SPECIFIC HUMIDITY AND ENCODE IT AS A STACKED EVENT IN
C>              THE PREPBUFR FILE UNLESS EITHER DOANLS IS TRUE OR
C>              RECALC_Q IS FALSE.)
C>      (NOTE2: DOES NOT APPLY TO ANY REPORT TYPES OTHER THAN THOSE
C>              MENTIONED ABOVE)
C>      (NOTE3: IF DOANLS=TRUE, THEN DOVTMP IS NOT ONLY FORCED TO BE
C>              FALSE, BUT ALSO SPECIFIC HUMIDITY IS NOT RE-CALCULATED.)
C>      (NOTE4: ADPUPA_VIRT DEFAULTS TO FALSE.)
C>      (NOTE5: IF DOVTMP=TRUE, THEN RECALC_Q IS MEANINGLESS.)
C>      (NOTE6: RECALC_Q DEFAULTS TO TRUE.)
C>
C>  -  DOFCST & SOME_FCST:
C>      DOFCST  - ENCODE FORECAST (FIRST GUESS) VALUES, INTERPOLATED FROM
C>                EITHER A SIGIO (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT
C>                GLOBAL FILE, INTO THE PREPBUFR FILE FOR ALL MESSAGE
C>                TYPES OR AT LEAST SOME MESSAGE TYPES?
C>                DOFCST = .TRUE.  ---> YES, ENCODE FORECST FOR ALL
C>                                       MESSAGE TYPES          (DEFAULT)
C>                DOFCST = .FALSE.
C>                  SOME_FCST = .FALSE. ---> NO, DO NOT ENCODE FORECAST
C>                                            FOR ANY MESSAGE TYPE
C>                                            (VALUES REMAIN MISSING)
C>                  SOME_FCST = .TRUE.  ---> YES, BUT ONLY FOR MESSAGE
C>                                            TYPES "ADPUPA", "AIRCFT",
C>                                            "AIRCAR", "PROFLR" OR
C>                                            "VADWND" (VALUES REMAIN
C>                                            MISSING FOR ALL OTHER
C>                                            MESSAGE TYPES)
C>      (NOTE1: THE CASE DOFCST=FALSE & SOME_FCST=TRUE WRITES THE
C>              FORECAST VALUES FOR THE TYPES MENTIONED ABOVE BECAUSE
C>              THEY ARE NEEDED BY SUBSEQUENT QUALITY CONTROL PROGRAMS.)
C>      (NOTE2: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE
C>              NON-GLOBAL VERSIONS SINCE ONLY THE GLOBAL REQUIRES A
C>              FIRST GUESS TO BE PRESENT FOR ALL CONVENTIONAL MESSAGE
C>              TYPES IN THE PREPBUFR FILE.)
C>      (NOTE3: IF DOANLS=TRUE, THEN DOFCST & SOME_FCST ARE FORCED TO BE
C>              FALSE, MEANING A GUESS WILL NOT BE ENCODED FOR ANY
C>              MESSAGE TYPE.)
C>      (NOTE4: IF DOFCST=TRUE, THEN SOME_FCST IS MEANINGLESS.)
C>      (NOTE5: SOME_FCST DEFAULTS TO TRUE.)
C>
C>  -  DOANLS  - ENCODE ANALYZED VALUES, INTERPOLATED FROM EITHER A SIGIO
C>    (SIGMA OR HYBRID)  INPUT OR NEMSIO INPUT GLOBAL FILE,
C>              INTO THE PREPBUFR FILE -
C>              POSTEVENTS MODE - ?
C>              DOANLS = .TRUE.  ---> YES, FOR ALL MESSAGE TYPES
C>              DOANLS = .FALSE. ---> NO,  FOR ALL MESSAGE TYPES
C>                                          - PREVENTS MODE -   (DEFAULT)
C>      (NOTE:  DOANLS=TRUE WILL OVERRIDE AND FORCE TO FALSE ALL OTHER
C>              SWITCHES. IN ADDITION, THE FORECAST VALUES WILL NOT
C>              BE ENCODED FOR ANY MESSAGE TYPE AND SPECIFIC HUMIDITY
C>              WILL NOT BE RE-CALCULATED.)
C>
C>  -  DOBERR  - ENCODE OBSERVATIONAL ERROR VALUES, AS READ FROM OBS.
C>              ERROR FILE, INTO THE PREPBUFR FILE?
C>              DOBERR = .TRUE.  ---> YES                       (DEFAULT)
C>              DOBERR = .FALSE. ---> NO (VALUES REMAIN MISSING)
C>      (NOTE1: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE
C>              RAP -AND PREVIOUS RUC- VERSION SINCE IT DOES NOT REQUIRE
C>              OBSERVATIONAL ERRORS TO BE PRESENT IN THE PREPBUFR FILE.)
C>      (NOTE2: IF DOANLS=TRUE, THEN DOBERR IS FORCED TO BE FALSE.)
C>
C>  -  QTOP_REJ  - THE PRESSURE LEVEL (IN MB) ABOVE WHICH ALL SPECIFIC
C>                HUMIDITY QUALITY MARKERS ARE "REJECTED" (THE QUALITY
C>                MARKER IS SET TO 9 ON ALL PRESSURE LEVELS LESS THAN
C>                THIS LEVEL)                              (DEFAULT=300.)
C>
C>  -  SATMQC  - PERFORM SPECIAL QUALITY CONTROL ON SATELLITE TEMPERATURE
C>              SOUNDINGS IN REPORT TYPES 160-179?
C>              SATMQC = .TRUE.  ---> YES
C>              SATMQC = .FALSE. ---> NO                        (DEFAULT)
C>      (NOTE: THIS APPLIES ONLY TO THE CDAS OR HISTORICAL RE-RUNS
C>             WITH TEMPERATURE SOUNDINGS IN THESE REPORT TYPES)
C>
C>  -  DOPMSL  - ENCODE DERIVED PMSL ("PMO") FOR ALL SURFACE REPORTS WHEN
C>              REPORTED PMSL IS MISSING - ?
C>              DOPMSL = .TRUE.  ---> YES
C>              DOPMSL = .FALSE. ---> NO ("PMO" REMAINS MISSING)(DEFAULT)
C>      {NOTE: THIS APPLIES ONLY WHEN DOVTMP=TRUE AND DOANLS=FALSE
C>             ("PREVENTS" MODE), VIRTUAL TEMPERATURE CAN BE CALCULATED,
C>             AND STATION PRESSURE AND SURFACE HEIGHT/ELEVATION ARE
C>             BOTH PRESENT.  THE DERIVED PMSL EITHER GETS A QUALITY
C>             MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE
C>             QUALITY MARK ("PQM") PQM, WHICHEVER IS GREATER. THE VALUE
C>             OF THE PMSL INDICATOR ("PMIN") IS SET TO 1 TO DENOTE PMSL
C>             WAS DERIVED RATHER THAN OBSERVED.}
C>
      SUBROUTINE gblevents(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET,
     $                     NEWTYP)
 
      USE sigio_module
      USE sigio_r_module
      USE nemsio_module
      USE nemsio_openclose
      USE nemsio_read
      USE gblevn_module
      use netcdf
 
      INTEGER, PARAMETER :: IM=384, jm=im/2+1
      integer :: idrt = 0
 
      CHARACTER*80 HEADR,OBSTR,QMSTR,FCSTR,OESTR,ANSTR
      CHARACTER*8  SUBSET
      REAL(8)      OBS_8,QMS_8,BAK_8,SID_8,HDR_8(10)
      REAL(8)      BMISS,GETBMISS
      LOGICAL      DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT,DOANLS,
     $             satmqc,adpupa_virt,recalc_q,doprev,dopmsl
      INTEGER*4 IRET,IRET1
      INTEGER*4    IUNITF(2),ncid
 
      CHARACTER*20 CFILE1
      TYPE(sigio_head) :: HEAD1
      TYPE(nemsio_gfile) :: GFILE1  
 
      COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
     $ xob,yob,dhr,typ,nlev
      COMMON /gbevbb/ pvcd,vtcd
      COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
     $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
      COMMON /gbevdd/ errs(300,33,6)
      COMMON /gbevff/ bmiss
 
      SAVE
 
      DATA ifirst/0/
 
      DATA headr /
     $ 'SID XOB YOB DHR TYP                                       '/
      DATA obstr /
     $ 'POB QOB TOB ZOB UOB VOB PWO  PW1O PW2O PW3O PW4O CAT PRSS '/
      DATA qmstr /
     $ 'PQM QQM TQM ZQM WQM PWQ PW1Q PW2Q PW3Q PW4Q NUL  NUL      '/
      DATA fcstr /
     $ 'PFC QFC TFC ZFC UFC VFC PWF  PW1F PW2F PW3F PW4F NUL      '/
      DATA anstr /
     $ 'PAN QAN TAN ZAN UAN VAN PWA  PW1A PW2A PW3A PW4A NUL      '/
      DATA oestr /
     $ 'POE QOE TOE ZOE WOE PWE PW1E PW2E PW3E PW4E NUL  NUL      '/
 
      namelist /prevdata/dovtmp,dofcst,some_fcst,doberr,doanls,
     $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
 
C----------------------------------------------------------------------
C----------------------------------------------------------------------
 
      IF(ifirst.EQ.0)  THEN
 
C -------------------------------
C FIRST TIME IN DO A FEW THINGS...
C -------------------------------
 
         ifirst = 1
         print 700
  700 FORMAT(/1x,100('#')/' =====> SUBROUTINE GBLEVENTS INVOKED FOR ',
     $ 'THE FIRST TIME - VERSION LAST UPDATED 2017-02-22'/)
 
         bmiss = getbmiss()
         print *
         print *, 'BUFRLIB value for missing passed into GBLEVENTS ',
     $    'is: ',bmiss
         print *
 
C  INITIALIZE NAMELIST SWITCHES TO DEFAULT VALUES
C  ----------------------------------------------
 
         dovtmp      = .true.
         doprev      = .true.
         recalc_q    = .true.
         dofcst      = .true.
         some_fcst   = .true.
         doberr      = .true.
         doanls      = .false.
         qtop_rej    = 300.
         satmqc      = .false.
         adpupa_virt = .false.
         dopmsl      = .false.
         READ(5,prevdata,err=101,END=102)
         GO TO 103
C-----------------------------------------------------------------------
  101    CONTINUE
 
C  ERROR READING STANDARD INPUT - THIS DEFAULTS TO POSTEVENTS MODE
C  ---------------------------------------------------------------
 
         print 7013
 7013 FORMAT(/' ##> GBLEVENTS: ERROR READING STANDARD INPUT DATA CARDS',
     $ ' -- DEFAULTS TO "POSTEVENTS" MODE'/)
         doanls = .true.
         GO TO 103
 
C-----------------------------------------------------------------------
  102    CONTINUE
 
C  STANDARD INPUT IS EMPTY - THIS DEFAULTS TO POSTEVENTS MODE
C  ----------------------------------------------------------
 
         print 7014
 7014 FORMAT(/' ##> GBLEVENTS: STANDARD INPUT DATA CARDS DO NOT ',
     $ 'EXIST -- DEFAULTS TO "POSTEVENTS" MODE'/)
         doanls = .true.
 
C-----------------------------------------------------------------------
  103    CONTINUE
         IF(doanls) THEN
            dovtmp      = .false.
            doprev      = .false.
            dofcst      = .false.
            some_fcst   = .false.
            doberr      = .false.
            adpupa_virt = .false.
            dopmsl      = .false.
         ENDIF
         IF(dovtmp)  recalc_q=.true. ! RECALC_Q must be T if DOVTMP is T
         WRITE (6,prevdata)
 
         fcst = dofcst
         virt = .false.
 
C  CHECK VALID-TIME DATE OF GUESS/ANALYSIS FILE(S) AGAINST THE CENTER
C   DATE FOR THE PREPBUFR FILE AND OBTAIN THE FIRST GUESS/ANALYSIS
C   UNLESS ALL OF DOFCST, SOME_FCST, DOANLS ARE FALSE
C  ------------------------------------------------------------------
 
         IF(.NOT.doanls)  THEN
            IF(.NOT.dofcst.AND..NOT.some_fcst)  THEN
               print 901
  901 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - FIRST GUESS NOT READ ',
     $ 'IN'/)
            ELSE
               print 701
  701 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - DATE CHECK AND ',
     $ 'TRANSFORM THE FIRST GUESS'/)
            ENDIF
         ELSE
            print 7701
 7701 FORMAT(/' --> GBLEVENTS: POSTEVENTS MODE - DATE CHECK AND ',
     $ 'TRANSFORM THE ANALYSIS'/)
         ENDIF
 
         IF(dofcst .OR. some_fcst .OR. doanls) THEN
            WRITE(cfile1,'("fort.",I2.2)') iunitf(1)
 
            iret = nf90_open(trim(cfile1),nf90_nowrite,ncid)
            if (iret == 0) then
               print *,' ===> GFS FCST/ANAL INPUT IS NETCDF'
               CALL gblevn13(iunitf,idatep,im,jm,idrt)
            else
               CALL sigio_rropen(iunitf(1),cfile1,iret)
               CALL sigio_srhead(iunitf(1),head1,iret1)
               IF(iret == 0 .AND. iret1 == 0) THEN
                  print *,' ===> GLOBAL FCST/ANAL INPUT IS SIGIO'
                  CALL sigio_sclose(iunitf(1),iret)
                  CALL gblevn10(iunitf,idatep,im,jm,idrt)
               ELSE
                  CALL nemsio_open(gfile1,trim(cfile1),'read',iret=iret)
                  IF(iret == 0) THEN
                     CALL nemsio_close(gfile1,iret=iret)
                     print *,' ===> GFS FCST/ANAL INPUT IS NEMSIO'
                     CALL gblevn12(iunitf,idatep,im,jm,idrt)
                  ENDIF
               ENDIF
            ENDIF 
         ENDIF
 
         print*,'after returning from GBLEVN10, GBLEVN12,',
     $    ' or GBLEVN13.   idrt=',idrt
 
         IF(doberr)  THEN
 
C  IF REQUESTED, READ ERROR FILES (ONLY POSSIBLE IN PREVENTS MODE)
C  ---------------------------------------------------------------
 
            print 702
  702 FORMAT(/' --> GBLEVENTS: READ ERROR FILES'/)
 
            CALL gblevn01(iunite)
 
         ELSE
 
            errs = 0
            IF(.NOT.doanls)  print 3702
 3702 FORMAT(/' --> GBLEVENTS: OBS. ERROR NOT ENCODED IN PREPBUFR ',
     $ '(BY CHOICE)'/)
 
         ENDIF
 
C  OBTAIN NECESSARY PROGRAM CODES (ONLY USED IN PREVENTS MODE)
C  -----------------------------------------------------------
 
         CALL ufbqcd(iunitp,'PREVENT',pvcd)
         CALL ufbqcd(iunitp,'VIRTMP ',vtcd)
 
         print 703
  703 FORMAT(/1x,100('#')/)
 
C  SET-UP OUTPUT "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE
C  (ONLY USED IN PREVENTS MODE)
C  ----------------------------------------------------------------
 
         IF(.NOT.doanls)  WRITE(iunits,1701) idatep
 1701 FORMAT(//130('#')//38x,'*** "PREVENT" EVENTS DATA FILTERING ',
     $ 'SUMMARY ***'/35x,'--> CENTER DATE FOR PREPBUFR FILE IS: ',i10,
     $ ' <--'//)
 
C----------------------------------------------------------------------
C----------------------------------------------------------------------
 
      ENDIF
 
      IF(.NOT.doanls)  THEN
 
         IF(newtyp.EQ.1)  WRITE(iunits,1702)  subset
 1702 FORMAT(130('-')/39x,'--> SUMMARY FOR TABLE A ENTRY "',a8,'" <--'/)
 
         IF(.NOT.dofcst .AND. some_fcst)  fcst = (subset.EQ.'ADPUPA  '
     $    .OR.subset.EQ.'PROFLR  '.OR.subset .EQ.'AIRCFT  '.OR.subset
     $    .EQ.'AIRCAR  '.OR.subset .EQ.'VADWND  ')
 
C   Will not subject ACARS reports to virtual temp. processing until
C    spec. humidity is used in production
 
ccccc    VIRT = (SUBSET.EQ.'ADPSFC  '.OR.SUBSET.EQ.'SFCSHP  '.OR.
ccccc$           SUBSET.EQ.'MSONET  '.OR.SUBSET.EQ.'AIRCAR  '.OR.
ccccc$           SUBSET.EQ.'RASSDA  '.OR.SUBSET.EQ.'SATEMP  '.OR.
ccccc$           (SUBSET.EQ.'ADPUPA  '.AND.ADPUPA_VIRT))
         virt = (recalc_q.AND.(subset.EQ.'ADPSFC  '.OR.
     $                         subset.EQ.'SFCSHP  '.OR.
     $                         subset.EQ.'MSONET  '.OR.
     $                         subset.EQ.'RASSDA  '.OR.
     $                         subset.EQ.'SATEMP  '.OR.
     $                       (subset.EQ.'ADPUPA  '.AND.adpupa_virt)))
 
 
         IF(.NOT.(fcst.OR.doberr.OR.virt.OR.doprev)) THEN
            IF(newtyp.EQ.1)  WRITE(iunits,1703)
 1703 FORMAT(/'  ==> DATA FILTERING NOT PERFORMED FOR THIS TABLE A ',
     $ 'ENTRY -- FORECAST, OBS ERROR, "VIRTMP", "PREVENT" PROCESSING ',
     $ 'NOT DONE'/)
            RETURN
         ENDIF
 
      ENDIF
 
C  READY TO RETRIEVE NECESSARY INFORMATION OUT OF THE NEXT REPORT WHICH
C  HAS BEEN "UFB" ENCODED INTO THE PREPBUFR FILE BY THE CALLING PROGRAM
C  (USE NEGATIVE UNIT NUMBER HERE SINCE FILE OPEN FOR OUTPUT)
C  (NOTE: THE CALLING PROGRAM HAS NOT YET WRIITEN THE REPORT INTO
C         THE PREPBUFR FILE VIA WRITSB!)
C  ----------------------------------------------------------------
 
      CALL ufbint(-iunitp,obs_8,13,255,nlev,obstr)
      CALL ufbint(-iunitp,qms_8,12,255,nlev,qmstr)
      CALL ufbint(-iunitp,hdr_8,10,  1,iret,headr)
      sid_8 = hdr_8(1)
      xob   = hdr_8(2)
      yob   = hdr_8(3)
      dhr   = hdr_8(4)
      typ   = hdr_8(5)
 
      IF(fcst.OR.doanls)  THEN
 
C  PREVENTS MODE:   ENCODE FIRST GUESS VALUES INTO PREPBUFR FILE
C  ------------------------------------------------------------
 
C  POSTEVENTS MODE: ENCODE ANALYSIS VALUES INTO REPORT AND RETURN TO
C                   CALLING PROGRAM TO WRITE GBL-EVENTED REPORT
C                   (SUBSET) INTO PREPBUFR FILE
C  -----------------------------------------------------------------
 
         CALL gblevn03(subset)
         IF(nlev.GT.0)  THEN
            IF(fcst)   THEN
               CALL ufbint(iunitp,bak_8,12,nlev,iret,fcstr)
            ELSE
               CALL ufbint(iunitp,bak_8,12,nlev,iret,anstr)
               RETURN
            ENDIF
         ENDIF
      ENDIF
 
C  --------------------------------------------------------------------
C  LOGIC FROM HERE ON PERTAINS ONLY TO PREVENTS MODE OF THIS SUBROUTINE
C  --------------------------------------------------------------------
 
C  ENCODE OBSERVATION ERRORS INTO REPORT
C  -------------------------------------
 
      IF(doberr)  THEN
         IF(newtyp.EQ.1)  WRITE(iunits,1710)
 1710 FORMAT(/'  ==> OBS ERROR VALUES ARE ENCODED FOR THIS TABLE A ',
     $ 'ENTRY'//'  ==> FILTERING VIA MISSING OBS ERROR TEST IS  ',
     $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES ARE ',
     $ 'PROCESSED/STORED'/)
         CALL gblevn04
         IF(nlev.GT.0)  CALL ufbint(iunitp,bak_8,12,nlev,iret,oestr)
      ELSE
         IF(newtyp.EQ.1)  WRITE(iunits,1705)
 1705 FORMAT(/'  ==> OBS ERROR VALUES NOT ENCODED FOR THIS TABLE A ',
     $ 'ENTRY'//'  ==> FILTERING VIA MISSING OBS ERROR TEST NOT ',
     $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES NOT ',
     $ 'PROCESSED/STORED'/)
      ENDIF
 
C  MAKE THE GBLEVENTS EVENTS AND ENCODE INTO REPORT
C  ------------------------------------------------
 
      IF(.NOT.fcst)  THEN
         IF(newtyp.EQ.1)  WRITE(iunits,1704)
 1704 FORMAT(/'  ==> FORECAST VALUES NOT ENCODED FOR THIS TABLE A ',
     $ 'ENTRY'//'  ==> FILTERING VIA POB VS. GESS PSFC TEST NOT ',
     $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES NOT ',
     $ 'PROCESSED/STORED'/)
      ELSE
         IF(newtyp.EQ.1)  WRITE(iunits,1708)
 1708 FORMAT(/'  ==> FORECAST VALUES ARE ENCODED FOR THIS TABLE A ',
     $ 'ENTRY'//'  ==> FILTERING VIA POB VS. GESS PSFC TEST IS  ',
     $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES ARE ',
     $ 'PROCESSED/STORED'/)
      ENDIF
 
      IF(doprev)  THEN
         IF(newtyp.EQ.1)  WRITE(iunits,1807)
 1807 FORMAT(/'  ==> "PREVENT" EVENT PROCESSING IS  PERFORMED FOR THIS',
     $ ' TABLE A ENTRY'/)
         CALL gblevn02(iunitp,iunits,newtyp,subset)
      ELSE
         IF(newtyp.EQ.1)  WRITE(iunits,1806)
 1806 FORMAT(/'  ==> "PREVENT" EVENT PROCESSING NOT PERFORMED FOR THIS',
     $ ' TABLE A ENTRY'/)
      ENDIF
 
C  MAKE THE VIRTUAL TEMPERATURE EVENTS AND ENCODE INTO REPORT
C  ----------------------------------------------------------
 
      IF(.NOT.virt)  THEN
         IF(newtyp.EQ.1)  WRITE(iunits,1706)
 1706 FORMAT(/'  ==> "VIRTMP" EVENT PROCESSING NOT PERFORMED FOR THIS ',
     $ 'TABLE A ENTRY'/)
      ELSE
         IF(newtyp.EQ.1)  WRITE(iunits,1707)
 1707 FORMAT(/'  ==> "VIRTMP" EVENT PROCESSING IS  PERFORMED FOR THIS ',
     $ 'TABLE A ENTRY'/)
         CALL gblevn08(iunitp,subset)
      ENDIF
 
C  RETURN TO CALLING PROGRAM TO WRITE GBL-EVENTED REPORT (SUBSET) INTO
C   PREPBUFR FILE
C  -------------------------------------------------------------------
 
      RETURN
 
      SUBROUTINE gblevents(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET, …
      END
C***********************************************************************
C***********************************************************************
C> Read observation error table
C> @param IUNITE   - UNIT NUMBER OF INPUT OBSERVATION ERROR FILE (USED ONLY IN PREVENTS MODE)
      SUBROUTINE gblevn01(IUNITE) ! FORMERLY SUBROUTINE ETABLE
 
      COMMON /gbevdd/ errs(300,33,6)
 
C  READ THE OBSERVATION ERROR TABLES
C  ---------------------------------
 
      rewind iunite
 
      irec = 0
 
   10 CONTINUE
      READ(iunite,'(1X,I3)',END=100) kx
      irec = irec + 1
      DO k=1,33
         READ(iunite,'(1X,6E12.5)') (errs(kx,k,m),m=1,6)
      ENDDO
      GO TO 10
 
  100 CONTINUE
      IF(irec.LE.0) THEN
         print'(" ##GBLEVENTS/GBLEVN01 - OBS. ERROR TABLE EMPTY OR ",
     $    "DOES NOT EXIST - STOP 60")'
         CALL errexit(60)
      ENDIF
 
      RETURN
 
      SUBROUTINE gblevn01(IUNITE) ! FORMERLY SUBROUTINE ETABLE …
      END
C***********************************************************************
C***********************************************************************
C> Filter data
C> @param IUNITP   - UNIT NUMBER OF OUTPUT PREPBUFR DATA SET
C> @param IUNITS   - UNIT NUMBER OF "PREVENT" EVENTS DATA FILTERING
C>              - SUMMARY PRINT FILE
C>              - (USED ONLY IN PREVENTS MODE)
C> @param SUBSET   - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
C>              - REPORT BEING PROCESSED
C> @param NEWTYP   - INDICATOR IF THE BUFR MESSAGE TABLE A ENTRY HAS
C>              - CHANGED FROM THAT OF THE PREVIOUS REPORT (=0 - NO,
C>              - =1 - YES)
C>                
      SUBROUTINE gblevn02(IUNITP,IUNITS,NEWTYP,subset)
                                         ! FORMERLY SUBROUTINE FILTAN
 
      dimension    nflgrt(100:299,12),oemin(2:6)
      CHARACTER*8  STNID,subset
      CHARACTER*40 PEVN,QEVN,TEVN,WEVN,PWVN,PW1VN,PW2VN,PW3VN,PW4VN
      REAL(8)      PEV_8(4,255),QEV_8(4,255),TEV_8(4,255),WEV_8(5,255),
     $             pwv_8(4,255),pw1v_8(4,255),pw2v_8(4,255),
     $             pw3v_8(4,255),pw4v_8(4,255),obs_8,qms_8,bak_8,sid_8,
     $             ufc_8,vfc_8
      LOGICAL      FCST,REJP_PS,REJPS,REJT,REJQ,REJW,REJPW,REJPW1,
     $             rejpw2,rejpw3,rejpw4,satmqc,satemp,soln60,sols60,
     $             moerr_p,moerr_t,adpupa_virt,doberr,dofcst,some_fcst,
     $             dovtmp,virt,recalc_q,doprev
      REAL(8)      BMISS
 
      COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
     $ xob,yob,dhr,typ,nlev
      COMMON /gbevbb/ pvcd,vtcd
      COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
     $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
      COMMON /gbevee/psg01,zsg01,tg01(500),ug01(500),vg01(500),
     x     qg01(500),zint(500),pint(500),pintlog(500),plev(500),
     x     plevlog(500)
      COMMON /gbevff/ bmiss
 
      equivalence(sid_8,stnid)
 
      DATA pevn  /'POB  PQM  PPC  PRC      '/
      DATA qevn  /'QOB  QQM  QPC  QRC      '/
      DATA tevn  /'TOB  TQM  TPC  TRC      '/
      DATA wevn  /'UOB  VOB  WQM  WPC  WRC '/
      DATA pwvn  /'PWO  PWQ  PWP  PWR      '/
      DATA pw1vn /'PW1O PW1Q PW1P PW1R     '/
      DATA pw2vn /'PW2O PW2Q PW2P PW2R     '/
      DATA pw3vn /'PW3O PW3Q PW3P PW3R     '/
      DATA pw4vn /'PW4O PW4Q PW4P PW4R     '/
 
      DATA nflgrt/2400*0/
 
      DATA oemin /0.5,0.1,1.0,0.5,1.0/
 
      ni  = mod((nint(typ)/10),10)
 
      IF(newtyp.EQ.1)  nflgrt = 0
 
C  LOGICAL SWITCHES FOR OBSERVATION LOCATION FILTERING
C  ---------------------------------------------------
 
      satemp = ((typ.GE.160.AND.typ.LE.179).AND.satmqc)
      soln60 = ((typ.GE.160.AND.typ.LE.163).AND.yob.GE.-60.AND.satmqc)
      sols60 = ((typ.EQ.160.OR.typ.EQ.162.OR.typ.EQ.163).AND.yob.LT.-60
     $ .AND.satmqc)
 
C  CLEAR THE EVENT ARRAYS
C  ----------------------
 
      pev_8  = bmiss
      qev_8  = bmiss
      tev_8  = bmiss
      wev_8  = bmiss
      pwv_8  = bmiss
      pw1v_8 = bmiss
      pw2v_8 = bmiss
      pw3v_8 = bmiss
      pw4v_8 = bmiss
 
      maxpev  = 0
      maxqev  = 0
      maxtev  = 0
      maxwev  = 0
      maxpwv  = 0
      maxpw1v = 0
      maxpw2v = 0
      maxpw3v = 0
      maxpw4v = 0
 
C  LOOP OVER LEVELS APPLYING UNDERGROUND FILTERING AND SPECIAL RULES
C  -----------------------------------------------------------------
 
      IF(nlev.GT.0)  THEN
      DO l=1,nlev
 
         pob  = obs_8( 1,l)
         qob  = obs_8( 2,l)
         tob  = obs_8( 3,l)
         uob  = obs_8( 5,l)
         vob  = obs_8( 6,l)
         pwo  = obs_8( 7,l)
         pw1o = obs_8( 8,l)
         pw2o = obs_8( 9,l)
         pw3o = obs_8(10,l)
         pw4o = obs_8(11,l)
         cat  = obs_8(12,l)
         prss = obs_8(13,l)
 
         pqm  = qms_8( 1,l)
         qqm  = qms_8( 2,l)
         tqm  = qms_8( 3,l)
         zqm  = qms_8( 4,l)
         wqm  = qms_8( 5,l)
         pwq  = qms_8( 6,l)
         pw1q = qms_8( 7,l)
         pw2q = qms_8( 8,l)
         pw3q = qms_8( 9,l)
         pw4q = qms_8(10,l)
 
         rejp_ps = .false.
         moerr_p = .false.
         moerr_t = .false.
         rcd = 99999
 
C  -------------------------------------------------------------------
C  RULES FOR PRESSURE (ON ANY LEVEL) -- ALL DATA (MASS AND WIND) ON
C   LEVEL REJECTED IF:
C    - PRESSURE MORE THAN 100 MB BELOW MODEL (GUESS) SURFACE PRESSURE
C       (AND SWITCH FCST=TRUE) -- "PREVENT" PGM REASON CODE 1
C    - PRESSURE IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON CODE 2
C    - SURFACE (MASS OR WIND) REPORT PRESSURE IS REPORTED ABOVE 450 MB
C       OR BELOW 1100 MB -- "PREVENT" PGM REASON CODE 2
C  REJECTION MEANS Q.M. SET TO 8
C  -------------------------------------------------------------------
 
         IF(pob.LT.bmiss)  THEN
            IF(.NOT.fcst)  psg01 = pob
            IF(pob-psg01.GE.100. .OR. pob.LE.0. .OR.
     $       ((pob.LE.450..OR.pob.GE.1100.) .AND. ni.EQ.8)) THEN
               IF(pob.LE.0..OR.pob.LE.450..OR.pob.GE.1100.) THEN
                  IF(ni.EQ.8)  THEN
                     WRITE(iunits,302) stnid,nint(typ),yob,xob,pob
  302 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, FAILS SANITY ',
     $ 'CHECK')
                  ELSE
                     WRITE(iunits,101) stnid,nint(typ),yob,xob,pob
  101 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $'E, REJECT ALL DATA ON LVL - POB=',f6.1,' MB, FAILS SANITY CHECK')
                  ENDIF
                  rcd = 2
               ELSE
                  IF(ni.EQ.8)  THEN
                     WRITE(iunits,303) stnid,nint(typ),yob,xob,pob,psg01
  303 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, > 100 MB ',
     $ 'BELOW GES PSFC(=',f6.1,'MB)')
                  ELSE
                     WRITE(iunits,102) stnid,nint(typ),yob,xob,pob,psg01
  102 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT ALL DATA ON LVL - POB=',f6.1,' MB, > 100 MB BELOW ',
     $ 'GES PSFC(=',f6.1,' MB)')
                  ENDIF
                  rcd = 1
               ENDIF
               rej = 8
               rejp_ps = .true.
               pev_8(1,l) = pob
               pev_8(2,l) = rej
               pev_8(3,l) = pvcd
               pev_8(4,l) = rcd
               maxpev = l
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR SURFACE PRESSURE -- ALL MASS DATA ON SURFACE LEVEL
C   REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C    - PRESSURE IS MORE THAN 100 MB ABOVE OR BELOW MODEL (GUESS)
C       SURFACE PRESSURE (AND SWITCH FCST=TRUE) -- 
C       "PREVENT" PGM REASON CODE 4
C    - PRESSURE IS REPORTED ABOVE 450 MB OR BELOW 1100 MB -- "PREVENT"
C       PGM REASON CODE 2 (NOTE: DOES NOT APPLY TO SURFACE REPORTS,
C       THESE WERE TESTED FOR THIS CRITERION IN ABOVE PRESSURE TEST)
C    - PRESSURE VIOLATES RULES FOR PRESSURE ON ANY LEVEL (SEE ABOVE)
C  REJECTION FOR FIRST RULE MEANS Q.M. SET TO 9 UNLESS:
C      - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pob.LT.bmiss .AND. cat.EQ.0) THEN
            IF(.NOT.fcst)  psg01 = pob
            rejps = oefg01(pob,typ,5,oemin(5)).GE.bmiss   .OR.
     $              abs(pob-psg01).GE.100.                .OR.
     $              pob.LE.450.                           .OR.
     $              pob.GE.1100.
            IF(rejps.OR.rejp_ps) THEN
               rej = 8
               IF(.NOT.rejp_ps)  THEN
                  IF(abs(pob-psg01).GE.100.)  THEN
                     WRITE(iunits,104) stnid,nint(typ),yob,xob,pob,psg01
  104 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, > 100 MB ',
     $ 'ABOVE GES PSFC(=',f6.1,'MB)')
                     rcd = 4
                  ELSE IF(pob.LE.450..OR.pob.GE.1100.)  THEN
                     WRITE(iunits,105) stnid,nint(typ),yob,xob,pob
  105 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, FAILS SANITY ',
     $ 'CHECK - this should never be printed since test now made in ',
     $ 'section above ')
                     rcd = 2
                  ELSE
                     IF(nflgrt(nint(typ),1).EQ.0)  THEN
                        WRITE(iunits,201) nint(typ)
  201 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'ALL DATA ON SURFACE LEVEL DUE TO MISSING SFC-P OBS ERROR'/)
                        nflgrt(nint(typ),1) = 1
                     ENDIF
CDAK CDAK CDAK CDAK  WRITE(IUNITS,103) STNID,NINT(TYP),YOB,XOB,POB
CD103 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
CDAK $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,'MB, MISSING OBS.',
CDAK $ ' ERROR')
                     rcd = 3
                     rej = 9
                  ENDIF
               ENDIF
               rejp_ps = .true.
               IF(rcd.EQ.3)  moerr_p = .true.
               IF(rej.EQ.9.AND.(pqm.GT.3.AND.pqm.LT.15))  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1401) STNID,NINT(TYP),YOB,XOB,PQM
 1401 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PQM =',f4.0,' -- DO NOT APPLY PSFC QM=9 EVENT')
               ELSE
                  pev_8(1,l) = pob
                  pev_8(2,l) = rej
                  pev_8(3,l) = pvcd
                  pev_8(4,l) = rcd
                  maxpev = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR TEMPERATURE -- TOB AND QOB ON LEVEL REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C    - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
C       MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
C    - REPORT IS TYPE 160-163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE
C       SOUNDINGS, ALL PATHS), AND IS AT OR NORTH OF 60 DEGREES SOUTH
C       LATITUDE, AND PRESSURE ON LEVEL IS AT OR BELOW 100 MB (AND
C       SWITCH SATMQC=TRUE)  -- "PREVENT" PGM REASON CODE 6
C    - REPORT IS TYPE 160,162,163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE
C       SOUNDINGS, ALL PATHS BUT CLEAR), AND IS SOUTH OF 60 DEGREES
C       SOUTH LATITUDE, AND PRESSURE ON LEVEL IS BELOW 100 MB (AND
C       SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 6
C    - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
C       (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
C       2 ABOVE) (SEE ABOVE)
C    - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
C  REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS:
C      - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(tob.LT.bmiss) THEN
            rejt = oefg01(pob,typ,2,oemin(2)).GE.bmiss  .OR.
     $             (soln60.AND.nint(pob*10.).GE.1000)   .OR.
     $             (sols60.AND.nint(pob*10.).GT.1000)
            IF(rejt.OR.rejp_ps) THEN
               rej = 8
               IF(.NOT.rejp_ps)  THEN
                  IF(soln60.AND.nint(pob*10.).GE.1000)  THEN
                     IF(nflgrt(nint(typ),6).EQ.0)  THEN
                        WRITE(iunits,7304) nint(typ)
 7304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $'TOB/QOB AT AND BELOW 100 MB IF REPORT IS NORTH OF 60S LATITUDE'/)
                        nflgrt(nint(typ),6) = 1
                     ENDIF
                     rcd = 6
                  ELSE  IF(sols60.AND.nint(pob*10.).GT.1000)  THEN
                     IF(nflgrt(nint(typ),7).EQ.0)  THEN
                        WRITE(iunits,7305) nint(typ)
 7305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'TOB/QOB BELOW 100 MB IF REPORT IS SOUTH OF 60S LATITUDE'/)
                        nflgrt(nint(typ),7) = 1
                     ENDIF
                     rcd = 6
                  ELSE
                     IF(nflgrt(nint(typ),2).EQ.0)  THEN
                        IF(ni.EQ.8)  THEN
                           WRITE(iunits,304) nint(typ)
  304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'TOB/QOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
                        ELSE
                           WRITE(iunits,202) nint(typ)
  202 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'TOB/QOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ',
     $ 'MISSING OBS ERROR'/)
                        ENDIF
                        nflgrt(nint(typ),2) = 1
cdak cdak cdak cdak cdakWRITE(IUNITS,106) STNID,NINT(TYP),YOB,XOB,TOB
cd106 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $ 'E, REJECT TOB/QOB ON LVL - TOB=',F5.1,'C, MISSING OBS. ERROR')
                     ENDIF
                     rcd = 3
                     rej = 9
                  ENDIF
               ELSE
                  IF(moerr_p)  THEN
                     rcd = 3
                     rej = 9
                  ENDIF
               ENDIF
               IF(rcd.EQ.3)  moerr_t = .true.
               IF(rej.EQ.9.AND.(tqm.GT.3.AND.tqm.LT.15))  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1402) STNID,NINT(TYP),YOB,XOB,TQM
 1402 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT TQM =',f4.0,' -- DO NOT APPLY TEMP QM=9 EVENT')
               ELSE
                  tev_8(1,l) = tob
                  tev_8(2,l) = rej
                  tev_8(3,l) = pvcd
                  tev_8(4,l) = rcd
                  maxtev = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR SPECIFIC HUMIDITY -- QOB ON LEVEL REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C    - PRESSURE ON LEVEL IS ABOVE "QTOP_REJ" MB {WHERE QTOP_REJ IS
C       READ IN FROM NAMELIST "PREVDATA" (SEE DOCBLOCK)} -- "PREVENT"
C       PGM REASON CODE 5
C    - OBSERVATION ERROR FOR TEMPERATURE ON LEVEL IS MISSING (AND
C       SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
C    - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
C       MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
C    - TEMPERATURE ON LEVEL IS MISSING OR IS LESS THAN -150 DEG. C --
C       "PREVENT" PGM REASON CODE 2
C    - SPECIFIC HUMIDITY IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON
C       CODE 2
C    - REPORT IS TYPE 160-179 (SATELLITE TEMPERATURE SOUNDINGS, ALL
C       TYPES, ALL PATHS, LAND AND SEA), ALL PRESSURE LEVELS (AND
C       SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 7
C    - TEMPERATURE ON LEVEL VIOLATES RULES FOR TEMPERATURE (EXCEPT FOR
C       MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE 2 ABOVE)
C       (SEE ABOVE)
C    - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
C       (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
C       4 ABOVE) (SEE ABOVE)
C    - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
C  REJECTION FOR FIRST FOUR RULES MEANS Q.M. SET TO 9 UNLESS:
C      - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(qob.LT.bmiss) THEN
            rejq = oefg01(pob,typ,3,oemin(3)).GE.bmiss        .OR.
     $             tob.GE.bmiss                               .OR.
     $             tob.LE.-150.                               .OR.
     $             nint(pob * 10.).LT.nint(qtop_rej * 10.)    .OR.
     $             qob.LE.0.                                  .OR.
     $             satemp                                     .OR.
     $             rejt
            IF(rejq.OR.rejp_ps) THEN
               rej = 8
               IF(.NOT.rejp_ps.AND..NOT.rejt)  THEN
                  IF(satemp)  THEN
                     IF(nflgrt(nint(typ),8).EQ.0)  THEN
                        WRITE(iunits,7306) nint(typ)
 7306 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'QOB ON ALL LEVELS'/)
                        nflgrt(nint(typ),8) = 1
                     ENDIF
                     rcd = 7
                  ELSE IF(qob.LE.0..OR.tob.GE.bmiss.OR.tob.LE.-150.)THEN
                     WRITE(iunits,111) stnid,nint(typ),yob,xob,
     $                qob/1000.,tob
  111 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT QOB ON LVL - QOB=',f6.3,' G/KG, FAILS SANITY CHECK ',
     $ '(TOB=',f5.1,' C)')
                     rcd = 2
                  ELSE IF(oefg01(pob,typ,3,oemin(3)).GE.bmiss)  THEN
                     IF(nflgrt(nint(typ),3).EQ.0)  THEN
                        IF(ni.EQ.8)  THEN
                           WRITE(iunits,305) nint(typ)
  305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'QOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
                        ELSE
                           WRITE(iunits,203) nint(typ)
  203 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'QOB ON AT LEAST ONE LEVEL (IF AVAILABLE ON THAT LEVEL) DUE TO ',
     $ 'MISSING OBS ERROR'/)
                        ENDIF
                        nflgrt(nint(typ),3) = 1
                     ENDIF
                     rcd = 3
                     rej = 9
cdak cdak cdak cdak  WRITE(IUNITS,108) STNID,NINT(TYP),YOB,XOB,QOB/1000.
cd108 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT QOB ON LVL - QOB=',F6.3,'G/KG, MISSING OBS. ERROR')
                  ELSE
                     WRITE(iunits,109) stnid,nint(typ),yob,xob,
     $                qob/1000.,qtop_rej,pob
  109 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, REJECT QOB ON LVL - QOB=',f6.3,' G/KG, ABOVE ',f6.1,
     $ 'MB (POB=',f6.1,' MB)')
                     rcd = 5
                     rej = 9
                  ENDIF
               ELSE
                  IF(moerr_p.OR.moerr_t)  THEN
                     rcd = 3
                     rej = 9
                  ENDIF
               ENDIF
               IF(rej.EQ.9.AND.(qqm.GT.3.AND.qqm.LT.15))  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1403) STNID,NINT(TYP),YOB,XOB,QQM
 1403 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT QQM =',f4.0,' -- DO NOT APPLY MSTR QM=9 EVENT')
               ELSE
                  qev_8(1,l) = qob
                  qev_8(2,l) = rej
                  qev_8(3,l) = pvcd
                  qev_8(4,l) = rcd
                  maxqev = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR WIND -- UOB AND VOB ON LEVEL REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C    - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
C       MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
C       (Note: This can currently never happen because earlier check
C              for missing obs error for sfc pressure is only done if
C              "surface" level is category 0 and this is not possible
C              for wind reports.)
C    - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
C       (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
C       2 ABOVE) (SEE ABOVE)
C    - This is a SFCSHP report with calm winds and non-missing
C      background u- or v-component wind is .GE. 5 m/sec --
C      "PREVENT" PGM REASON CODE 8, Q.M. set to 8
C    - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
C  REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS:
C      - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(min(uob,vob).LT.bmiss) THEN
            rejw = oefg01(pob,typ,4,oemin(4)).GE.bmiss
            IF(rejw.OR.rejp_ps) THEN
               rej = 8
               IF(.NOT.rejp_ps)  THEN
                  IF(nflgrt(nint(typ),4).EQ.0)  THEN
                     IF(ni.EQ.8)  THEN
                        WRITE(iunits,1304) nint(typ)
 1304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'UOB/VOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
                     ELSE
                        WRITE(iunits,204) nint(typ)
  204 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'UOB/VOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ',
     $ 'MISSING OBS ERROR'/)
                     ENDIF
                     nflgrt(nint(typ),4) = 1
                  ENDIF
cdak cdak cdak    WRITE(IUNITS,112) STNID,NINT(TYP),YOB,XOB
cd112 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT UOB/VOB ON LVL - MISSING OBS. ERROR')
                  rcd = 3
                  rej = 9
               ELSE
                  IF(moerr_p)  THEN  ! This currently can never be TRUE
                                     ! since CAT is never 0 for "sfc"
                                     ! level of wind reports
                     rcd = 3
                     rej = 9
                  ENDIF
               ENDIF
               IF(rej.EQ.9.AND.(wqm.GT.3.AND.wqm.LT.15))  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1404) STNID,NINT(TYP),YOB,XOB,WQM
 1404 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT WQM =',f4.0,' -- DO NOT APPLY WIND QM=9 EVENT')
               ELSE
                  wev_8(1,l) = uob
                  wev_8(2,l) = vob
                  wev_8(3,l) = rej
                  wev_8(4,l) = pvcd
                  wev_8(5,l) = rcd
                  maxwev = l
               ENDIF
            ENDIF
c  Reject calm sfc marine winds if background u- or v- wind .ge. 5 m/s
            if(subset.eq."SFCSHP".and.uob.eq.0..and.vob.eq.0.) then
              call ufbint(-iunitp,ufc_8,1,1,iret,'UFC')
              call ufbint(-iunitp,vfc_8,1,1,iret,'VFC')
cccccccc      if(ibfms(ufc_8).eq.0.or.ibfms(vfc_8).eq.0) then
                          ! DAK: changed to only allow this test if
                          !      both UFC and VFC are non-missing
              if(ibfms(ufc_8).eq.0.and.ibfms(vfc_8).eq.0) then
                if(abs(ufc_8).ge.5..or.abs(vfc_8).ge.5.) then
                  if(wqm.le.3) then
                    write(iunits,1504) stnid,nint(typ),ufc_8,vfc_8
 1504 FORMAT(/' --> ID ',a8,', (RTP ',i3,' UFC=',f5.2,' VFC=',f5.2,') ',
     $ 'NEW WIND EVENT, UOB/VOB CALM (0 M/S) WHILE |UFC| AND/OR |VFC| ',
     $ '>= 5 M/S, QM SET TO 8'/)
                    wev_8(1,1) = uob
                    wev_8(2,1) = vob
                    wev_8(3,1) = 8
                    wev_8(4,1) = pvcd
                    wev_8(5,1) = 8
                    maxwev = 1
                  end if
                end if
              end if
            end if
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR TOTAL COLUMN PRECIPITABLE WATER -- PWO REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C  REJECTION MEANS Q.M. SET TO 9 UNLESS:
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pwo.LT.bmiss) THEN
            rejpw = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
            IF(rejpw) THEN
               IF(nflgrt(nint(typ),5).EQ.0)  THEN
                  WRITE(iunits,205) nint(typ)
  205 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'PWO DUE TO MISSING OBS ERROR'/)
                  nflgrt(nint(typ),5) = 1
               ENDIF
cdakcdakcdak   WRITE(IUNITS,113) STNID,NINT(TYP),YOB,XOB,PWO
cd113 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT PWO ON LVL - PWO=',F5.1,'MM, MISSING OBS. ERROR')
               IF(pwq.GT.3.AND.pwq.LT.15)  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1405) STNID,NINT(TYP),YOB,XOB,PWQ
 1405 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PWQ =',f4.0,' -- DO NOT APPLY PWtO QM=9 EVENT')
               ELSE
                  pwv_8(1,l) = pwo
                  pwv_8(2,l) = 9
                  pwv_8(3,l) = pvcd
                  pwv_8(4,l) = 3
                  maxpwv = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR LAYER 1 PRECIPITABLE WATER -- PW1O REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C  REJECTION MEANS Q.M. SET TO 9 UNLESS:
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pw1o.LT.bmiss) THEN
            rejpw1 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
            IF(rejpw1) THEN
               IF(nflgrt(nint(typ),9).EQ.0)  THEN
                  WRITE(iunits,206) nint(typ)
  206 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'PW1O DUE TO MISSING OBS ERROR'/)
                  nflgrt(nint(typ),9) = 1
               ENDIF
cdakcdakcdak   WRITE(IUNITS,114) STNID,NINT(TYP),YOB,XOB,PW1O
cd114 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT PW1O ON LVL - PW1O=',F5.1,'MM, MISSING OBS. ERROR')
               IF(pw1q.GT.3.AND.pw1q.LT.15)  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1406) STNID,NINT(TYP),YOB,XOB,PW1Q
 1406 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PW1Q =',f4.0,' -- DO NOT APPLY PW1O QM=9 EVENT')
               ELSE
                  pw1v_8(1,l) = pw1o
                  pw1v_8(2,l) = 9
                  pw1v_8(3,l) = pvcd
                  pw1v_8(4,l) = 3
                  maxpw1v = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR LAYER 2 PRECIPITABLE WATER -- PW2O REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C  REJECTION MEANS Q.M. SET TO 9 UNLESS:
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pw2o.LT.bmiss) THEN
            rejpw2 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
            IF(rejpw2) THEN
               IF(nflgrt(nint(typ),10).EQ.0)  THEN
                  WRITE(iunits,207) nint(typ)
  207 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'PW2O DUE TO MISSING OBS ERROR'/)
                  nflgrt(nint(typ),10) = 1
               ENDIF
cdakcdakcdak   WRITE(IUNITS,115) STNID,NINT(TYP),YOB,XOB,PW2O
cd115 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT PW2O ON LVL - PW2O=',F5.1,'MM, MISSING OBS. ERROR')
               IF(pw2q.GT.3.AND.pw2q.LT.15)  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1407) STNID,NINT(TYP),YOB,XOB,PW2Q
 1407 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PW2Q =',f4.0,' -- DO NOT APPLY PW2O QM=9 EVENT')
               ELSE
                  pw2v_8(1,l) = pw2o
                  pw2v_8(2,l) = 9
                  pw2v_8(3,l) = pvcd
                  pw2v_8(4,l) = 3
                  maxpw2v = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR LAYER 3 PRECIPITABLE WATER -- PW3O REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C  REJECTION MEANS Q.M. SET TO 9 UNLESS:
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pw3o.LT.bmiss) THEN
            rejpw3 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
            IF(rejpw3) THEN
               IF(nflgrt(nint(typ),11).EQ.0)  THEN
                  WRITE(iunits,208) nint(typ)
  208 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'PW3O DUE TO MISSING OBS ERROR'/)
                  nflgrt(nint(typ),11) = 1
               ENDIF
cdakcdakcdak   WRITE(IUNITS,116) STNID,NINT(TYP),YOB,XOB,PW3O
cd116 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT PW3O ON LVL - PW3O=',F5.1,'MM, MISSING OBS. ERROR')
               IF(pw3q.GT.3.AND.pw3q.LT.15)  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1408) STNID,NINT(TYP),YOB,XOB,PW3Q
 1408 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PW3Q =',f4.0,' -- DO NOT APPLY PW3O QM=9 EVENT')
               ELSE
                  pw3v_8(1,l) = pw3o
                  pw3v_8(2,l) = 9
                  pw3v_8(3,l) = pvcd
                  pw3v_8(4,l) = 3
                  maxpw3v = l
               ENDIF
            ENDIF
         ENDIF
 
C  -------------------------------------------------------------------
C  RULES FOR LAYER 4 PRECIPITABLE WATER -- PW4O REJECTED IF:
C    - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
C       "PREVENT" PGM REASON CODE 3
C  REJECTION MEANS Q.M. SET TO 9 UNLESS:
C      - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
C        AS IS
C  -------------------------------------------------------------------
 
         IF(pw4o.LT.bmiss) THEN
            rejpw4 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
            IF(rejpw4) THEN
               IF(nflgrt(nint(typ),12).EQ.0)  THEN
                  WRITE(iunits,209) nint(typ)
  209 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
     $ 'PW4O DUE TO MISSING OBS ERROR'/)
                  nflgrt(nint(typ),12) = 1
               ENDIF
cdakcdakcdak   WRITE(IUNITS,117) STNID,NINT(TYP),YOB,XOB,PW4O
cd117 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
cdak $'E, REJECT PW4O ON LVL - PW4O=',F5.1,'MM, MISSING OBS. ERROR')
               IF(pw4q.GT.3.AND.pw4q.LT.15)  THEN
CDAKCDAKCDAKCDAK  WRITE(IUNITS,1409) STNID,NINT(TYP),YOB,XOB,PW4Q
 1409 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
     $ 'E, INPUT PW4Q =',f4.0,' -- DO NOT APPLY PW4O QM=9 EVENT')
               ELSE
                  pw4v_8(1,l) = pw4o
                  pw4v_8(2,l) = 9
                  pw4v_8(3,l) = pvcd
                  pw4v_8(4,l) = 3
                  maxpw4v = l
               ENDIF
            ENDIF
         ENDIF
 
      ENDDO
      ENDIF
 
C  APPLY THE PROPER EVENTS
C  -----------------------
 
      IF(maxpev .GT.0) CALL ufbint(iunitp,pev_8, 4,maxpev, iret,pevn)
      IF(maxqev .GT.0) CALL ufbint(iunitp,qev_8, 4,maxqev, iret,qevn)
      IF(maxtev .GT.0) CALL ufbint(iunitp,tev_8, 4,maxtev, iret,tevn)
      IF(maxwev .GT.0) CALL ufbint(iunitp,wev_8, 5,maxwev, iret,wevn)
      IF(maxpwv .GT.0) CALL ufbint(iunitp,pwv_8, 4,maxpwv, iret,pwvn)
      IF(maxpw1v.GT.0) CALL ufbint(iunitp,pw1v_8,4,maxpw1v,iret,pw1vn)
      IF(maxpw2v.GT.0) CALL ufbint(iunitp,pw2v_8,4,maxpw2v,iret,pw2vn)
      IF(maxpw3v.GT.0) CALL ufbint(iunitp,pw3v_8,4,maxpw3v,iret,pw3vn)
      IF(maxpw4v.GT.0) CALL ufbint(iunitp,pw4v_8,4,maxpw4v,iret,pw4vn)
 
      RETURN
      SUBROUTINE gblevn02(IUNITP,IUNITS,NEWTYP,subset) …
      END
C***********************************************************************
C***********************************************************************
C>  GBLEVN03 - INTERPOLATE MODEL DATA (FIRST GUESS OR ANALYSIS) TO OB LOCATIONS
C> @param SUBSET   - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR REPORT BEING PROCESSED
C-----------------------------------------------------------------------
      SUBROUTINE gblevn03(SUBSET) ! FORMERLY SUBROUTINE GETFC
 
      USE gblevn_module
 
      REAL(8)   OBS_8,QMS_8,BAK_8,SID_8
      REAL(8)   BMISS
      CHARACTER*8  SUBSET
 
      COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
     $ xob,yob,dhr,typ,nlev
      COMMON /gbevee/psg01,zsg01,tg01(500),ug01(500),vg01(500),
     x     qg01(500),zint(500),pint(500),pintlog(500),plev(500),
     x     plevlog(500)
      COMMON /gbevff/ bmiss
 
 
      DATA tzero / 273.15 /
      DATA betap / .0552  /
      DATA beta  / .00650 /
      DATA rog   / 29.261 /
 
C  CLEAR THE BACKGROUND EVENT ARRAY
C  --------------------------------
 
      bak_8 = bmiss
 
C  GET GUESS PROFILE AT OB LOCATION
C  --------------------------------
      CALL gblevn06(xob,yob)
 
      
C  INTERPOLATE GUESS PROFILES TO OB PRESSURES
C  ------------------------------------------
 
      IF(nlev.GT.0)  THEN
      DO 10 l=1,nlev
 
         pob  = obs_8( 1,l)
         qob  = obs_8( 2,l)
         tob  = obs_8( 3,l)
         zob  = obs_8( 4,l)
         uob  = obs_8( 5,l)
         vob  = obs_8( 6,l)
         pwo  = obs_8( 7,l)
         pw1o = obs_8( 8,l)
         pw2o = obs_8( 9,l)
         pw3o = obs_8(10,l)
         pw4o = obs_8(11,l)
         cat  = obs_8(12,l)
 
         IF(pob.LE.0. .OR. pob.GE.bmiss) GOTO 10
 
         poblog = log(pob)
         
         la = -999
         lb = -999
         do k=1,kmax-1
            if (poblog<=plevlog(k) .and. poblog>plevlog(k+1)) then
               la = k
               lb = k+1
               exit
            endif
         enddo
         if (la > 0) then
            wt = (poblog-plevlog(lb)) / (plevlog(la)-plevlog(lb))
         else
            la = 1
            lb = la+1
            wt = 0.0
         endif
         
         li=0
         do k=1,kmax-1
            if (poblog<=pintlog(k) .and. poblog>pintlog(k+1)) then
               li = k
               exit
            endif
         enddo
 
C  SURFACE PRESSURE
C  ----------------
 
         IF(cat.EQ.0 .AND. zob.LT.bmiss) THEN
            ts = tg01(1) + (psg01-plev(1))*betap
            dz  = zob-zsg01
            tm  = ts - dz*beta*.5
            pfc = psg01*exp(-dz/(tm*rog))
         ELSE
            pfc = bmiss
         ENDIF
 
C  SPECIFIC HUMIDITY
C  -----------------
 
         IF(qob.LT.bmiss.OR.tob.LT.bmiss.OR.typ.EQ.111) THEN
 
C  (QFC NEEDED BY SYNDATA PROGRAM BUT ONLY FOR REPORT TYPE 111)
 
            qob = qg01(lb) + (qg01(la)-qg01(lb))*wt
         ENDIF
 
C  TEMPERATURE
C  -----------
 
         IF(tob.LT.bmiss.OR.subset.EQ.'VADWND  '.OR.typ.EQ.111) THEN
 
C  (TFC NEEDED BY CQCVAD AND SYNDATA PROGRAMS, LATTER ONLY FOR REPORT
C   TYPE 111)
 
            IF(pob.GT.plev(1)) THEN
               tob = tg01(1) + (pob-plev(1))*betap
            ELSE
               tob = tg01(lb) + (tg01(la)-tg01(lb))*wt
            ENDIF
            tob = tob - tzero
         ENDIF
 
C  HEIGHT
C  ------
 
         IF(zob.LT.bmiss) THEN
            IF(pob.GT.plev(1)) THEN
               tm = tg01(1) + (.5*(pint(1)+pob)-plev(1))*betap
               zob = zint(1) - rog*tm*log(pob/pint(1))
            ELSE
               tm = tg01(lb) + (tg01(la)-tg01(lb))*wt
               zob = zint(li) - rog*tm*log(pob/pint(li))
            ENDIF
         ENDIF
 
C  U AND V COMPONENTS
C  ------------------
 
         IF(uob.LT.bmiss .OR. vob.LT.bmiss) THEN
            uob = ug01(lb) + (ug01(la)-ug01(lb))*wt
            vob = vg01(lb) + (vg01(la)-vg01(lb))*wt
         ENDIF
 
 
C  PRECIPITABLE WATER
C  ------------------
 
         pwo  = bmiss
         pw1o = bmiss
         pw2o = bmiss
         pw3o = bmiss
         pw4o = bmiss
 
C  RELATIVE HUMIDITY
C  -----------------
 
         rho = bmiss
 
C  SCATTER THE PROPER FIRST GUESS/ANALYSIS VALUES
C  ----------------------------------------------
 
         bak_8(1,l)  = pfc
         bak_8(2,l)  = qob
         bak_8(3,l)  = tob
         bak_8(4,l)  = zob
         bak_8(5,l)  = uob
         bak_8(6,l)  = vob
         bak_8(7,l)  = pwo
         bak_8(8,l)  = pw1o
         bak_8(9,l)  = pw2o
         bak_8(10,l) = pw3o
         bak_8(11,l) = pw4o
         bak_8(12,l) = rho
 
   10 ENDDO
      ENDIF
 
      RETURN
      SUBROUTINE gblevn03(SUBSET) ! FORMERLY SUBROUTINE GETFC …
      END
C***********************************************************************
C***********************************************************************
C>  Get observation error
      SUBROUTINE gblevn04 ! FORMERLY SUBROUTINE GETOE
 
      dimension oemin(2:6)
      REAL(8)   OBS_8,QMS_8,BAK_8,SID_8
      REAL(8)   BMISS
 
      COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
     $ xob,yob,dhr,typ,nlev
      COMMON /gbevff/ bmiss
 
      DATA oemin /0.5,0.1,1.0,0.5,1.0/
 
C  CLEAR THE EVENT ARRAY
C  ---------------------
 
      bak_8 = bmiss
 
C  LOOP OVER LEVELS LOOKING UP THE OBSERVATION ERROR
C  -------------------------------------------------
 
      IF(nlev.GT.0)  THEN
      DO l=1,nlev
 
         pob  = obs_8( 1,l)
         qob  = obs_8( 2,l)
         tob  = obs_8( 3,l)
         wob  = max(obs_8(5,l),obs_8(6,l))
         pwo  = obs_8( 7,l)
         pw1o = obs_8( 8,l)
         pw2o = obs_8( 9,l)
         pw3o = obs_8(10,l)
         pw4o = obs_8(11,l)
         cat  = obs_8(12,l)
 
         IF(cat .EQ.0    ) bak_8( 1,l) = oefg01(pob,typ,5,oemin(5))
         IF(qob .LT.bmiss) bak_8( 2,l) = oefg01(pob,typ,3,oemin(3))
         IF(tob .LT.bmiss) bak_8( 3,l) = oefg01(pob,typ,2,oemin(2))
         IF(wob .LT.bmiss) bak_8( 5,l) = oefg01(pob,typ,4,oemin(4))
         IF(pwo .LT.bmiss) bak_8( 6,l) = oefg01(pob,typ,6,oemin(6))
         IF(pw1o.LT.bmiss) bak_8( 7,l) = oefg01(pob,typ,6,oemin(6))
         IF(pw2o.LT.bmiss) bak_8( 8,l) = oefg01(pob,typ,6,oemin(6))
         IF(pw3o.LT.bmiss) bak_8( 9,l) = oefg01(pob,typ,6,oemin(6))
         IF(pw4o.LT.bmiss) bak_8(10,l) = oefg01(pob,typ,6,oemin(6))
 
      ENDDO
      ENDIF
 
      RETURN
      SUBROUTINE gblevn04 ! FORMERLY SUBROUTINE GETOE …
      END
C***********************************************************************
C***********************************************************************
C>  SUBROUTINE GBLEVN06 - 2D LINEAR HORIZONTAL INTERPOLATION
C-----------------------------------------------------------------------
      SUBROUTINE gblevn06(XOB,YOB) ! FORMERLY SUBROUTINE HTERP
 
      USE gblevn_module
 
      COMMON /gbevee/ psi,zsi,ti(500),ui(500),vi(500),qi(500),
     x     zint(500),pint(500),pintlog(500),plev(500),plevlog(500)
 
      DATA rog   / 29.261 /
 
 
C  CALCULATE HORIZONTAL WEIGHTS AND INTERPOLATE
C  --------------------------------------------
 
      wx = xob/dlon + 1.0
      i0 = wx
      i1 = mod(i0,imax) + 1
      wx = wx-i0
 
      wy = (yob+90.)/dlat + 1.0
      j0 = wy
      j1 = min(j0+1,jmax)
      wy = wy-j0
 
C  HTERP FOR SURFACE HEIGHT
C  ------------------------
 
      p1  = iar12z(i0,j0)
      p2  = iar12z(i0,j1)
      p3  = iar12z(i1,j0)
      p4  = iar12z(i1,j1)
      p5  = p1+(p2-p1)*wy
      p6  = p3+(p4-p3)*wy
      zsi = p5+(p6-p5)*wx
 
C  HTERP FOR SURFACE PRESSURE
C  --------------------------
 
      p1  = iar13p(i0,j0)
      p2  = iar13p(i0,j1)
      p3  = iar13p(i1,j0)
      p4  = iar13p(i1,j1)
      p5  = p1+(p2-p1)*wy
      p6  = p3+(p4-p3)*wy
      psi = p5+(p6-p5)*wx
 
C  HTERP FOR UPA T,U,V,Q
C  ---------------------
 
      DO k=1,kmax
 
         p1 = iar14t(i0,j0,k)
         p2 = iar14t(i0,j1,k)
         p3 = iar14t(i1,j0,k)
         p4 = iar14t(i1,j1,k)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         ti(k) = p5+(p6-p5)*wx
 
         p1 = iar15u(i0,j0,k)
         p2 = iar15u(i0,j1,k)
         p3 = iar15u(i1,j0,k)
         p4 = iar15u(i1,j1,k)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         ui(k) = p5+(p6-p5)*wx
 
         p1 = iar16v(i0,j0,k)
         p2 = iar16v(i0,j1,k)
         p3 = iar16v(i1,j0,k)
         p4 = iar16v(i1,j1,k)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         vi(k) = p5+(p6-p5)*wx
 
         p1 = iar17q(i0,j0,k)
         p2 = iar17q(i0,j1,k)
         p3 = iar17q(i1,j0,k)
         p4 = iar17q(i1,j1,k)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         qi(k) = p5+(p6-p5)*wx
 
C  Layer Pressure
C  --------------
 
         p1 = iarpsl(i0,j0,k)
         p2 = iarpsl(i0,j1,k)
         p3 = iarpsl(i1,j0,k)
         p4 = iarpsl(i1,j1,k)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         plev(k) = p5+(p6-p5)*wx
 
C  Interface Pressure
C  ------------------
 
         p1 = iarpsi(i0,j0,k+1)
         p2 = iarpsi(i0,j1,k+1)
         p3 = iarpsi(i1,j0,k+1)
         p4 = iarpsi(i1,j1,k+1)
         p5 = p1+(p2-p1)*wy
         p6 = p3+(p4-p3)*wy
         pint(k+1) = p5+(p6-p5)*wx
 
      ENDDO
 
C  Compute interface heights
C  -------------------------
 
      zint(1) = zsi
      pint(1) = psi
      pintlog(1) = log(pint(1))
      do k=2,kmax
         k0 = k-1
         zint(k) = zint(k0) - rog*ti(k0)*log(pint(k)/pint(k0))
         pintlog(k) = log(pint(k))
      enddo
      pint(kmax+1) = 0.0
 
C  Compute log(pressure) at layer midpoints
C  ----------------------------------------
 
      do k=1,kmax
         plevlog(k) = log(plev(k))
      enddo
 
ccccc print *,' pint=',pint(1:kmax)
ccccc print *,' zint=',zint(1:kmax)
 
      RETURN
      SUBROUTINE gblevn06(XOB,YOB) ! FORMERLY SUBROUTINE HTERP …
      END
 
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C>
C> SUBPROGRAM:    OEFG01
C>   PRGMMR: D. A. KEYSER     ORG: NP22       DATE: 2007-09-14
C>
C> ABSTRACT: FUNCTION WHICH RETURNS THE OBSERVATION ERROR FOR A
C>   REQUESTED VARIABLE INTERPOLATED TO A DEFINED PRESSURE LEVEL FOR A
C>   DEFINED REPORT TYPE.  IT IS OBTAINED FROM AN INPUT ARRAY CONTAINING
C>   OBSERVATION ERRORS ON FIXED PRESSURE LEVELS BY VARIABLE AND REPORT
C>   TYPE (READ EARLIER FROM THE EXTERNAL OBSERVATION ERROR TABLE)
C>
C> PROGRAM HISTORY LOG:
C> 1995-05-17  J. WOOLLEN (NP20) - ORIGINAL AUTHOR (FUNCTION OEF)
C> 2007-09-14 D. A. KEYSER -- MODIFIED TO USE EXACT LOGIC AS IN GSI
C>    (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE TYPE,
C>     LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO MB, A
C>     FEW OTHER MINOR CHANGES)
C>
C> USAGE:    XX = OEFG01(P,TYP,IE,OEMIN)
C>   INPUT ARGUMENT LIST:
C> @param P        - REAL PRESSURE LEVEL (MB) TO INTERPOLATE OBS ERROR TO
C> @param TYP      - REAL PREPBUFR REPORT TYPE
C> @param IE       - VARIABLE TYPE BEING INTERPOLATED (=2 - TEMPERATURE,
C>              - =3 - MOISTURE, =4 - WIND, =5 - SURFACE PRESSURE, =6 -
C>              - PRECIPITABLE WATER)
C>              - (USED ONLY IN PREVENTS MODE)
C> @param OEMIN    - REAL MINIMUM VALUE FOR OBS ERROR (FOR VARIABLE BEING
C>              - INTERPOLATED)
C> @return OEFG01 OBSERVATION ERROR
C>
C> REMARKS: 'OEFG01' RETURNED IS OBSERVATION ERROR FOR VARIABLE "IE" IN
C>     REPORT TYPE "TYP", INTERPOLATED TO PRESSURE LEVEL "P".
C>
C>
C> ATTRIBUTES:
C>   LANGUAGE: FORTRAN 90
C>   MACHINE:  NCEP WCOSS
C>
      FUNCTION oefg01(P,TYP,IE,OEMIN)
 
      REAL(8) bmiss
 
      COMMON /gbevdd/errs(300,33,6)
      COMMON /gbevff/ bmiss
 
      oefg01 = bmiss
 
C  LOOK UP ERRORS FOR PARTICULAR OB TYPES
C  --------------------------------------
 
      kx = typ
 
      p = max(0.,min(p,2000.))
 
      IF(p.GE.errs(kx,1,1)) k1 = 1
 
      DO kl = 1,32
         IF(p.GE.errs(kx,kl+1,1).AND.p.LE.errs(kx,kl,1))  k1 = kl
      ENDDO
 
      IF(p.LE.errs(kx,33,1)) k1 = 5
 
      k2 = k1 + 1
 
      ediff = errs(kx,k2,1) - errs(kx,k1,1)
      IF(abs(ediff).GT.0.0)  THEN
         del = (p - errs(kx,k1,1))/ediff
      ELSE
         del = bmiss
      ENDIF
 
      del = max(0.,min(del,1.0))
 
      oefg01 = ((1.0 - del) * errs(kx,k1,ie)) + (del * errs(kx,k2,ie))
      oefg01 = max(oefg01,oemin)
 
C  SET MISSING ERROR VALUE TO "BMISS"
C  ----------------------------------
 
      IF(oefg01.GE.5e5) oefg01 = bmiss
 
      RETURN
 
      FUNCTION oefg01(P,TYP,IE,OEMIN) …
      END
 
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C>
C> SUBPROGRAM:    GBLEVN08
C>   PRGMMR: S. MELCHIOR      ORG: NP22       DATE: 2014-03-25
C>
C> ABSTRACT: CREATE VIRTUAL TEMPERATURE EVENTS WITHIN GBLEVENTS
C>   SUBROUTINE.  FOR ALL TYPES EXCEPT RASS, THIS CONSISTS OF FIRST RE-
C>   CALCULATING THE SPECIFIC HUMIDITY FROM THE REPORTED DEWPOINT
C>   TEMPERATURE AND PRESSURE, FOLLOWED BY THE CALCULATION OF VIRTUAL
C>   TEMPERATURE FROM THE JUST-CALCULATED SPECIFIC HUMIDITY AND THE
C>   REPORTED (SENSIBLE) TEMPERATURE.  THE RE-CALCULATED SPECIFIC
C>   HUMIDITY IS THEN ENCODED AS A STACKED EVENT TO BE LATER WRITTEN
C>   INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP", REASON CODE 0).
C>   IF THE NAMELIST SWITCH DOVTMP IS TRUE, THEN THE JUST-CALCULATED
C>   VIRTUAL TEMPERATURE IS THEN ALSO ENCODED AS A STACKED EVENT TO BE
C>   LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP",
C>   REASON CODE 0, 2 OR 6).  FOR RASS DATA, SPECIFIC HUMIDITY IS
C>   MISSING HOWEVER IF THE NAMELIST SWITCH DOVTMP IS TRUE, A SIMPLE
C>   COPY OF THE REPORTED (VIRTUAL) TEMPERATURE IS ENCODED AS A STACKED
C>   EVENT TO BE LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM
C>   "VIRTMP", REASON CODE 3).  FOR SURFACE DATA WITH A MISSING PMSL, IF
C>   DOVTMP=T AND DOPMSL=T AND A VIRTUAL TEMPERATURE HAS BEEN COMPUTED,
C>   CALCULATE AN ESTIMATED PMSL AND ENCODE IT INTO PREPBUFR FILE ALONG
C>   WITH AN INDICATOR THAT IS WAS DERIVED HERE. THIS SUBROUTINE IS
C>   CURRENTLY ONLY CALLED FOR SURFACE LAND ("ADPSFC"), MARINE
C>   ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA") OR SATELLITE
C>   TEMPERATURE RETRIEVAL ("SATEMP") DATA TYPES WHEN SWITCH
C>   "ADPUPA_VIRT" IS FALSE AND ONLY FOR SURFACE LAND ("ADPSFC"), MARINE
C>   ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA"), SATELLITE
C>   TEMPERATURE RETRIEVAL ("SATEMP") OR RAOB/DROP/MULTI-LVL RECCO
C>   ("ADPUPA") DATA TYPES WHEN SWITCH "ADPUPA_VIRT" IS TRUE.  IT IS
C>   ALSO ONLY CALLED IN THE PREVENTS MODE.  THIS ROUTINE IS CALLED ONCE
C>   FOR EACH VALID REPORT IN THE PREPBUFR FILE.
C>
C> PROGRAM HISTORY LOG:
C> 1995-05-17  J. WOOLLEN (NP20) - ORIGINAL AUTHOR
C> 1997-06-01  D.A. KEYSER - STREAMLINED, ADDED SWITCH DOVTMP
C> 1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE
C>     NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT
C>     TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD
C>     QUAL. MARKER (8)
C> 2004-08-30 D. A. KEYSER -- FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY
C>     OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF
C>     DOVTMP IS TRUE, GETS NEW REASON CODE 3
C> 2006-07-14 D. A. KEYSER --  PROCESSES REPORTS IN MESSAGE TYPE ADPUPA
C>     (I.E., RAOBS, DROPS, MULTI-LEVEL RECCOS) WITH SAME LOGIC AS IN
C>     SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR WHEN NEW NAMELIST
C>     SWITCH "ADPUPA_VIRT" IS TRUE {NORMALLY "ADPUPA_VIRT" IS FALSE
C>     (DEFAULT) BECAUSE SUBSEQUENT PROGRAM PREPOBS_CQCBUFR PERFORMS
C>     THIS FUNCTION}
C> 2007-09-14 D. A. KEYSER -- FOR NON-"ADPUPA" TYPES, Q.M. 9 IS NOW
C>     ASSIGNED TO CALCULATED VIRT. TEMPS IF THE MOISTURE Q.M. IS 9 OR
C>     15 AND ORIG. TEMP NOT "BAD", THESE "VIRTMP" EVENTS RECEIVE NEW
C>     REASON CODE 4, HAD RECEIVED Q.M. 8 WITH REASON CODE 2 LIKE VIRT.
C>     TEMPS CALCULATED FROM "BAD" MOISTURE - THIS MEANS ONLY TRULY
C>     "BAD" VIRT. TEMPS WILL NOW GET Q.M. 8 AND VIRT. TEMPS FLAGGED FOR
C>     NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M. 9
C>     (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES NOT
C>     EVEN CONSIDER OBS WITH Q.M. 8); FOR "ADPUPA" TYPES, Q.M. 3 IS NOW
C>     ASSIGNED TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS
C>     TRULY BAD (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM
C>     IS 1 OR 2 AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN
C>     QQM WAS 9 OR 15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES,
C>     ENCODES A SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A
C>     "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO
C>     WHAT IS ALREADY DONE FOR "RASSDA" TYPES)
C> 2013-04-12  D. A. KEYSER -- DON'T ALLOW CALCULATED Q TO BE < 0 WHICH
C>     CAN OCCUR ON WCOSS FOR CASES OF HORRIBLY BAD MESONET DATA
C> 2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH,
C>     WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION
C>     PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST-
C>     COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN
C>     REPORTED PMSL IS MISSING.  DOVTMP MUST BE TRUE AND DOANLS MUST BE
C>     FALSE ("PREVENTS" MODE). THE DERIVED PMSL EITHER GETS A QUALITY
C>     MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE QUALITY MARK
C>     ("PQM"), WHICHEVER IS GREATER. THE VALUE OF THE PMSL INDICATOR
C>     (NEW MNEMONIC "PMIN") IS SET TO 1 TO DENOTE PMSL WAS DERIVED
C>     RATHER THAN OBSERVED.  THE DEFAULT FOR "DOPMSL" IS FALSE
C>     (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS).  IT IS FORCED TO BE
C>     FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE).
C>
C> USAGE:    CALL GBLEVN08(IUNITP)
C>   INPUT ARGUMENT LIST:
C> @param IUNITP   - BUFR OUTPUT FILE UNIT
C> @param SUBSET   - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
C>              - REPORT BEING PROCESSED
C>
C> REMARKS: WILL IMMEDIATELY RETURN TO CALLING PROGRAM IF ANY OF THE
C>   FOLLOWING CONDITIONS EXIST: THERE ARE NO LEVELS OF VALID DEWPOINT,
C>   OBS, TEMPERATURE Q.M. OR SPEC. HUMIDITY Q.M. IN THE INPUT PREPBUFR
C>   FILE FOR THE REPORT.  WILL NOT ATTEMPT EITHER SPEC. HUMIDITY NOR
C>   VIRT. TEMP CALC. ON A GIVEN LEVEL IF ANY OF THE FOLLOWING
C>   CONDITIONS EXIST: REPORTED PRESSURE OBS IS MISSING, REPORTED
C>   (SENSIBLE) TEMPERATURE OBS IS MISSING, OR REPORTED DEWPOINT OBS IS
C>   MISSING.
C> 
C> ATTRIBUTES:
C>   LANGUAGE: FORTRAN 90
C>   MACHINE:  NCEP WCOSS
C>
      SUBROUTINE gblevn08(IUNITP,SUBSET) ! FORMERLY SUBROUTINE VTPEVN
 
      parameter(rd=287.04, g=9.81)
      CHARACTER*80 EVNSTQ,EVNSTV,evnstp
      CHARACTER*8  SUBSET,stnid
      REAL(8)      TDP_8(255),TQM_8(255),QQM_8(255),BAKQ_8(4,255),
     $             bakv_8(4,255),bakp_8(3),obs_8,qms_8,bak_8,
     $             sid_8,pqm_8
      real(8)      pmo_8,zob_8,pmsl_8
      REAL(8)      BMISS
 
      LOGICAL      EVNQ,EVNV,DOVTMP,TROP,ADPUPA_VIRT,DOBERR,DOFCST,
     $             some_fcst,fcst,virt,satmqc,recalc_q,doprev,
     $             evnp,dopmsl,surf
 
      COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
     $ xob,yob,dhr,typ,nlev
      COMMON /gbevbb/ pvcd,vtcd
      COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
     $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
      COMMON /gbevff/ bmiss
 
      DATA evnstq /'QOB QQM QPC QRC'/
      DATA evnstv /'TOB TQM TPC TRC'/
      data evnstp /'PMO PMQ PMIN'/
 
      equivalence(sid_8,stnid)
 
C-----------------------------------------------------------------------
C FCNS BELOW CONVERT TEMP/TD (K) & PRESS (MB) INTO SAT./ SPEC. HUM.(G/G)
C-----------------------------------------------------------------------
      es(t) = 6.1078*exp((17.269*(t - 273.16))/((t - 273.16)+237.3))
      qs(t,p) = (0.622*es(t))/(p-(0.378*es(t)))
C-----------------------------------------------------------------------
C FCN BELOW derives mean sea-level pressure (mb)from P (station
C  pressure, mb), Tv (virtual temperature, K), and Z (station height/
C  elevation, m) for surface reports
C-----------------------------------------------------------------------
      pmsl_fcn(p,tv,z) = p*exp((g*z)/(rd*tv))
C-----------------------------------------------------------------------
 
C  CLEAR TEMPERATURE AND SPECIFIC HUMIDITY EVENTS
C  ----------------------------------------------
 
      evnq   = .false.
      evnv   = .false.
      evnp   = .false.
      bakq_8 = bmiss
      bakv_8 = bmiss
      bakp_8 = bmiss
      trop   = .false.
 
      surf = (subset.eq.'ADPSFC'.or.subset.eq.'SFCSHP'.or.
     $           subset.eq.'MSONET')
 
C  GET DEWPOINT TEMPERATURE AND CURRENT T,Q QUALITY MARKERS
C  Also get mean sea level pressure, station pressure quality mark,
C   and station height (elevation) for surface reports if dopmsl=T
C  ----------------------------------------------------------------
 
      CALL ufbint(-iunitp,tdp_8,1,255,nltd,'TDO')
      CALL ufbint(-iunitp,qqm_8,1,255,nlqq,'QQM')
      CALL ufbint(-iunitp,tqm_8,1,255,nltq,'TQM')
      if(surf.and.dopmsl) then
         call ufbint(-iunitp,pmo_8,1,1,nlpm,'PMO')
         call ufbint(-iunitp,zob_8,1,1,nlzo,'ZOB')
         call ufbint(-iunitp,pqm_8,1,1,nlpq,'PQM')
      endif
      IF(subset.NE.'RASSDA  '.AND.subset.NE.'SATEMP  ') THEN
         IF(nltd.EQ.0) RETURN
         IF(nlqq.EQ.0) RETURN
      ENDIF
      IF(nltq.EQ.0) RETURN
      IF(subset.NE.'RASSDA  '.AND.subset.NE.'SATEMP  ') THEN
         IF(nltd.NE.nlev) THEN
            print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLTD  NLEV - STOP 61")'
            CALL errexit(61)
         ENDIF
         IF(nlqq.NE.nlev) THEN
            print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLQQ  NLEV - STOP 63")'
            CALL errexit(63)
         ENDIF
      ENDIF
      IF(nltq.NE.nlev) THEN
         print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLTQ  NLEV - STOP 62")'
         CALL errexit(62)
      ENDIF
 
C  COMPUTE SPECIFIC HUMIDITY AND VIRTUAL TEMPERATURE USING REPORTED DEWP
C  For surface reports, also calculate PMSL if it is missing, dopmsl=T,
C   and a virtual temperature has been computed
C  ---------------------------------------------------------------------
 
      IF(nlev.GT.0)  THEN
      DO l=1,nlev
         pob = obs_8(1,l)
         tdo = tdp_8(l)
         tob = obs_8(3,l)
         cat = obs_8(12,l)
         IF(dovtmp) THEN
            IF(subset.EQ.'RASSDA  '.OR.subset.EQ.'SATEMP  ') THEN
               IF(tob.LT.bmiss) THEN
                  bakv_8(1,l) = tob
                  bakv_8(2,l) = tqm_8(l)
                  bakv_8(3,l) = vtcd
                  bakv_8(4,l) = 3
                  evnv = .true.
                  cycle
               ENDIF
            ENDIF
         ENDIF
         IF(pob.LT.bmiss .AND. tob.LT.bmiss
     $                   .AND. tdo.LT.bmiss) THEN
            IF(qqm_8(l).GT.3) THEN
C  Don't update q or calculate Tv if bad moisture obs fails sanity check
C   (also don't calc. PMSL if it is missing and dopmsl=T for sfc rpts)
cdak           IF(TDO.LT.-103.15 .OR. TDO.GT.46.83 .OR. POB.LT.0.1 .OR.
cdak $          POB.GT.1100.)
cdak $ print *, '&&& bad QM fails sanity check'
               IF(tdo.LT.-103.15 .OR. tdo.GT.46.83 .OR. pob.LT.0.1 .OR.
     $          pob.GT.1100.)  cycle
            ENDIF
            qob = qs(tdo+273.16,pob)
ccccc       BAKQ_8(1,L) = QOB*1E6 ! dak fix 2/27/13: can't be > bmiss
                                  !  else flting pt overflow in BUFRLIB
ccc         IF(QOB*1E6.LT.BMISS) BAKQ_8(1,L) = QOB*1E6
                                  ! dak add'l fix 4/12/13: don't allow
                                  !  calc. q to be < 0 which can occur
                                  !  on WCOSS for cases of horribly bad
                                  !  mesonet data
            IF(qob*1e6.LT.bmiss .AND. qob.GT.0.) bakq_8(1,l) = qob*1e6
            bakq_8(2,l) = qqm_8(l)  ! Moist qm same as before for
                                    !  re-calc. q
            bakq_8(3,l) = vtcd
            bakq_8(4,l) = 0     ! Re-calc. q gets unique reason code 0
            evnq = .true.
C  If message type ADPUPA, test this level to see if at or above trop
C   (trop must be above 500 mb to pass test; if no trop level found
C   assume it's at 80 mb)
C  Don't calculate Tv on this level if at or above trop (doesn't affect
C   q calculation)
            trop = (subset.EQ.'ADPUPA  ' .AND.
     $       ((cat.EQ.5 .AND. pob.LT.500.) .OR. pob.LT. 80. .OR. trop))
            IF(dovtmp .AND. .NOT.trop) THEN
               bakv_8(1,l) = (tob+273.16)*(1.+.61*qob)-273.16
               bakv_8(3,l) = vtcd
               IF(subset.EQ.'ADPUPA  ') THEN
C  Message type ADPUPA comes here
                  IF((qqm_8(l).LT.4.OR.qqm_8(l).EQ.9.OR.qqm_8(l).EQ.15)
     $             .OR. tqm_8(l).EQ.0 .OR. tqm_8(l).GT.3
     $             .OR. pob.LE.700.) THEN
                     bakv_8(2,l) = tqm_8(l) ! Tv qm same as for T when
                                            !  q ok or q flagged by
                                            !  PREPRO (but not bad) 
                     bakv_8(4,l) = 0    ! Tv gets unique reason code 0
                  ELSE
                     bakv_8(2,l) = 3    !Tv qm susp for bad moist below
                                        ! 700 mb
                     bakv_8(4,l) = 6    !Tv gets unique reason code 6
                  ENDIF
               ELSE
C  All other message types come here
                  IF(qqm_8(l).LT.4) THEN
                     bakv_8(2,l) = tqm_8(l) ! Tv qm same as for T when
                                            !  q ok
                     bakv_8(4,l) = 0    ! Tv gets unique reason code 0
                  ELSE IF((qqm_8(l).EQ.9.OR.qqm_8(l).EQ.15).AND.
     $                    (tqm_8(l).LE.3.OR.tqm_8(l).GE.15.OR.
     $                    tqm_8(l).EQ.9)) THEN
cdak  print'(" %%% process tvirt on lvl ",I0," for missing moist obs ",
cdak $ "error/high-up moist case when orig. T not ""bad"" (set TQM_8=",
cdak $ "9)")', l
                     bakv_8(2,l) = 9 ! Tv qm 9 for moist w/ missing obs
                                     !  error or moist flagged by
                                     !  PREPRO (but not bad) and T qm
                                     !  orig not "bad"
                     bakv_8(4,l) = 4 ! Tv gets unique reason code 4
                  ELSE
cdak  print'(" %%% process tvirt on lvl ",I0," for ""bad"" QQM_8 case ",
cdak $ "or missing moist obs error/high-up moist w/ ""bad"" TQM_8 case",
cdak $ " (set TQM_8=8)")', l
                     bakv_8(2,l) = 8 ! Tv qm 8 (bad) for "bad" moist or
                                     !  moist w/ missing obs error or
                                     !  moist flagged by PREPRO (but
                                     !  not bad) and T qm orig "bad"
                     bakv_8(4,l) = 2 ! Tv gets unique reason code 2
                  ENDIF
                  if(surf) then
c  For sfc rpts & dopmsl=T, derive Pmsl in cases when it is not reported
c  DAK: Note right now this can only happen for sfc rpts where a Tv is
c       set to be calculated and where it is successfully calculated.
c       Eventually this may be able to be relaxed such that PMSL can be
c       derived even if, e.g., no moisture were reported (in this case
c       Ts would have to be used, still a decent estimate of PMSL could
c       liekly be made). Might also be able to be derived if no T rpted.
                    if(dopmsl) then
                      if(ibfms(pmo_8).ne.0) then
                        tv = bakv_8(1,1) + 273.16
                        zob=zob_8
                        pmsl_8=pmsl_fcn(pob,tv,zob)
                        bakp_8(1) = pmsl_8
                        bakp_8(2) = max(3,int(pqm_8))! qm>=3 b/c derived
                        bakp_8(3) = 1. ! pmin=1 for derived value
cccccccccc              write(*,4000) stnid,bakp_8(3),bakp_8(2)
 4000 format('--> ID ',a8,' Pmsl missing - derived from Pstn; ',
     $ 'PMIN = ',f4.1,' PMQ = ',f4.1,'')
                        evnp = .true.
c Diagnostics for Pmsl values that are suspiciously high
                        if(pmsl_8.gt.1060) then
                          write(*,4001) stnid,pob,bakp_8(1)
 4001 format('--> ID ',a8,' Derived PMSL unrealistic; FLAG; ',
     $ 'POB = ',f7.2,' PMO = ',f7.2,'')
                        end if
                      end if
                    end if
                  end if
               ENDIF
               evnv = .true.
            ENDIF
         ENDIF
      ENDDO
      ENDIF
 
C  ENCODE EVENTS INTO REPORT
C  -------------------------
 
      IF(nlev.GT.0)  THEN
         IF(evnq) CALL ufbint(iunitp,bakq_8,4,nlev,iret,evnstq)
         IF(evnv) CALL ufbint(iunitp,bakv_8,4,nlev,iret,evnstv)
         if(nlev.eq.1.and.evnp)
     $    call ufbint(iunitp,bakp_8,3,nlev,iret,evnstp)
      ENDIF
 
      RETURN
      SUBROUTINE gblevn08(IUNITP,SUBSET) ! FORMERLY SUBROUTINE VTPEVN …
      END
C#######################
C#######################
C#######################
C#######################
C#######################
C#######################
C***********************************************************************
C***********************************************************************
C> Do something
C> @param IUNITF   - 2-WORD ARRAY:
C>        For SIGIO input:
C>              - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
C>                IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
C>                AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
C>                IS A MULTIPLE OF 3
C>              - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
C>                OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
C>                THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
C>        For NEMSIO input:
C>              - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
C>                (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
C>                THE DATE IN IDATEP
C>              - WORD 2 - NOT USED, SHOULD BE A NULL FILE
C>                
C> @param IDATEP   - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
C> @param im
C> @param jm
C> @param IDRT INTEGER GRID IDENTIFIER
C>             (IDRT=4 FOR GAUSSIAN GRID,
C>             IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C>             IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
      SUBROUTINE gblevn10(IUNITF,IDATEP,IM,JM,IDRT) ! FORMERLY
                                                    ! SUBROUTINE GESRES
      USE gblevn_module
      USE sigio_module
      USE sigio_r_module
 
      IMPLICIT NONE
      INTEGER IUNITF(2), IDATEP, IM, JM, IDRT
      REAL, PARAMETER :: PI180=.0174532
      INTEGER*4, PARAMETER :: ONE=1, ten=10
 
      TYPE(sigio_head) :: HEAD(2)
      TYPE(sigio_dats) :: DATS
      TYPE(sigio_datm) :: DATM
 
      INTEGER*4 IRET,IRET1,IRETS,IMJM4,KM4,IDVM,NTRAC,IUNIT4(2)
 
      INTEGER KFILES,IFILE,JFILE,IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,
     $ mdima,mdimb,mdimc,iromb,maxwv,idir,ns,i,j,k,l,ii,jj,ib,ie
 
      INTEGER IDATE(8,2),JDATE(8,2),KDATE(8,2),KINDX(2)
 
      CHARACTER*6  COORD(3)
      CHARACTER*20 CFILE
 
      REAL FHOUR,RINC(5)
 
      DATA coord /'SIGMA ','HYBRID','GENHYB'/
 
      REAL,         ALLOCATABLE :: cofs(:,:),     cofv(:,:,:)
      REAL,         ALLOCATABLE :: cofs_f(:,:,:), cofv_f(:,:,:,:)
      REAL (KIND=4),allocatable :: grds(:,:,:),   grdv(:,:,:,:),
     $                             wrk1(:,:),     wrk2(:,:)
 
      imax  = im
      jmax  = jm
      imjm4 = im*jm
      iunit4(:) = iunitf(:)
 
      IF(mod(mod(idatep,100),3).EQ.0)  THEN
         kfiles = 1
         kindx = 0
         print 331, mod(idatep,100)
  331    FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',i2.2,
     $    ') IS A MULTIPLE OF 3 - ONLY ONE SIGIO (SIGMA OR HYBRID) ',
     $    'INPUT GLOBAL FILE IS READ,'/16x,'NO TIME INTERPOLATION OF ',
     $    'SPECTRAL COEFFICIENTS IS PERFORMED'/)
      ELSE
         kfiles = 2
         kindx(1) = mod(mod(idatep,100),3)
         kindx(2) = kindx(1) - 3
         print 332, mod(idatep,100)
  332    FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',i2.2,
     $    ') IS NOT A MULTIPLE OF 3 - TWO SPANNING SIGIO (SIGMA OR ',
     $    'HYBRID) INPUT GLOBAL FILES'/16x,'ARE READ AND THE SPECTRAL ',
     $    'COEFFICIENTS ARE INTERPOLATED TO THE PREPBUFR CENTER TIME'/)
      ENDIF
 
C  GET VALID-TIME DATE OF SIGMA OR HYBRID FILE(S), ALSO READ HEADERS
C  -----------------------------------------------------------------
 
      jfile = 0
      rinc  = 0
      DO ifile=1,kfiles
         jfile = ifile
 
         WRITE(cfile,'("fort.",I2.2)') iunitf(ifile)
 
         CALL sigio_rropen(iunit4(ifile),cfile,iret)
         CALL sigio_rrhead(iunit4(ifile),head(ifile),iret1)
         print *,' cfile_sig=',cfile,'sigio_rropen: iret=',iret,
     $   'sigio_rrhead: iret1=',iret1
 
         IF(iret.NE.0)  GO TO 903
         IF(iret1.NE.0) GO TO 904
 
         idate(:,ifile)   = 0
 
         idate(1,ifile)   = head(ifile)%IDATE(4)
         idate(2:3,ifile) = head(ifile)%IDATE(2:3)
         idate(5,ifile)   = head(ifile)%IDATE(1)
 
         fhour = head(ifile)%FHOUR
         print'("  idate=",I5,7I3.2,"  fhour=",F7.3)', idate(:,ifile),
     $    head(ifile)%fhour
 
         IF(idate(1,ifile).LT.100)  THEN
 
C IF 2-DIGIT YEAR FOUND IN GLOBAL SIMGA FILE INITIAL DATE
C  (IDATE(1,IFILE)), MUST USE "WINDOWING" TECHNIQUE TO CREATE A 4-DIGIT
C  YEAR (NOTE: THE T170 IMPLEMENTATION IN JUNE 1998 WAS TO INCLUDE THE
C  WRITING OF A 4-DIGIT YEAR HERE.  PRIOR TO THIS, THE YEAR HERE WAS
C  2-DIGIT.)
 
            print'(" ##GBLEVENTS/GBLEVN10 - 2-DIGIT YEAR FOUND IN ",
     $       "SIGIO (SIGMA OR HYBRID) INPUT GLOBAL FILE ",I0,
     $       "; INITIAL DATE (YEAR IS: ",I0,")")', ifile,idate(1,ifile)
            print'("      - USE WINDOWING TECHNIQUE TO OBTAIN 4-DIGIT",
     $       " YEAR")'
            IF(idate(1,ifile).GT.20)  THEN
               idate(1,ifile) = 1900 + idate(1,ifile)
            ELSE
               idate(1,ifile) = 2000 + idate(1,ifile)
            ENDIF
            print'(" ##GBLEVENTS/GBLEVN10 - CORRECTED 4-DIGIT YEAR IS",
     $       " NOW: ",I0)', idate(1,ifile)
         ENDIF
 
         rinc(2) = fhour
         CALL w3movdat(rinc,idate(:,ifile),jdate(:,ifile))
 
         print 1, ifile,head(ifile)%FHOUR,
     $    (idate(ii,ifile),ii=1,3),idate(5,ifile),(jdate(ii,ifile),
     $    ii=1,3),jdate(5,ifile)
    1    FORMAT(' --> GBLEVENTS: SIGIO (SIGMA OR HYBRID) INPUT GLOBAL ',
     $    'FILE',i2,' HERE IS A ',f5.1,' HOUR FORECAST FROM ',i5.4,
     $    3i3.2,' VALID AT ',i5.4,3i3.2)
 
         kdate(:,ifile) = jdate(:,ifile)
 
         IF(kfiles.EQ.2) THEN
            rinc(2) = real(kindx(ifile))
            CALL w3movdat(rinc,jdate(:,ifile),kdate(:,ifile))
         ENDIF
 
         idatgs_cor = (kdate(1,ifile) * 1000000) + (kdate(2,ifile) *
     $    10000) + (kdate(3,ifile) * 100) + kdate(5,ifile)
 
C  VALID DATES MUST MATCH
C  ----------------------
 
         IF(idatep.NE.idatgs_cor)  GO TO 901
 
      ENDDO
 
 
C  EXTRACT HEADER INFO
C  -------------------
 
      jcap    = head(1)%JCAP
      kmax    = head(1)%LEVS
      km4     = kmax
      idvc    = head(1)%IDVC
      idsl    = head(1)%IDSL
      idvm    = head(1)%IDVM
      ntrac   = head(1)%NTRAC
      nvcoord = head(1)%NVCOORD
      ALLOCATE (vcoord(kmax+1,nvcoord))
      vcoord  = 0.0
      vcoord  = head(1)%VCOORD
 
      sfcpress_id  = mod(head(1)%IDVM,ten)
      thermodyn_id = mod(head(1)%IDVM/ten,ten)
      IF(idvc == 3 .AND. thermodyn_id == 3) THEN
         kmaxs = (ntrac+1)*kmax + 2
      ELSE
         kmaxs = 2*kmax + 2
         ntrac = 1
      ENDIF
 
      ALLOCATE (iar12z(im,jm), iar13p(im,jm))
      ALLOCATE (iar14t(im,jm,kmax),  iar15u(im,jm,kmax),
     $          iar16v(im,jm,kmax),  iar17q(im,jm,kmax),
     $          iarpsl(im,jm,kmax),  iarpsi(im,jm,kmax+1))
 
 
      if(idvc.eq.0)  idvc = 1  ! Reset IDVC=0 to 1 (sigma coord.)
      IF(idvc < 0 .or. idvc > 3) THEN
         print *, '##GBLEVENTS/GBLEVN10: INVALID VERT COORD ID (=',idvc
      ENDIF
 
 
C  DEFINE THE OTHER RESOLUTION PARAMETERS
C  --------------------------------------
 
      jcap1   = jcap+1
      jcap2   = jcap+2
      jcap1x2 = jcap1*2
      mdima   = jcap1*jcap2
      mdimb   = mdima/2+jcap1
      mdimc   = mdimb*2
      imax    = 384
      jmax    = imax/2+1
 
      dlat  = 180./(jmax-1)
      dlon  = 360./imax
 
      print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,coord(idvc)
    2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
     $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
     $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
     $ 'COORD: ',a)
 
      GO TO 902
 
  901 CONTINUE
      print 9901, jfile,(jdate(ii,jfile),ii=1,3),jdate(5,jfile),idatep
 9901 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,' DATE',
     $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH -OR SPAN- PREPBUFR FILE ',
     $ 'CENTER DATE (',i10,') -STOP 68'/)
      CALL errexit(68)
  903 CONTINUE
      print 9903, jfile,iret
 9903 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,
     $ ' RETURNED FROM SIGIO_RROPEN WITH R.C.',i3,' -STOP 70'/)
      CALL errexit(70)
  904 CONTINUE
      print 9904, jfile,iret1
 9904 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,
     $ ' RETURNED FROM SIGIO_RRHEAD WITH R.C.',i3,' -STOP 71'/)
      CALL errexit(71)
  902 CONTINUE
      IF(kmax.GT.500) then
         print'(" ##GBLEVENTS/GBLEVN10 - KMAX TOO BIG = ",I0,
     $    " - UNABLE TO TRANSFORM GLOBAL SIGMA FILE(S) - STOP 69")',
     $    kmax
         CALL errexit(69)
      ENDIF
 
C***********************************************************************
C***********************************************************************
 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
C USAGE:    CALL SPTEZM(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,WAVE,GRID,IDIR)
C   INPUT ARGUMENTS:
C     IROMB    - INTEGER SPECTRAL DOMAIN SHAPE
C                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C     MAXWV    - INTEGER SPECTRAL TRUNCATION
C     IDRT     - INTEGER GRID IDENTIFIER
C                (IDRT=4 FOR GAUSSIAN GRID,
C                 IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C                 IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
C     IMAX     - INTEGER EVEN NUMBER OF LONGITUDES
C     JMAX     - INTEGER NUMBER OF LATITUDES
C     KMAX     - INTEGER NUMBER OF FIELDS TO TRANSFORM
C     WAVE     - REAL (2*MX,KMAX) WAVE FIELD IF IDIR>0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     GRID     - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR<0
C     IDIR     - INTEGER TRANSFORM FLAG
C                (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
C   OUTPUT ARGUMENTS:
C     WAVE     - REAL (2*MX,KMAX) WAVE FIELD IF IDIR<0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     GRID     - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR>0
 
 
C USAGE:    CALL SPTEZMV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
C    &                   WAVED,WAVEZ,GRIDU,GRIDV,IDIR)
C   INPUT ARGUMENTS:
C     WAVED    - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR>0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     WAVEZ    - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     GRIDU    - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR<0
C     GRIDV    - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR<0
C   OUTPUT ARGUMENTS:
C     WAVED    - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR<0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     WAVEZ    - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0
C                WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
C     GRIDU    - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR>0
C     GRIDV    - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR>0
 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
      allocate (cofs_f(mdima,kmaxs,2), cofv_f(mdima,kmax,2,2))
 
      iromb = 0
      maxwv = jcap
      if (idrt < 0 .or. idrt > 256) idrt  = 0
      idir  = 1
 
      IF(kindx(1).EQ.0)  THEN
         kfiles = 1
      ELSE
         kfiles = 2
      ENDIF
 
C  Allocate for sigio read
C  -----------------------
 
      sfcpress_id  = mod(head(1)%IDVM,ten)
      thermodyn_id = mod(head(1)%IDVM/ten,ten)
 
      print *,'in sig sfcpress_id=',sfcpress_id,' thermodyn_id=',
     $ thermodyn_id,' ntrac=',ntrac
      print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,
     $     ' nvcoord=',nvcoord
 
      DO ifile=1,kfiles
         CALL sigio_aldats(head(ifile),dats,irets)
         CALL sigio_aldatm(head(ifile),one,km4,datm,irets)
                                          ! Read surface fields
         CALL sigio_rrdats(iunit4(ifile),head(ifile),dats,irets)
 
         IF(irets.NE.0) THEN
            print *,' irets from sigio_rrdats = ', irets
            RETURN
         ENDIF
 
         DO i=1,mdima
            cofs_f(i,1,ifile) = dats%HS(i)
            cofs_f(i,2,ifile) = dats%PS(i)
         ENDDO
 
                                 ! Read fields on levels 1 through kmax
         CALL sigio_rrdatm(iunit4(ifile),head(ifile),datm,irets)
         IF(irets.NE.0) THEN
            print *,' irets from sigio_rrdatm = ', irets
            RETURN
         ENDIF
ccccc    print *,' aft sigio_rrdatm irets=',irets
 
         ie = kmax + 2
         cofs_f(:,3:ie,ifile)        = datm%T
         DO i=1,ntrac
           ib = ie + 1
           ie = ib + kmax - 1
           cofs_f(:,ib:ie,ifile) = datm%Q(:,1:kmax,i)
         ENDDO
         cofv_f(:,:,1,ifile)             = datm%D
         cofv_f(:,:,2,ifile)             = datm%Z
 
         CALL sigio_axdats(dats,irets)
         CALL sigio_axdatm(datm,irets)
      ENDDO
 
ccccc print *,' after sigio_axdatm'
 
      ALLOCATE (cofs(mdima,kmaxs),     cofv(mdima,kmax,2))
      ALLOCATE (grds(imax,jmax,kmaxs), grdv(imax,jmax,kmax,2))
      ALLOCATE (wrk1(imax*jmax,kmax),  wrk2(imax*jmax,kmax+1))
 
      IF(kfiles.EQ.1)  THEN
         DO i = 1,mdima
            cofs(i,1:kmaxs)    = cofs_f(i,1:kmaxs,1)
            cofv(i,1:kmax,1:2) = cofv_f(i,1:kmax,1:2,1)
         ENDDO
 
      ELSE
         cofs=
     $   ((abs(kindx(2))*cofs_f(:,:,1))  +(kindx(1)*cofs_f(:,:,2)))/3.
         cofv=
     $   ((abs(kindx(2))*cofv_f(:,:,:,1))+(kindx(1)*cofv_f(:,:,:,2)))/3.
      ENDIF
 
      DEALLOCATE (cofs_f, cofv_f)
 
      CALL sptezm(iromb,maxwv,idrt,imax,jmax,kmaxs,cofs,grds,idir)
      CALL sptezmv(iromb,maxwv,idrt,imax,jmax,kmax,
     &         cofv(1,1,1),cofv(1,1,2),grdv(1,1,1,1),grdv(1,1,1,2),idir)
 
      IF( sfcpress_id == 2 ) THEN   ! for enthalpy version
         grds(:,:,2) = 1000.0*grds(:,:,2)         ! Now in Pa
      ELSE
         grds(:,:,2) = 1000.0*exp(grds(:,:,2))    ! Now in Pa
      ENDIF
 
      DO ns=1, kmaxs
         CALL gblevn11(imax,jmax,grds(1,1,ns))
      ENDDO
      DO j=1,jmax
         DO i=1,imax
            iar12z(i,j) = grds(i,j,1)         ! Orography
            iar13p(i,j) = grds(i,j,2) * 0.01  ! Surface Pressure in hPa
         ENDDO
      ENDDO
 
      IF(thermodyn_id == 3 .AND. idvc == 3) THEN
         grds(:,:,3:2+kmax) = grds(:,:,3:2+kmax) / head(1)%CPI(1)
         print *,' cpi(0)=',head(1)%cpi(1)
      ENDIF
      CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
     $                 grds(1,1,2),grds(1,1,3),pm=wrk1,pd=wrk2(1,2))
 
      DO j=1,jmax
         jj = (j-1)*imax
         DO i=1,imax
            wrk2(i+jj,1) = grds(i,j,2)                     ! in Pa
         ENDDO
      ENDDO
      DO l=1,kmax
         wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1)             ! in Pa
      ENDDO
 
ccccc print *,' wrk1=',wrk1(1001,:)
ccccc print *,' wrk2=',wrk2(1001,:)
 
ccccc CALL GBLEVN11(IMAX,JMAX,WRK2(1,KMAX+1))
ccccc DO L=1,KMAX
ccccc    CALL GBLEVN11(IMAX,JMAX,WRK1(1,L))
ccccc    CALL GBLEVN11(IMAX,JMAX,WRK2(1,L))
ccccc ENDDO
 
      IF(thermodyn_id == 3 .AND. idvc == 3) THEN
                                 ! Convert from enthalpy to temperature
         grds(:,:,3:2+kmax) = grds(:,:,3:2+kmax) * head(1)%CPI(1)
         CALL sigio_cnvtdv(imjm4,imjm4,km4,idvc,idvm,ntrac,iret,
     $                     grds(1,1,3),grds(1,1,3+kmax),head(1)%CPI,1_4)
                                 ! Convert back to virtual temperature
         grds(:,:,3:kmax+2) = grds(:,:,3:kmax+2) *
     $            (1.+(461.50/287.05-1)*grds(:,:,3+kmax:2+kmax*2))
      ENDIF
 
      DO l=1,kmax
         DO k=1,2
            CALL gblevn11(imax,jmax,grdv(1,1,l,k))
         ENDDO
         DO j=1,jmax
            jj = (j-1)*imax
            DO i=1,imax
               iar14t(i,j,l) = grds(i,j,2+l)        ! Temp (virtual)
               iar15u(i,j,l) = grdv(i,j,l,1)        ! U component
               iar16v(i,j,l) = grdv(i,j,l,2)        ! V component
                                                    ! specific humidity
               iar17q(i,j,l) = max(0.0,grds(i,j,2+kmax+l)*1.0e6)
               iarpsl(i,j,l) = wrk1(i+jj,l)*0.01    ! 3D layer pres(hPa)
            ENDDO
         ENDDO
      ENDDO
      DO l=1,kmax+1
         DO j=1,jmax
            jj = (j-1)*imax
            DO i=1,imax
               iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
                                                 ! (hPa)
            ENDDO
         ENDDO
      ENDDO
 
 
ccccc print *,' iar14t=',iar14t(1,80,:)
ccccc print *,' iar15u=',iar15u(1,80,:)
ccccc print *,' iar16v=',iar16v(1,80,:)
ccccc print *,' iarpsi=',iarpsi(1,80,:)
ccccc print *,' iarpsl=',iarpsl(1,80,:)
 
      DEALLOCATE (cofs, cofv)
      DEALLOCATE (grds, grdv, wrk1, wrk2)
 
      print *,' RETURNING from GBLENV10'
 
      RETURN
 
      SUBROUTINE gblevn10(IUNITF,IDATEP,IM,JM,IDRT) ! FORMERLY …
      END
C***********************************************************************
C***********************************************************************
C>    North-south swap
      subroutine gblevn11(imax,jmax,grid) ! formerly subroutine n_s_swap
      implicit none
      integer imax, jmax
      real grid(imax,jmax)
      real temp (imax)
      integer i, j, jj
 
      do j=1,jmax/2
        jj = jmax-j+1
        do i=1,imax
          temp(i)    = grid(i,j)
          grid(i,j)  = grid(i,jj)
          grid(i,jj) = temp(i)
        enddo
      enddo
      return
      subroutine gblevn11(imax,jmax,grid) ! formerly subroutine n_s_swap …
      end
C***********************************************************************
C***********************************************************************
C>    Does something   
      subroutine gblevn11d(imax,jmax,grid)
      implicit none
      integer imax, jmax
      real(kind=8) grid(imax,jmax)
      real(kind=8) temp (imax)
      integer i, j, jj
 
      do j=1,jmax/2
        jj = jmax-j+1
        do i=1,imax
          temp(i)    = grid(i,j)
          grid(i,j)  = grid(i,jj)
          grid(i,jj) = temp(i)
        enddo
      enddo
      return
      subroutine gblevn11d(imax,jmax,grid) …
      end
!---------------------------------------------------------
C> Does something
C> @param iunitf 
C> @param idatep 
C> @param im 
C> @param jm
C> @param IUNITF   - 2-WORD ARRAY:
C>        For SIGIO input:
C>              - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
C>                IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
C>                AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
C>                IS A MULTIPLE OF 3
C>              - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
C>                OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
C>                THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
C>        For NEMSIO input:
C>              - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
C>                (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
C>                THE DATE IN IDATEP
C>              - WORD 2 - NOT USED, SHOULD BE A NULL FILE
C> @param IDATEP   - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH      
C> @param im
C> @param jm
C> @param IDRT INTEGER GRID IDENTIFIER
C>              (IDRT=4 FOR GAUSSIAN GRID,
C>              IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C>              IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
      SUBROUTINE gblevn12(IUNITF,IDATEP,IM,JM,IDRT)
!---------------------------------------------------------
!--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NEMSIO
!--AND SUBSEQUENTLY MODIFIED BY FANGLIN YANG AND RUSS TREADON
 
      USE gblevn_module
      USE nemsio_module
      USE nemsio_openclose
      USE nemsio_read
      USE nemsio_write
      use sigio_module
 
      IMPLICIT NONE
      INTEGER IUNITF(2), IDATEP, IM, JM, IDRT
      REAL(KIND=8),parameter:: con_rv=4.6150e+2,con_rd=2.8705e+2,
     $          fv=con_rv/con_rd-1.,oner=1.0,qmin=1.0e-10
      INTEGER*4, PARAMETER :: TEN=10
 
      INTEGER*4 IDRTNEMS
      INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC
 
      INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC,
     $ iromb,maxwv,idir,i,j,k,l,ii,jj
 
      INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT
      INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4
      REAL(8) tfac
      REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:)
      REAL,ALLOCATABLE :: tmp(:)
      TYPE(nemsio_gfile) :: GFILE
 
      INTEGER IDATE(8),JDATE(8)
      CHARACTER*20 CFILE2
      REAL FHOUR,RINC(5)
      REAL (KIND=4),allocatable :: psfc(:,:), tv(:,:,:),
     $                             wrk1(:,:), wrk2(:,:)
 
      imax  = im
      jmax  = jm
      imjm4 = im*jm
 
! no time interpolation needed for nemsio input since files are
!  available every hour - original logic to perform time interpolation
!  is flawed and has been removed in this revision
! --------------------------------------------------------------------
 
      print 331, mod(idatep,100)
  331 format(/' --> GBLEVENTS: ONLY ONE NEMSIO INPUT GLOBAL FILE IS ',
     $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ',
     $ 'INTERPOLATE'/16x,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR',
     $ ' (',i2.2,') IS NOT A MULTIPLE OF 3)'/)
 
C  GET VALID-TIME DATE OF NEMSIO INPUT FILE, ALSO READ HEADERS
C  -----------------------------------------------------------
      CALL nemsio_init()
 
      rinc  = 0
      idate = 0
      WRITE(cfile2,'("fort.",I2.2)') iunitf(1)
 
      CALL nemsio_open(gfile,trim(cfile2),'read',iret=iret)
      print *,' cfile_nems2=',cfile2,'nemsio open,iret=',iret
 
         idrtnems=idrt
         CALL nemsio_getfilehead(gfile,idate=idate7,
     &        jcap=jcap4,dimx=im4,dimy=jm4,
     &         dimz=dimz4,idvc=idvc4,ntrac=ntrac,idrt=idrtnems,
     &         nfhour=nfhour,nfminute=nfminute,nfsecondn=nfsecondn,
     &         nfsecondd=nfsecondd,idsl=idsl4,idvm=idvm4,iret=iret)
      jcap=jcap4
      idvc=idvc4
      idsl=idsl4
      idvm=idvm4
      imax=im4
      jmax=jm4
      kmax=dimz4
      if(idrt==-9999) idrt=4 ! default for gaussian grid
      idate(1:3)=idate7(1:3)
      idate(5)=idate7(4)
 
      ALLOCATE(vcoord(kmax+1,3))
      vcoord=0.0
 
      IF(nfsecondd/=0.AND. nfsecondd/=-9999) THEN
         fhour=real(nfhour,8)+real(nfminute/60.,8)+
     $    real(nfsecondn*1./(nfsecondd*360.),8)
      ELSE
         fhour=real(nfhour,8)+real(nfminute/60.,8)
      ENDIF
      print'("  idate=",I5,7I3.2,"  fhour=",F7.3)', idate,fhour
 
      ALLOCATE(vcoord4(kmax+1,3,2))
      CALL nemsio_getfilehead(gfile,vcoord=vcoord4,iret=iret)
      nvcoord=3
      IF(maxval(vcoord4(:,3,1))==0..AND.minval(vcoord4(:,3,1))==0.) THEN
         nvcoord=2
         IF(maxval(vcoord4(:,2,1))==0..AND.minval(vcoord4(:,2,1))==0.)
     $    nvcoord=1
      ENDIF
 
      vcoord(1:kmax+1,1:nvcoord)=vcoord4(1:kmax+1,1:nvcoord,1)
      DEALLOCATE(vcoord4)
 
      ALLOCATE(cpi(ntrac+1))
      CALL nemsio_getheadvar(gfile,'CPI',cpi,iret=iret)
 
      rinc(2) = fhour
      CALL w3movdat(rinc,idate,jdate)
 
      print 1, fhour,(idate(ii),ii=1,3),idate(5),(jdate(ii),ii=1,3),
     $ jdate(5)
    1 FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',f5.1,
     $ ' HOUR FORECAST FROM ',i5.4,3i3.2,' VALID AT ',i5.4,3i3.2)
 
      idatgs_cor = (jdate(1) * 1000000) + (jdate(2) * 10000) +
     $ (jdate(3) * 100) + jdate(5)
 
C  VALID DATES MUST MATCH
C  ----------------------
 
      IF(idatep.NE.idatgs_cor)  GO TO 901
 
C------------------------------------------
 
      sfcpress_id  = mod(idvm,ten)
      thermodyn_id = mod(idvm/ten,ten)
 
      IF(idvc == 3 .AND. thermodyn_id == 3) THEN
         kmaxs = (ntrac+1)*kmax + 2
      ELSE
         kmaxs = 2*kmax + 2
         ntrac = 1
      ENDIF
 
      ALLOCATE (iar12z(imax,jmax), iar13p(imax,jmax))
      ALLOCATE (iar14t(imax,jmax,kmax),  iar15u(imax,jmax,kmax),
     $          iar16v(imax,jmax,kmax),  iar17q(imax,jmax,kmax),
     $          iarpsl(imax,jmax,kmax),  iarpsi(imax,jmax,kmax+1),
     $          iarpsd(imax,jmax,kmax))
 
!!! DAK: is this right???????
      if(idvc.eq.0)  idvc = 1  ! Reset IDVC=0 to 1 (sigma coord.)
      IF(idvc < 0 .or. idvc > 3) THEN
         print *, '##GBLEVENTS/GBLEVN12: INVALID VERT COORD ID (=',idvc
      ENDIF
 
 
C  DEFINE THE OTHER RESOLUTION PARAMETERS
C  --------------------------------------
 
      jcap1   = jcap+1
      jcap2   = jcap+2
      jcap1x2 = jcap1*2
      mdima   = jcap1*jcap2
      mdimb   = mdima/2+jcap1
      mdimc   = mdimb*2
 
      dlat  = 180./(jmax-1)
      dlon  = 360./imax
 
      print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,idvc
 
! DAK: is below correct??????????
    2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
     $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
     $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
     $ 'COORD ID IS: ',i0,' (not sure what this means!')
 
C***********************************************************************
 
      iromb = 0
      maxwv = jcap
      if (idrt < 0 .or. idrt > 256) idrt  = 0
      idir  = 1
 
      print *,'in nem sfcpress_id=',sfcpress_id,' thermodyn_id=',
     $ thermodyn_id,' ntrac=',ntrac
      print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=',
     $ nvcoord
 
      call w3kind(kindreal,kindint)
 
!     print *,'in gblevent, get fields,imax,jmax
      allocate(tmp(imax*jmax))
!
! hgt
      k4=1
      if(kindreal==4) then
         CALL nemsio_readrecvw34(gfile,'hgt','sfc',k4,tmp,iret=iret)
      else
         CALL nemsio_readrecv(gfile,'hgt','sfc',k4,tmp,iret=iret)
      endif
!     print *,'in getnemsio,hgtsfc,',maxval(tmp),minval(tmp),iret
      if(iret==0) THEN
         iar12z(:,:)=reshape(tmp,(/imax,jmax/))
         call gblevn11d(imax,jmax,iar12z)
      ENDIF
 
!sfc pres, input in Pa, out in hPa
      if(kindreal==4) then
         CALL nemsio_readrecvw34(gfile,'pres','sfc',k4,tmp,iret=iret)
      else
         CALL nemsio_readrecv(gfile,'pres','sfc',k4,tmp,iret=iret)
      endif
!     print *,'in getnemsio,pressfc,',maxval(tmp),minval(tmp),iret
      if(iret==0) THEN
         iar13p(:,:)=reshape(tmp*0.01,(/imax,jmax/))
         call gblevn11d(imax,jmax,iar13p)
      ENDIF
!dpres
!     DO K4=1,KMAX
!        if(kindreal==4) then
!           CALL NEMSIO_READRECVw34(GFILE,'dpres','mid layer',K4,tmp,
!    $       iret=iret)
!        else
!           CALL NEMSIO_READRECV(GFILE,'dpres','mid layer',K4,tmp,
!    $       iret=iret)
!        endif
!        print *,'in getnemsio,dpres,k4=',k4,maxval(tmp),minval(tmp),
!    $    iret
!        if(iret==0) THEN
!           IARPSI(:,:,K4+1)=reshape(tmp*0.01,(/imax,jmax/))
!        ENDIF
!     ENDDO
 
!pres
!     DO K4=1,KMAX
!        if(kindreal==4) then
!           CALL NEMSIO_READRECVw34(GFILE,'pres','mid layer',K4,tmp,
!    $       iret=iret)
!        else
!           CALL NEMSIO_READRECV(GFILE,'pres','mid layer',K4,tmp,
!    $       iret=iret)
!        endif
!        print *,'in getnemsio,pres,k4=',k4,maxval(tmp),minval(tmp),iret
!        if(iret==0) THEN
!           IARPSL(:,:,k4)=reshape(tmp*0.01,(/imax,jmax/))
!        ENDIF
!     ENDDO
! tmp
      DO k4=1,kmax
         if(kindreal==4) then
            CALL nemsio_readrecvw34(gfile,'tmp','mid layer',k4,tmp,
     $       iret=iret)
         else
            CALL nemsio_readrecv(gfile,'tmp','mid layer',k4,tmp,
     $       iret=iret)
         endif
!        print *,'in getnemsio,tmp,k4=',k4,maxval(tmp),minval(tmp),iret
         if(iret==0) THEN
            iar14t(:,:,k4)=reshape(tmp,(/imax,jmax/))
            call gblevn11d(imax,jmax,iar14t(1,1,k4))
         ENDIF
      ENDDO
!u
      DO k4=1,kmax
         if(kindreal==4) then
            CALL nemsio_readrecvw34(gfile,'ugrd','mid layer',k4,tmp,
     $       iret=iret)
         else
            CALL nemsio_readrecv(gfile,'ugrd','mid layer',k4,tmp,
     $       iret=iret)
         endif
!        print *,'in getnemsio,utmp,k4=',k4,maxval(tmp),minval(tmp),iret
         if(iret==0) THEN
            iar15u(:,:,k4)=reshape(tmp,(/imax,jmax/))
            call gblevn11d(imax,jmax,iar15u(1,1,k4))
         ENDIF
      ENDDO
!     do k=1,kmax
!        print *,'in getnemsio,i15u,k=',k,maxval(iar15u(:,:,k)),
!    $    minval(iar15u(:,:,k)),iret
!     end do
!     print *,'in getnemsio,IAR15U.1,',maxval(IAR15U),minval(IAR15U)
 
!v
      DO k4=1,kmax
         if(kindreal==4) then
            CALL nemsio_readrecvw34(gfile,'vgrd','mid layer',k4,tmp,
     $       iret=iret)
         else
            CALL nemsio_readrecv(gfile,'vgrd','mid layer',k4,tmp,
     $       iret=iret)
         endif
         if(iret==0) THEN
            iar16v(:,:,k4)=reshape(tmp,(/imax,jmax/))
            call gblevn11d(imax,jmax,iar16v(1,1,k4))
         ENDIF
      ENDDO
!     print *,'in getnemsio,IAR15U.2,',maxval(IAR15U),minval(IAR15U)
!q
      DO k4=1,kmax
         if(kindreal==4) then
            CALL nemsio_readrecvw34(gfile,'spfh','mid layer',k4,tmp,
     $       iret=iret)
         else
            CALL nemsio_readrecv(gfile,'spfh','mid layer',k4,tmp,
     $       iret=iret)
         endif
         if(iret==0) THEN
            iar17q(:,:,k4)=max(0.0,reshape(tmp*1.0e6,(/imax,jmax/)) )
            call gblevn11d(imax,jmax,iar17q(1,1,k4))
         endif
      ENDDO
!     print *,'in getnemsio,IAR15U.3,',maxval(IAR15U),minval(IAR15U)
 
!
      deallocate(tmp)
      CALL nemsio_close(gfile,iret=iret)
!     print *,'end of nemsio'
!
 
!     print *,'in getnemsio,IAR15U.4,',maxval(IAR15U),minval(IAR15U)
 
!--compute Tv from temp
      DO k=1,kmax
         DO j=1,jmax
            DO i=1,imax
               tfac=oner+fv*max(iar17q(i,j,k)*1.0e-6,qmin)
               iar14t(i,j,k)=iar14t(i,j,k)*tfac
            ENDDO
         ENDDO
      ENDDO
 
!     print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U)
 
!--COMPUTE PM AND PD
      ALLOCATE (psfc(imax,jmax), tv(imax,jmax,kmax))
      ALLOCATE (wrk1(imax*jmax,kmax),  wrk2(imax*jmax,kmax+1))
 
      imjm4=imax*jmax
      km4=kmax
      psfc(:,:) = iar13p(:,:)*100.
      tv(:,:,:) = iar14t(:,:,:)
 
      IF(thermodyn_id == 3 .AND. idvc == 3) THEN
         tv(:,:,:) = tv(:,:,:) / cpi(1)
         print *,' cpi(1)=',cpi(1)
      ENDIF
 
      CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
     $ psfc,tv,pm=wrk1,pd=wrk2(1,2))
      DO j=1,jmax
         jj = (j-1)*imax
         DO i=1,imax
            wrk2(i+jj,1) = psfc(i,j)                       ! in Pa
         ENDDO
      ENDDO
      DO l=1,kmax
         wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1)             ! in Pa
      ENDDO
 
      DO l=1,kmax
         DO j=1,jmax
            jj = (j-1)*imax
            DO i=1,imax
               iarpsl(i,j,l) = wrk1(i+jj,l)*0.01    ! 3D layer pres(hPa)
            ENDDO
         ENDDO
      ENDDO
      DO l=1,kmax+1
         DO j=1,jmax
            jj = (j-1)*imax
            DO i=1,imax
               iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
                                                 !  (hPa)
            ENDDO
         ENDDO
      ENDDO
 
!     print*,'GBLEVN12 IARPSI,',maxval(IARPSI),minval(IARPSI)
!     print*,'GBLEVN12 IARPSL,',maxval(IARPSL),minval(IARPSL)
 
      CALL nemsio_finalize()
!
      CALL getlats(idrt)
 
      RETURN
 
  901 CONTINUE
      print 9901, (jdate(ii),ii=1,3),jdate(5),idatep
 9901 FORMAT(/' ##GBLEVENTS/GBLEVN12 - NEMSIO INPUT GLOBAL FILE DATE',
     $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ',
     $ 'DATE (',i10,') -STOP 68'/)
      CALL errexit(68)
 
      SUBROUTINE gblevn12(IUNITF,IDATEP,IM,JM,IDRT) …
      END
! --------------------------------------------------------------------
C> NetCDF Input
C> @param IUNITF   - 2-WORD ARRAY:
C>        For SIGIO input:
C>              - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
C>                IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
C>                AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
C>                IS A MULTIPLE OF 3
C>              - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
C>                (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
C>                ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
C>                THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
C>                OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
C>                THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
C>        For NEMSIO input:
C>              - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
C>                (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
C>                THE DATE IN IDATEP
C>              - WORD 2 - NOT USED, SHOULD BE A NULL FILE
C> @param IDATEP   - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH     
C> @param IMX
C> @param JMX
C> @param IDRT INTEGER GRID IDENTIFIER
C>                (IDRT=4 FOR GAUSSIAN GRID,
C>                 IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C>                 IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
      SUBROUTINE gblevn13(IUNITF,IDATEP,IMX,JMX,IDRT) 
!---------------------------------------------------------
!--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NETCDF INPUT
                                                   
      USE gblevn_module
      use sigio_module
      use netcdf
 
      IMPLICIT NONE
      INTEGER IUNITF(2), IDATEP, IDRT,IMX,JMX
      INTEGER yyyy,mm,dd,hh
      INTEGER*4 ncid
      integer*4  error, id_var, dimid, len
      integer*4  im,jm,km, lm,n, k,nargs
      REAL(KIND=8),parameter:: con_rv=4.6150e+2,con_rd=2.8705e+2,
     $          fv=con_rv/con_rd-1.,oner=1.0,qmin=1.0e-10         
      INTEGER*4, PARAMETER :: TEN=10
 
      INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC
 
      INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC,
     $ iromb,maxwv,idir,i,j,l,ii,jj
 
      INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT
      INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4,kr
      REAL(8) tfac, time
      REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:),ak(:),bk(:)
      REAL,ALLOCATABLE :: temp(:,:), temp3d(:,:,:)
!      CHARACTER,ALLOCATABLE :: attone(:)
      character (len = 80) :: attone
      INTEGER nt1, nt2
      character(len=10) :: dim_nam, grid
      CHARACTER*20 gfname
 
      INTEGER(4) IDATE(8),JDATE(8)
      REAL(4) FHOUR,RINC(5)
 
      REAL (KIND=4),allocatable :: psfc(:,:), tv(:,:,:),
     $                             wrk1(:,:), wrk2(:,:)
 
      imax  = im
      jmax  = jm
 
! no time interpolation needed for netcdf input since files are
!  available every hour - original logic to perform time interpolation
!  is flawed and has been removed in this revision
! --------------------------------------------------------------------
 
      print 331, mod(idatep,100)
  331 format(/' --> GBLEVENTS: ONLY ONE NETCDF INPUT GLOBAL FILE IS ',
     $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ',
     $ 'INTERPOLATE'/16x,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR',
     $ ' (',i2.2,') IS NOT A MULTIPLE OF 3)'/)
 
C  GET VALID-TIME DATE OF NETCDF INPUT FILE, ALSO READ HEADERS
C  -----------------------------------------------------------
 
      WRITE(gfname,'("fort.",I2.2)') iunitf(1)
 
       error=nf90_open(trim(gfname),nf90_nowrite,ncid)
       error=nf90_inq_dimid(ncid,"grid_xt",dimid)
       error=nf90_inquire_dimension(ncid,dimid,dim_nam,im)
       error=nf90_inq_dimid(ncid,"grid_yt",dimid)
       error=nf90_inquire_dimension(ncid,dimid,dim_nam,jm)
       error=nf90_inq_dimid(ncid,"pfull",dimid)
       error=nf90_inquire_dimension(ncid,dimid,dim_nam,km)
       error=nf90_inq_dimid(ncid,"phalf",dimid)
       error=nf90_inquire_dimension(ncid,dimid,dim_nam,lm)
       print*, "im,jm,kmi,lm:",im,jm,km,lm
       print*, "IMX, JMX:", imx, jmx
       IF (imx .NE. im .OR. jmx .NE. jm) print*, "Different Resolutions"
       imax = im
       jmax = jm
       kmax = km
 
!      error=nf90_inq_varid(ncid, 'pfull', id_var)
!      error=nf90_get_var(ncid, id_var, )      
 
 
!      CALL W3MOVDAT(RINC,IDATE,JDATE)
 
 
 
C  VALID DATES MUST MATCH
C  ----------------------
      error=nf90_inq_varid(ncid, "time", id_var)
      error=nf90_inquire_attribute(ncid, id_var, "units", len=len)
!      ALLOCATE (attone(len))
      error=nf90_get_att(ncid, id_var, "units", attone)
        read(attone(len-18:len-15),'(i4)') yyyy
        read(attone(len-13:len-12),'(i2)') mm
        read(attone(len-10:len-9),'(i2)') dd
        read(attone(len-7:len-6),'(i2)') hh
      IF(attone(1:5) .NE. "hours") THEN
      print*, "checkout the time unit, not hourly",attone
      ELSE
      print*, "base time", yyyy,mm,dd,hh
      ENDIF
!      DEALLOCATE (attone)
       error=nf90_get_var(ncid, id_var, time)
      fhour=time
 
      idate = 0
      idate(1)=yyyy
      idate(2)=mm
      idate(3)=dd
      idate(5)=hh
 
      rinc=0.0
      rinc(2)=fhour
 
      CALL w3movdat(rinc,idate,jdate)
 
      print 1, fhour,(idate(ii),ii=1,3),idate(5),(jdate(ii),ii=1,3),
     $ jdate(5)
 1     FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',f5.1,
     $ ' HOUR FORECAST FROM ',i5.4,3i3.2,' VALID AT ',i5.4,3i3.2)
 
      idatgs_cor = (jdate(1) * 1000000) + (jdate(2) * 10000) +
     $ (jdate(3) * 100) + jdate(5)
 
C  VALID DATES MUST MATCH
C  ----------------------
 
      IF(idatep.NE.idatgs_cor)  GO TO 901
 
C------------------------------------------
 
      ALLOCATE (iar12z(im,jm), iar13p(im,jm))
      ALLOCATE (iar14t(im,jm,km),  iar15u(im,jm,km),
     $          iar16v(im,jm,km),  iar17q(im,jm,km),
     $          iarpsl(im,jm,km),  iarpsi(im,jm,km+1),
     $          iarpsd(im,jm,km))
 
      error=nf90_get_att(ncid, nf90_global, "grid", grid)
      IF (grid == "gaussian")THEN
        idrt=4
      ENDIF
 
      nvcoord=2     !for idvc=2, nvcoord=2: hybrid interface a and b
      idvc=2        !(1 for sigma and 2 for hybrid)
      idsl=1        !(1 for phillips or 2 for mean)
      idvm=1
      jcap = -9999
 
      sfcpress_id  = mod(idvm,ten)
      thermodyn_id = mod(idvm/ten,ten)
 
      IF(idvc == 3 .AND. thermodyn_id == 3) THEN
         kmaxs = (ntrac+1)*kmax + 2
      ELSE
         kmaxs = 2*kmax + 2
         ntrac = 1
      ENDIF
 
C  DEFINE THE OTHER RESOLUTION PARAMETERS
C  --------------------------------------
 
      dlat  = 180./(jmax-1)
      dlon  = 360./imax
 
      print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,idvc
 
 2     FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
     $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
     $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
     $ 'COORD ID IS: ',i0)
 
      print *,'in netcdf sfcpress_id=',sfcpress_id,' thermodyn_id=',
     $ thermodyn_id,' ntrac=',ntrac
      print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=',
     $ nvcoord
 
      ALLOCATE(vcoord(km+1,3))
      vcoord=0.0
      allocate(ak(km+1))
      allocate(bk(km+1))
 
      error=nf90_get_att(ncid, nf90_global, "ak", ak)
      error=nf90_get_att(ncid, nf90_global, "bk", bk)
 
      do k=1,km+1
         kr=km+2-k
         vcoord(k,1) = ak(kr)
         vcoord(k,2) = bk(kr)
      end do
      
      deallocate(ak,bk)
 
      allocate(temp(im,jm))
      error=nf90_inq_varid(ncid, 'hgtsfc', id_var)
      error=nf90_get_var(ncid, id_var, temp)
      
      iar12z(:,:)=temp(:,:)
      call gblevn11d(imax,jmax,iar12z)
      
 
      error=nf90_inq_varid(ncid, 'pressfc', id_var)
      error=nf90_get_var(ncid, id_var, temp)
 
      iar13p(:,:)=temp(:,:)*0.01
      call gblevn11d(imax,jmax,iar13p)
 
      deallocate(temp)
 
!!      write(6,111)'hgtsfc', 1,minval(iar12z),maxval(iar12z)
!!      write(6,111)'pressfc',1,minval(iar13p),maxval(iar13p)
!! 111  format(a8,1x,i3,1x,2(g13.6,1x))
 
 
      allocate(temp3d(im,jm,km))
      error=nf90_inq_varid(ncid, 'tmp', id_var)
      error=nf90_get_var(ncid, id_var, temp3d)
      DO k4=1,km
         kr=km+1-k4
         iar14t(:,:,k4)=temp3d(:,:,kr)
         call gblevn11d(imax,jmax,iar14t(1,1,k4))
!!         write(6,111)'tmp',k4,minval(iar14t(:,:,k4)),
!!     &        maxval(iar14t(:,:,k4))
      ENDDO
 
      error=nf90_inq_varid(ncid, 'ugrd', id_var)
      error=nf90_get_var(ncid, id_var, temp3d)
      DO k4=1,km
         kr=km+1-k4
         iar15u(:,:,k4)=temp3d(:,:,kr)
         call gblevn11d(imax,jmax,iar15u(1,1,k4))
!!         write(6,111)'ugrd',k4,minval(iar15u(:,:,k4)),
!!     &        maxval(iar15u(:,:,k4))
      ENDDO
 
      error=nf90_inq_varid(ncid, 'vgrd', id_var)
      error=nf90_get_var(ncid, id_var, temp3d)
      DO k4=1,km
         kr=km+1-k4
         iar16v(:,:,k4)=temp3d(:,:,kr)
         call gblevn11d(imax,jmax,iar16v(1,1,k4))
!!         write(6,111)'vgrd',k4,minval(iar16v(:,:,k4)),
!!     &        maxval(iar16v(:,:,k4))
      ENDDO
 
      error=nf90_inq_varid(ncid, 'spfh', id_var)
      error=nf90_get_var(ncid, id_var, temp3d)
      DO k4=1,km
         kr=km+1-k4
         iar17q(:,:,k4)=max(0.0,temp3d(:,:,kr)*1.e6)
         call gblevn11d(imax,jmax,iar17q(1,1,k4))
!!         write(6,111)'spfh',k4,minval(iar17q(:,:,k4)),
!!     &        maxval(iar17q(:,:,k4))
      ENDDO
 
      deallocate(temp3d)
!
 
 
!--compute Tv from temp
      DO k=1,km
         DO j=1,jm
            DO i=1,im
               tfac=oner+fv*max(iar17q(i,j,k)*1.0e-6,qmin)
               iar14t(i,j,k)=iar14t(i,j,k)*tfac
            ENDDO
         ENDDO
      ENDDO
 
!     print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U)
 
 
 
!--COMPUTE PM AND PD
      ALLOCATE (psfc(im,jm), tv(im,jm,km))
      ALLOCATE (wrk1(im*jm,km),  wrk2(im*jm,km+1))
 
      imjm4=im*jm
      km4=km
      psfc(:,:) = iar13p(:,:)*100.
      tv(:,:,:) = iar14t(:,:,:)
 
      IF(thermodyn_id == 3 .AND. idvc == 3) THEN
         tv(:,:,:) = tv(:,:,:) / cpi(1)
         print *,' cpi(1)=',cpi(1)
      ENDIF
 
      CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
     $ psfc,tv,pm=wrk1,pd=wrk2(1,2))
      DO j=1,jm
         jj = (j-1)*im
         DO i=1,im
            wrk2(i+jj,1) = psfc(i,j)                       ! in Pa
         ENDDO
      ENDDO
      DO l=1,km
         wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1)             ! in Pa
      ENDDO
 
      DO l=1,km
         DO j=1,jm
            jj = (j-1)*im
            DO i=1,im
               iarpsl(i,j,l) = wrk1(i+jj,l)*0.01    ! 3D layer pres(hPa)
            ENDDO
         ENDDO
      ENDDO
      DO l=1,km+1
         DO j=1,jm
            jj = (j-1)*im
            DO i=1,im
               iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
                                                 !  (hPa)
            ENDDO
         ENDDO
      ENDDO
      
      error=nf90_close(ncid)        
!     print*,'GBLEVN13 IARPSI,',maxval(IARPSI),minval(IARPSI)
!     print*,'GBLEVN13 IARPSL,',maxval(IARPSL),minval(IARPSL)
 
!
      CALL getlats(idrt)
 
      RETURN
 
  901 CONTINUE
      print 9901, idatep,idatgs_cor
 9901 FORMAT(/' ##GBLEVENTS/GBLEVN13 - NETCDF INPUT GLOBAL FILE DATE',
     $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ',
     $ 'DATE (',i10,') -STOP 68'/)
      CALL errexit(68)
 
      SUBROUTINE gblevn13(IUNITF,IDATEP,IMX,JMX,IDRT)  …
      END
 
!-----------------------------------------------------------------
!-----------------------------------------------------------------
C> get latitudes
C> @param IDRT INTEGER GRID IDENTIFIER
C>               (IDRT=4 FOR GAUSSIAN GRID,
C>               IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C>               IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
      subroutine getlats(idrt)
 
      USE gblevn_module
!      integer :: jmax
!      real,allocatable :: slat(:),wlat(:)
!      real :: rad2deg
      real,allocatable :: slats(:)
      real(4) slat(jmax),wlat(jmax),rad2deg
!get gaussian or regular latitude array based on idrt
 
!      print *,'in getlats,idrt=',idrt
      call splat(idrt,jmax,slat,wlat)
      rad2deg=180./acos(-1.)
      allocate(slats(jmax));
      rad2deg=180./acos(-1.)
      slats(:)=-asin(slat(:))*rad2deg
      dlat=180./float(jmax-1)
 
      return
      subroutine getlats(idrt) …
      end