gblevents  1.8.0
•All Files Functions Pages
gblevents.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Pre/Post processing of of prepbufr events
3  MODULE gblevn_module
4 
5  IMPLICIT NONE
6 
7  SAVE
8 
9  INTEGER IMAX,JMAX,KMAX,KMAXS
10  INTEGER*4 IDVC,IDSL,NVCOORD,SFCPRESS_ID,THERMODYN_ID
11  REAL (KIND=8), ALLOCATABLE :: IAR14T(:,:,:), IAR15U(:,:,:),
12  $ IAR16V(:,:,:), IAR17Q(:,:,:),
13  $ IAR12Z(:,:), IAR13P(:,:),
14  $ IARPSI(:,:,:), IARPSL(:,:,:),
15  $ IARPSD(:,:,:)
16  REAL (KIND=4), ALLOCATABLE :: VCOORD(:,:)
17  REAL DLAT,DLON
18 
19  END MODULE gblevn_module
20 
21 C> SUBPROGRAM: GBLEVENTS PRE/POST PROCESSING OF PREPBUFR EVENTS
22 C> PRGMMR: DENNIS KEYSER ORG: EMC DATE: 2017-02-22
23 C>
24 C> RUNS IN TWO MODES: "PREVENTS" AND "POSTEVENTS". IN THE
25 C> PREVENTS MODE, PREPARES OBSERVATIONAL PREPBUFR REPORTS FOR
26 C> SUBSEQUENT QUALITY CONTROL AND ANALYSIS PROGRAMS. THIS IS DONE
27 C> THROUGH THE FOLLOWING: INTERPOLATION OF GLOBAL FIRST GUESS {FROM
28 C> EITHER SIGIO (SIGMA OR HYBRID) OR NEMSIO INPUT} TO PREPBUFR
29 C> OBSERVATION LOCATIONS WITH ENCODING OF FIRST GUESS VALUES INTO
30 C> PREPBUFR REPORTS; ENCODING OF "PREVENT" AND/OR "VIRTMP" EVENTS INTO
31 C> PREPBUFR REPORTS; IN CERTAIN CASES, ENCODING OF A DERIVED PMSL INTO
32 C> SURFACE PREPBUFR REPORTS; AND ENCODING OF OBSERVATION ERRORS FROM
33 C> THE ERROR SPECIFICATION FILE INTO PREPBUFR REPORTS. IN THE
34 C> POSTEVENTS MODE, AFTER ALL QUALITY CONTROL AND ANALYSIS PROGRAMS
35 C> HAVE RUN, INTERPOLATES THE GLOBAL ANALYSIS {FROM EITHER SIGIO
36 C> (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT} TO PREPBUFR OBSERVATION
37 C> LOCATIONS AND ENCODES THESE ANALYZED VALUES INTO PREPBUFR REPORTS.
38 C> THE REMAINDER OF THIS ABSTRACT APPLIES ONLY TO THE PREVENTS MODE.
39 C> THE "PREVENT" EVENT CAN CHANGE A QUALITY MARKER TO FLAG AN
40 C> OBSERVATION DATUM FOR NON-USE BY SUBSEQUENT QC AND ANALYSIS
41 C> PROGRAMS (FILTERING). EXAMPLES WHERE THIS SUBROUTINE WILL WRITE
42 C> AN EVENT TO FLAG A DATUM INCLUDE: THE OBSERVATION ERROR FOR THAT
43 C> DATUM IS READ IN AS MISSING IN THE INPUT ERROR FILE, THE DATUM
44 C> ITSELF VIOLATES A GROSS OR "SANITY" CHECK, OR THE OBSERVED
45 C> PRESSURE DATUM IS MORE THAN 100 MB BELOW THE GUESS SURFACE
46 C> PRESSURE. THE "VIRTMP" EVENT CAN CHANGE THE SPECIFIC HUMIDITY
47 C> OBSERVATION (RE-CALCULATED) AS WELL AS THE TEMPERATURE
48 C> OBSERVATION (FROM SENSIBLE TO VIRTUAL TEMPERATURE, BASED ON
49 C> JUST-CALCULATED SPECIFIC HUMIDITY). CURRENTLY THIS APPLIES ONLY
50 C> TO SURFACE (LAND, MARINE AND MESONET) DATA TYPES, POSSIBLY TO
51 C> RAOB, DROP AND MULTI-LEVEL RECCO DATA TYPES IF THE SWITCH
52 C> "ADPUPA_VIRT" IS TRUE (NORMALLY, HOWEVER IT IS FALSE) [OTHER DATA
53 C> TYPES WITH REPORTED SENSIBLE TEMPERATURE EITHER HAVE MISSING
54 C> MOISTURE (E.G., ALL AIRCRAFT TYPES EXCEPT FOR SOME ACARS, SATELLITE
55 C> WIND TYPES), FLAGGED MOISTURE (E.G., SOME ACARS) OR CALCULATE
56 C> SPECIFIC HUMIDITY/VIRTUAL TEMPERATURE IN SUBSEQUENT PROGRAMS (E.G.,
57 C> RAOBS, DROPS AND MULTI-LEVEL RECCOS WHICH CALCULATE THESE IN
58 C> PROGRAM "CQCBUFR", IN WHICH CASE THE SWITCH "ADPUPA_VIRT" HERE MUST
59 C> BE FALSE!)]. FOR CASES WHERE THE SWITCH "DOBERR" IS FALSE, THE
60 C> OBSERVATION ERROR FOR ALL DATA REMAINS MISSING IN THE PREPBUFR
61 C> FILE. IN THIS CASE, THE INPUT ERROR FILE IS USUALLY A NULL FILE
62 C> AND THE "PREVENT" EVENT TO FLAG THE DATUM IS NOT INVOKED. FOR
63 C> CASES WHERE THE SWITCH "DOFCST" IS FALSE, IF THE SWITCH "SOME_FCST"
64 C> IS ALSO FALSE, THEN FORECAST VALUES ARE NOT ENCODED FOR ANY MESSAGE
65 C> TYPE; IF "SOME_FCST" IS TRUE THEN FORECAST VALUES ARE ENCODED, BUT
66 C> ONLY FOR REPORTS IN THOSE MESSAGE TYPES FOR WHICH A GUESS VALUE IS
67 C> NEEDED BY SUBSEQUENT QC PROGRAMS. IT SHOULD BE NOTED THAT THE
68 C> FILTERING OF DATA ASSOCIATED WITH THE "PREVENT" EVENT PROCESSING IS
69 C> NOT INVOKED IF ALL THREE ARE TRUE: DOBERR= FALSE, THE FORECAST
70 C> VALUES ARE MISSING (DOFCST=FALSE & SOME_FCST=TRUE & MESSAGE TYPE IS
71 C> NOT "ADPUPA", "AIRCFT", "AIRCAR", "PROFLR", OR "VADWND" -- OR --
72 C> DOFCST=FALSE & SOME_FCST=FALSE), AND "VIRTMP" EVENT PROCESSING IS
73 C> NOT INVOKED (EITHER MESSAGE TYPE IS NOT "ADPSFC", "SFCSHP" OR
74 C> "MSONET" WHEN "ADPUPA_VIRT" IS FALSE, OR MESSAGE TYPE IS NOT
75 C> "ADPSFC", "SFCSHP", "MSONET" OR "ADPUPA" WHEN "ADPUPA_VIRT" IS
76 C> TRUE). ALSO, IF VIRTUAL TEMPERATURE PROCESSING IS PERFORMED, ALL
77 C> SURFACE REPORTS WITH MISSING PMSL WILL ENCODE A DERIVED PMSL INTO
78 C> PREPBUFR IF THE SWITCH DOPMSL IS TRUE AND A VIRTUAL TEMPERATURE WAS
79 C> SUCCESSFULLY CALCULATED.
80 C>
81 C> PROGRAM HISTORY LOG:
82 C> -1999-07-01 D. A. KEYSER -- ORIGINAL AUTHOR (ADAPTED FROM PREVENTS
83 C> SUBROUTINE IN PREPDATA PROGRAM, BUT NOW GENERALIZED FOR
84 C> POSTEVENTS MODE)
85 C> -1999-07-12 D. A. KEYSER -- MODIFIED TO INTERPOLATE MODEL SPECIFIC
86 C> HUMIDITY TO OBSERVATION LOCATION WHEN OBS. SPECIFIC HUMIDITY IS
87 C> MISSING AS LONG AS OBS. TEMPERATURE IS NON-MISSING
88 C> -1999-09-09 D. A. KEYSER -- ADDED "VADWND" TO THE LIST OF MESSAGE
89 C> TYPES FOR WHICH FORECAST VALUES MUST BE ENCODED, EVEN WHEN
90 C> DOFCST=FALSE (NECESSARY BECAUSE THE NEW PROGRAM CQCVAD NEEDS THE
91 C> BACKGROUND DATA)
92 C> -1999-09-09 D. A. KEYSER -- CHANGES TO MAKE CODE MORE PORTABLE;
93 C> 'TFC' NOW GENERATED FOR VADWND MESSAGE TYPES EVEN THOUGH TOB IS
94 C> MISSING (NEEDED BY CQCVAD PROGRAM)
95 C> -1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE
96 C> NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT
97 C> TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD
98 C> QUAL. MARKER (8)
99 C> -2000-09-21 D. A. KEYSER -- THE PRESSURE LEVEL ABOVE WHICH ALL SPEC.
100 C> HUMIDITY QUAL. MARKERS ARE "REJECTED" (Q.M. SET TO 9) IS NOW READ
101 C> IN AS A N-LIST SWITCH (QTOP_REJ), BEFORE IT WAS HARDWIRED TO 300
102 C> MB
103 C> -2000-12-13 D. A. KEYSER -- WILL NO LONGER PERFORM VIRTUAL TEMPERATURE
104 C> PROCESSING FOR ACARS DATA SINCE MOISTURE IS FLAGGED RIGHT NOW
105 C> (ACARS MOISTURE ONLY WRITTEN INTO PREPBUFR FILE FOR STATISTICAL
106 C> REASONS)
107 C> -2001-02-02 D. A. KEYSER -- RESTORED LEGACY LOGIC TO FLAG CERTAIN
108 C> SATELLITE TEMPERATURE SOUNDINGS EITHER BELOW 100 MB (TEMP. OBS)
109 C> OR ON ALL LEVELS (SPEC. HUM. OBS), CONTROLLED BY NEW NAMELIST
110 C> SWITCH "SATMQC"
111 C> -2001-09-27 D. A. KEYSER -- 'TFC' AND 'QFC' NOW GENERATED FOR REPORT
112 C> TYPE 111 (SYNDAT REPORTS AT STORM CENTER) EVEN THOUGH "TOB" AND
113 C> "QOB" ARE MISSING (NEEDED BY SYNDATA PROGRAM); IN PREPARATION FOR
114 C> CHANGE FROM T170L42 TO T254L64 SGES, NOW MAKES COEFFICIENT ARRAYS
115 C> ALLOCATABLE TO ALLOW THEM TO OBTAIN MEMORY FROM "HEAP" RATHER
116 C> THAN FROM "STACK", ALSO HAVE INCREASED THE MAX NUMBER OF LEVELS
117 C> IN ARRAYS FROM 42 TO 64, FINALLY ALSO NO LONGER STOPS WITH C.
118 C> CODE 70 IF EVEN NUMBER OF LONGITUDES IN SIGMA GUESS (IMAX,
119 C> HARDWIRED TO 384) IS .LT. SPECTRAL RESOLUTION (JCAP) * 2
120 C> -2001-10-10 D. A. KEYSER -- AT PREPBUFR CENTER DATES WITH AN HOUR THAT
121 C> IS NOT A MULTIPLE OF 3 (WHEN A GLOBAL SIGMA GUESS/ANAL FILE IS
122 C> NOT AVAILABLE; E.G., IN RUC2A RUNS) NOW PERFORMS A LINEAR
123 C> INTERPOLATION BETWEEN SPECTRAL COEFFICIENTS IN 2 SPANNING SIGMA
124 C> GUESS/ANAL FILES 3-HRS APART TO CENERATE A GUESS/ANAL FILE VALID
125 C> AT THE PREPBUFR CENTER TIME
126 C> -2002-05-10 D. A. KEYSER -- ADDED "AIRCAR" TO THE LIST OF TABLE A
127 C> MESSAGE TYPES THAT WILL STILL HAVE THE BACKGROUND ENCODED WHEN
128 C> DOFCST IS FALSE (BECAUSE ACARS ARE NOW Q.C.'d IN PREPOBS_ACARSQC
129 C> PROGRAM)
130 C> -2003-09-02 D. A. KEYSER -- ADDED "MSONET" TO THE LIST OF TABLE A
131 C> MESSAGE TYPES THAT WILL HAVE THE VIRTUAL TEMPERATURE CALCULATED;
132 C> DOES NOT CALL UFBINT FOR OUTPUTTING DATA IF "NLEV" (4'TH
133 C> ARGUMENT) IS ZERO (NOW CAN ONLY HAPPEN FOR GOESND FORECAST DATA
134 C> WHEN ONLY RADIANCES ARE PRESENT)
135 C> -2004-08-30 D. A. KEYSER -- NOW INCLUDES THE 4 LAYER PWATERS, THESE
136 C> GET AN OBS. ERROR (EACH THE SAME AS TOTAL PWATER) AND AN EVENT
137 C> IS GENERATED WITH A REJECTED Q.M. FOR THE 4 LAYER PWATERS IF THE
138 C> PWATER OBS. ERROR READ IN IS MISSING (THIS CHANGE ALLOWS THE ETA/
139 C> GSI TO PROCESS OBS. ERRORS IN THE PREPBUFR FILE THE SAME AS THE
140 C> ETA/3DVAR DID WHEN READING THE OBS. ERRORS FROM AN EXTERNAL
141 C> FILE); FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY OF THE REPORTED
142 C> (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS
143 C> NEW REASON CODE 3
144 C> -2004-09-10 D. T. KLEIST -- ADDED CAPABILITY TO READ GUESS FIELDS FROM
145 C> EITHER HYBRID OR, AS BEFORE, SIGMA GLOBAL FORECAST FILES
146 C> -2005-01-03 D. A. KEYSER -- FIXED ERROR READING CDAS SGES FILE WHICH
147 C> STILL HAS A 207-WORD HEADER (T62) {2004-09-10 CHANGE ASSUMED ALL
148 C> SGES FILES HAD A 226-WORD HEADER (T254), BUT THIS IS VALID ONLY
149 C> FOR GFS SGES)
150 C> -2006-05-05 R. E. TREADON -- CHANGE VERTICAL INTERPOLATION TO DIRECTLY
151 C> USE PRESSURE PROFILE, NOT PRESSURE PROFILE CONVERTED TO SIGMA.
152 C> THIS CHANGE IS IN SUBROUTINE GBLEVN03. AS A RESULT OF THIS
153 C> CHANGE, SUBROUTINE GBLEVN07 WAS REMOVED.
154 C> -2006-07-14 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "SOME_FCST"
155 C> WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOFCST" IS FALSE: IF
156 C> DOFCST=F AND SOME_FCST=T THEN, JUST AS BEFORE WHEN DOFCST=F, A
157 C> FORECAST WILL STILL BE ENCODED FOR REPORTS IN CERTAIN MESSAGE
158 C> TYPES USED IN SUBSEQUENT Q.C. PROGRAMS (I.E, "ADPUPA", "AIRCFT",
159 C> "AIRCAR", "PROFLR" OR "VADWND") (THE DEFAULT FOR SOME_FCST IS
160 C> TRUE); HOWEVER IF DOFCST=F AND SOME_FCST=F THEN A FORECAST WILL
161 C> NOT BE ENCODED INTO REPORTS IN ANY MESSAGE TYPE (THIS ALLOWS
162 C> THIS PROGRAM TO ENCODE OBS ERRORS AND/OR VIRTUAL TEMPERATURE
163 C> EVENTS INTO A PREPBUFR FILE WITHOUT ENCODING A FORECAST); ADDED
164 C> NEW NAMELIST SWITCH "ADPUPA_VIRT" WHICH, WHEN TRUE, INCLUDES
165 C> REPORTS IN MESSAGE TYPE ADPUPA (I.E., RAOBS, DROPS, MULTI-LEVEL
166 C> RECCOS) IN THE "VIRTMP" PROCESSING (PROCESSING THEM WITH SAME
167 C> LOGIC AS IN SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR)
168 C> {NORMALLY "ADPUPA_VIRT" IS FALSE (DEFAULT) BECAUSE SUBSEQUENT
169 C> PROGRAM PREPOBS_CQCBUFR PERFORMS THIS FUNCTION}
170 C> -2007-09-14 S. MOORTHI -- ADDED CAPABILITY TO READ GENERALIZED SIGMA/
171 C> HYBRID FILES FROM THE GFS USING "SIGIO" UTILITY; ALSO, CLEANED UP
172 C> SOME CODE; NEW ERROR CONDITION CODES 70 AND 71 ADDED
173 C> -2007-09-14 D. A. KEYSER -- FUNCTION OEFG01, WHICH RETURNS THE OBS
174 C> ERROR FOR A REQUESTED VARIABLE INTERP. TO A DEFINED PRESSURE
175 C> LEVEL FOR A DEFINED REPORT TYPE, MODIFIED TO USE EXACT LOGIC AS
176 C> IN GSI (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE
177 C> TYPE, LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO
178 C> MB, A FEW OTHER MINOR CHANGES) - THIS WILL ALLOW GSI TO READ OBS
179 C> ERROR DIRECTLY OUT OF PREPBUFR FILE RATHER THAN OUT OF AN
180 C> EXTERNAL FILE; FOR PW TYPES, NOW PASSES REPORTED SURFACE PRESSURE
181 C> (PRSS * 0.01) INTO FUNCTION OEFG01 RATHER THAN VERTICAL
182 C> COORDINATE PRESSURE (POB), SINCE LATTER IS ALWAYS MISSING FOR
183 C> THESE TYPES (DOESN'T CHANGE VALUE COMING OUT OF OEFG01 SINCE IT
184 C> IS CONSTANT ON ALL LEVELS ANYWAY FOR PW); IN SUBR. GBLEVN02, Q.M.
185 C> 9 IS NOW ASSIGNED TO A VARIABLE ONLY IF ITS OBS ERROR IS MISSING,
186 C> OR IN THE CASE OF MOISTURE IF THE LEVEL IS ABOVE PRESSURE LEVEL
187 C> "QTOP_REJ" OR IF ITS TEMPERATURE OBS ERROR IS MISSING, ALL OTHER
188 C> EVENT (E.G., GROSS CHECK ERRORS) ASSIGN Q.M. 8 (EVEN IF OBS ERROR
189 C> IS MISSING), PRIOR TO THIS ONLY REJECTION OF PRESSURE ON LEVEL
190 C> RESULTED IN Q.M. 8, ALL OTHER REJECTIONS GOT Q.M. 9 - THIS MEANS
191 C> TRULY "BAD" OBS WILL NOW ALWAYS GET Q.M. 8 AND ONLY OBS FLAGGED
192 C> FOR NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M.
193 C> 9 (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES
194 C> NOT EVEN CONSIDER OBS WITH Q.M. 8); CORRECTED ERROR WHICH
195 C> MISTAKENLY ASSIGNED REASON CODE OF 9 INSTEAD OF 3 TO MOISTURE
196 C> WITH MISSING OBS ERROR; IN SUBR. GBLEVN02, Q.M. 9 WILL NOT BE
197 C> ASSIGNED TO A VARIABLE IF THAT VARIABLE ALREADY HAS A "BAD" Q.M.
198 C> (I.E., > 3 BUT < 15), IN FACT THE "PREVENT" EVENT WHICH WOULD
199 C> ASSIGN Q.M. 9 IS SKIPPED ENTIRELY (DO NOT WANT THE GSI TO MONITOR
200 C> THE OBS WHICH REALLY ARE ARE "BAD"); IN SUBR. GBLEVN08, FOR NON-
201 C> "ADPUPA" TYPES, Q.M. 9 IS NOW ASSIGNED TO CALCULATED VIRT. TEMPS
202 C> IF THE MOISTURE Q.M. IS 9 OR 15 AND ORIG. TEMP NOT "BAD", THESE
203 C> "VIRTMP" EVENTS RECEIVE NEW REASON CODE 4, HAD RECEIVED Q.M. 8
204 C> WITH REASON CODE 2 LIKE VIRT. TEMPS CALCULATED FROM "BAD"
205 C> MOISTURE - THIS MEANS ONLY TRULY "BAD" VIRT. TEMPS WILL NOW GET
206 C> Q.M. 8 AND VIRT. TEMPS FLAGGED FOR NON-USE BY ASSIMILATION (BUT
207 C> STILL "GOOD") WILL NOW GET Q.M. 9 (GSI MONITORS, BUT DOES NOT
208 C> USE, OBS WITH Q.M. 9, BUT IT DOES NOT EVEN CONSIDER OBS WITH Q.M.
209 C> 8); IN SUBR. GBLEVN08, FOR "ADPUPA" TYPES, Q.M. 3 IS NOW ASSIGNED
210 C> TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS TRULY BAD
211 C> (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM IS 1 OR 2
212 C> AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN QQM WAS 9 OR
213 C> 15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES, ENCODES A
214 C> SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP"
215 C> EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO WHAT IS
216 C> ALREADY DONE FOR "RASSDA" TYPES)
217 C> -2010-01-29 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "RECALC_Q"
218 C> WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOVTMP" IS FALSE: IF
219 C> DOVTMP=F AND RECALC_Q=T THEN, JUST AS BEFORE WHEN DOVTMP=F, SPEC.
220 C> HUMIDITY IS STILL RE-CALCULATED AND THE EVENT IS ENCODED INTO THE
221 C> PREPBUFR FILE (BUT VIRTUAL TEMP. IS NOT ENCODED) (THE DEFAULT FOR
222 C> RECALC_Q IS TRUE), HOWEVER IF DOVTMP=F AND RECALC_Q=F THEN SPEC.
223 C> HUMIDITY IS NOT RE-CALCULATED (AND NEITHER IS VIRTUAL
224 C> TEMPERATURE) (THIS ALLOWS THIS PROGRAM TO BYPASS ALL "VIRTMP"
225 C> EVENT PROCESSING); ADDED NEW NAMELIST SWITCH "DOPREV" WHICH,
226 C> WHEN TRUE, WRITES "PREVENT" EVENTS INTO THE PREPBUFR FILE (IT
227 C> ALWAYS DID THIS BEFORE) (DEFAULT), BUT NOW ALLOWS THE PROGRAM
228 C> TO BYPASS "PREVENT" EVENT PROCESSING WHEN DOPREV=F; INITIALIZED
229 C> ARRAY IDATE AS ZERO IN SUBR. GBLEVN10, CORRECTED BUG WHICH
230 C> EXPOSED PREVIOUSLY HIDDEN MEMORY CLOBBERING WHEN CALLING PROGRAMS
231 C> WERE LINKED TO NEW BUFRLIB; RULES IN SUBROUTINE GBLEVN02 REFINED
232 C> TO INCLUDE FULL SFC PRESSURE SANITY CHECK FOR ALL SFC REPORTS
233 C> (MASS, 18x, & WIND, 28x), BEFORE ONLY DONE FOR SFC MASS REPORTS
234 C> (18x) AND STILL NOT DONE FOR NON-SFC WIND REPORTS SINCE LOWEST
235 C> LEVEL PRESSURE NOT NECESSARILY AT THE SFC), AS A RESULT 28x WINDS
236 C> WILL NOW GET QM=8 IF PRESSURE FAILS SANITY CHECK (OFTEN HAPPENS
237 C> IN MESONET REPORTS) (GSI WAS ALREADY NOT USING THESE WINDS SINCE
238 C> PRESSURE QM SET TO 8 ALL ALONG)
239 C> -2012-11-20 J. WOOLLEN INITIAL PORT TO WCOSS. ADDED CALL TO BUFRLIB
240 C> ROUTINE GETBMISS TO ADAPT BMISS TO LINUX ENVIRONMENT IF NEED BE
241 C> {I.E., OBTAINS BUFRLIB MISSING (BMISS) VIA CALL TO GETBMISS
242 C> RATHER THAN HARDWIRING IT TO 10E10 (10E10 CAN CAUSE INTEGER
243 C> OVERFLOW ON WCOSS - SEE CALLING PROGRAM FOR MORE INFO)}
244 C> -2013-02-13 D. A. KEYSER -- FINAL CHANGES TO RUN ON WCOSS: USE
245 C> FORMATTED PRINT STATEMENTS WHERE PREVIOUSLY UNFORMATTED PRINT WAS
246 C> > 80 CHARACTERS; RENAME ALL REAL(8) VARIABLES AS *_8
247 C> -2013-04-12 D. A. KEYSER -- IN SUBROUTINE GBLEVN08, DON'T ALLOW
248 C> CALCULATED Q TO BE < 0 WHICH CAN OCCUR ON WCOSS FOR CASES OF
249 C> HORRIBLY BAD MESONET DATA
250 C> -2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH,
251 C> WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION
252 C> PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST-
253 C> COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN
254 C> REPORTED PMSL IS MISSING (DONE IN SUBROUTINE GBLEVN08). DOVTMP
255 C> MUST BE TRUE AND DOANLS MUST BE FALSE ("PREVENTS" MODE). THE
256 C> DERIVED PMSL EITHER GETS A QUALITY MARK ("PMQ") OF 3 OR INHERITS
257 C> THE STATION PRESSURE QUALITY MARK ("PQM"), WHICHEVER IS GREATER.
258 C> THE VALUE OF THE PMSL INDICATOR (NEW MNEMONIC "PMIN") IS SET TO 1
259 C> TO DENOTE PMSL WAS DERIVED RATHER THAN OBSERVED. THE DEFAULT FOR
260 C> "DOPMSL" IS FALSE (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS). IT
261 C> IS FORCED TO BE FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE). IN
262 C> SUBROUTINE GBLEVN02, SFCSHP REPORTS WITH CALM WINDS AND NON-
263 C> MISSING BACKGROUND U- OR V-COMPONENT WIND .GE. 5 M/SEC ARE
264 C> FLAGGED WITH Q.M. 8 (EVENT PGM "PREVENT", REASON CODE 8).
265 C> -2014-05-08 JWhiting -- altered print statement (2 format) in GBLEVN10
266 C> subroutine; increased field width for spectral resolution to
267 C> accommodate models w/ up to 5-digit resolution (I3 to I5).
268 C> -2016-06-13 FANGLIN YANG AND RUSS TREADON -- HANG LEI ADDED NEMSIO
269 C> TO SUBROUTINE GBLEVN10 AND REMOVED ALL SIGIO CAPABILITY.
270 C> THIS UPDATE RESTORES GBLEVN10 FOR PROCESSING SIGIO INPUT, AND
271 C> ADDS A NEW GBLEVN12 FOR PROCESSING NEMSIO INPUT. THE INPUT GFS
272 C> FILE TYPE SIGIO VS NEMSIO IS NOW DETERMINED IN THE MAIN PROGRAM.
273 C> THE CODE IS ALSO UPDATED TO REMOVE BUGS. SUBROUTINE SIGIO_MODPR
274 C> IS USED TO COMPUTE LAYER AND INTERFACE PRESSURE FOR NEMSIO INPUT.
275 C> -2017-02-17 D Keyser & J Whiting -- In subroutine GBLEVN12, removed
276 C> references to multiple input files, since only 1 nemsio formatted
277 C> input is needed (no interpolation is attempted; c.f.; what is
278 C> done with sigio formatted inputs); GBLEVN12 routine still retains
279 c some structures (esp array references) left over from multi file
280 C> input. Updated comments and docblock to account for new NEMSIO
281 C> input.
282 C> -2017-02-22 D. Keyser -- Further changes to subr. GBLEVN12 to remove
283 C> array and logic references to multiple input files.
284 C> -2019-10-31 Hang Lei -- Add GBLEVN13 to process netcdf input.
285 C>
286 C>
287 C> USAGE: CALL GBLEVENTS(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET,
288 C> $ NEWTYP)
289 C>
290 C> INPUT ARGUMENT LIST:
291 C> @param IDATEP - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
292 C> @param IUNITF - 2-WORD ARRAY:
293 C> For SIGIO input:
294 C> - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
295 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
296 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
297 C> THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
298 C> IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
299 C> AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
300 C> IS A MULTIPLE OF 3
301 C> - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
302 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
303 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
304 C> THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
305 C> OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
306 C> THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
307 C> For NEMSIO input:
308 C> - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
309 C> (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
310 C> THE DATE IN IDATEP
311 C> - WORD 2 - NOT USED, SHOULD BE A NULL FILE
312 C> @param IUNITE - UNIT NUMBER OF INPUT OBSERVATION ERROR FILE
313 C> - (USED ONLY IN PREVENTS MODE)
314 C> @param IUNITP - UNIT NUMBER OF OUTPUT PREPBUFR DATA SET
315 C> @param IUNITS - UNIT NUMBER OF "PREVENT" EVENTS DATA FILTERING
316 C> - SUMMARY PRINT FILE
317 C> - (USED ONLY IN PREVENTS MODE)
318 C> @param SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
319 C> - REPORT BEING PROCESSED
320 C> @param NEWTYP - INDICATOR IF THE BUFR MESSAGE TABLE A ENTRY HAS
321 C> - CHANGED FROM THAT OF THE PREVIOUS REPORT (=0 - NO,
322 C> - =1 - YES)
323 C>
324 C>
325 C> INPUT FILES:
326 C> - UNIT 05 - STANDARD INPUT (DATA CARDS - SEE NAMELIST
327 C> DOCUMENTATION BELOW)
328 C> (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS
329 C> SUBROUTINE RUNS IN POSTEVENTS MODE)
330 C> - UNIT AA - PREPBUFR DATA SET
331 C> (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN
332 C> INPUT ARGUMENT LIST)
333 C> - UNIT BB - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE)
334 C> FILE
335 C> (WHERE BB IS UNIT NUMBER DEFINED AS IUNITF(1) IN
336 C> INPUT ARGUMENT LIST)
337 C> - UNIT CC - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE)
338 C> FILE
339 C> (WHERE CC IS UNIT NUMBER DEFINED AS IUNITF(2) IN
340 C> INPUT ARGUMENT LIST)
341 C> NOTE: only valid for SIGIO input
342 C> - UNIT DD - OBSERVATION ERROR FILE (WHERE DD IS UNIT NUMBER
343 C> DEFINED AS IUNITE IN INPUT ARGUMENT LIST)
344 C> (USED ONLY IN PREVENTS MODE)
345 C>
346 C> OUTPUT FILES:
347 C> - UNIT 06 - STANDARD OUTPUT PRINT
348 C> - UNIT AA - PREPBUFR DATA SET
349 C> (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN
350 C> INPUT ARGUMENT LIST)
351 C> - UNIT DD - "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE
352 C> (WHERE DD IS UNIT NUMBER DEFINED AS IUNITS IN
353 C> INPUT ARGUMENT LIST)
354 C> (USED ONLY IN PREVENTS MODE)
355 C>
356 C> SUBPROGRAMS CALLED:
357 C> UNIQUE: GBLEVN02 GBLEVN03 GBLEVN04 GBLEVN06 OEFG01
358 C> GBLEVN08 GBLEVN10 GBLEVN11 GBLEVN11D GBLEVN12
359 C> GBLEVN13
360 C> GETLATS
361 C> MODULES: GBLEVN_MODULE SIGIO_MODULE SIGIO_R_MODULE
362 C> NEMSIO_MODULE NEMSIO_OPENCLOSE NEMSIO_READ
363 C> NEMSIO_WRITE
364 C> LIBRARY:
365 C> SIGIO - SIGIO_RROPEN SIGIO_RRHEAD SIGIO_SCLOSE SIGIO_ALDATS
366 C> SIGIO_ALDATM SIGIO_RRDATS SIGIO_RRDATM SIGIO_AXDATS
367 C> SIGIO_AXDATM SIGIO_MODPR SIGIO_CNVTDV
368 C> SPLIB - SPTEZM SPTEZMV SPLAT
369 C> W3NCO - W3MOVDAT ERREXIT
370 C> BUFRLIB - UFBINT UFBQCD GETBMISS IBFMS
371 C> NEMSIO - NEMSIO_OPEN NEMSIO_CLOSE NEMSIO_INIT
372 C> NEMSIO_GETFILEHEAD NEMSIO_READRECV NEMSIO_FINALIZE
373 C> NEMSIO_GETHEADVAR NEMSIO_READRECVw34
374 C>
375 C> EXIT STATES:
376 C> - COND = 0 - SUCCESSFUL RUN
377 C> - COND = 60 - OBSERVATION ERROR TABLE EMPTY OR DOES NOT EXIST
378 C> - COND = 61 - VARIABLE NLTD .NE. VARIABLE NLEV
379 C> - COND = 62 - VARIABLE NLTQ .NE. VARIABLE NLEV
380 C> - COND = 63 - VARIABLE NLQQ .NE. VARIABLE NLEV
381 C> - COND = 68 - DATE OF FIRST GUESS/ANALYSIS FILE(S) DOES NOT MATCH,
382 C> OR AT LEAST SPAN, THE CENTER DATE FOR THE PREPBUFR
383 C> FILE
384 C> - COND = 69 - FOR SIGIO INPUT GLOBAL FILES, VARIABLE KMAX TOO BIG
385 C> - UNABLE TO TRANSFORM FIRST GUESS OR ANALYSIS FILE(S)
386 C> - COND = 70 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RROPEN
387 C> RETURNED WITH NON-ZERO R.C.
388 C> - COND = 71 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RRHEAD
389 C> RETURNED WITH NON-ZERO R.C.
390 C>
391 C>
392 C> REMARKS: THIS SUBROUTINE MAY NOT WORK CORRECTLY IN THE EIGHT BYTE
393 C> INTEGER W3NCO (_8) LIBRARY. PLEASE COMPILE APPLICATION CODE
394 C> USING A FOUR BYTE REAL W3NCO LIBRARY (_4 OR _d).
395 C>
396 C> THIS ROUTINE PROCESSES ONE REPORT AT A TIME. IT EXPECTS THAT THE
397 C> CALLING PROGRAM HAS ALREADY ENCODED THE REPORT INTO THE PREPBUFR
398 C> FILE VIA THE UFBINT OR UFBCPY ROUTINES. THE CALLING PROGRAM
399 C> SHOULD THEN CALL THIS ROUTINE AND, UPON ITS RETURN, THE CALLING
400 C> PROGRAM SHOULD CALL WRITSB TO ACTUALLY WRITE THE UPDATED SUBSET
401 C> (REPORT) INTO THE BUFR MESSAGE.
402 C>C
403 C> ***** VARIABLES IN NAMELIST PREVDATA READ IN BY THIS SUBROUTINE *****
404 C> (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS
405 C> SUBROUTINE RUNS IN POSTEVENTS MODE - DOANLS=TRUE
406 C> AND ALL OTHER VARIABLES ARE SET TO FALSE)
407 C>C
408 C>C
409 C> - DOPREV - WRITE "PREVENT" EVENT INTO THE PREPBUFR FILE?
410 C> DOPREV = .TRUE. ---> YES (DEFAULT)
411 C> DOPREV = .FALSE. ---> NO
412 C> - DOVTMP, ADPUPA_VIRT & RECALC_Q:
413 C> DOVTMP - WRITE VIRTUAL TEMPERATURE EVENT ("VIRTMP") INTO THE
414 C> PREPBUFR FILE (I.E., RE-CALCULATE SPECIFIC HUMIDITY
415 C> THEN CALCULATE VIRTUAL TEMPERATURE) FOR THE FOLLOWING
416 C> TYPES OF REPORTS:
417 C> ADPUPA_VIRT = .FALSE. ---> SURFACE LAND, MARINE,
418 C> MESONET AND RASS REPORTS?
419 C> ADPUPA_VIRT = .TRUE. ---> SURFACE LAND, MARINE,
420 C> MESONET RASS, RAOB, DROP
421 C> AND MULTI-LEVEL RECCO
422 C> REPORTS?
423 C> FOR ALL TYPES EXCEPT RASS, THIS WILL ATTEMPT TO
424 C> CALCULATE VIRTUAL TEMPERATURE FROM SENSIBLE TEMPERATURE
425 C> AND THE JUST RE-CALCULATED SPECIFIC HUMIDITY AND ENCODE
426 C> IT AS A STACKED EVENT IN THE PREPBUFR FILE. FOR RASS
427 C> REPORTS THIS WILL JUST ENCODE THE REPORTED TEMPERATURE
428 C> AS A STACKED EVENT IN THE PREPBUFR FILE SINCE THE
429 C> REPORTED TEMPERATURE IS ALREADY VIRTUAL (NO MOISTURE IS
430 C> PRESENT SO Q IS NOT RE-CALCULATED FOR RASS REPORTS).
431 C> DOVTMP = .TRUE. ---> YES (DEFAULT)
432 C> DOVTMP = .FALSE.
433 C> RECALC_Q = .TRUE. ---> RE-CALCULATE SPECIFIC
434 C> HUMIDITY BUT DO NOT THEN
435 C> CALCULATE VIRTUAL
436 C> TEMPERATURE (DEFAULT)
437 C> RECALC_Q = .FALSE. ---> NO, DO NOT RE-CALCULATE
438 C> SPECIFIC HUMIDITY AND DO
439 C> NOT CALCULATE VIRTUAL
440 C> TEMPERATURE
441 C> {NOTE1: FOR SURFACE LAND, MARINE AND MESONET REPORTS, (AND
442 C> RAOB, DROP AND MULTI-LEVEL RECCO REPORTS IF
443 C> "ADPUPA_VIRT"=TRUE) DOVTMP=FALSE WILL STILL RE-CALCULATE
444 C> SPECIFIC HUMIDITY AND ENCODE IT AS A STACKED EVENT IN
445 C> THE PREPBUFR FILE UNLESS EITHER DOANLS IS TRUE OR
446 C> RECALC_Q IS FALSE.)
447 C> (NOTE2: DOES NOT APPLY TO ANY REPORT TYPES OTHER THAN THOSE
448 C> MENTIONED ABOVE)
449 C> (NOTE3: IF DOANLS=TRUE, THEN DOVTMP IS NOT ONLY FORCED TO BE
450 C> FALSE, BUT ALSO SPECIFIC HUMIDITY IS NOT RE-CALCULATED.)
451 C> (NOTE4: ADPUPA_VIRT DEFAULTS TO FALSE.)
452 C> (NOTE5: IF DOVTMP=TRUE, THEN RECALC_Q IS MEANINGLESS.)
453 C> (NOTE6: RECALC_Q DEFAULTS TO TRUE.)
454 C>
455 C> - DOFCST & SOME_FCST:
456 C> DOFCST - ENCODE FORECAST (FIRST GUESS) VALUES, INTERPOLATED FROM
457 C> EITHER A SIGIO (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT
458 C> GLOBAL FILE, INTO THE PREPBUFR FILE FOR ALL MESSAGE
459 C> TYPES OR AT LEAST SOME MESSAGE TYPES?
460 C> DOFCST = .TRUE. ---> YES, ENCODE FORECST FOR ALL
461 C> MESSAGE TYPES (DEFAULT)
462 C> DOFCST = .FALSE.
463 C> SOME_FCST = .FALSE. ---> NO, DO NOT ENCODE FORECAST
464 C> FOR ANY MESSAGE TYPE
465 C> (VALUES REMAIN MISSING)
466 C> SOME_FCST = .TRUE. ---> YES, BUT ONLY FOR MESSAGE
467 C> TYPES "ADPUPA", "AIRCFT",
468 C> "AIRCAR", "PROFLR" OR
469 C> "VADWND" (VALUES REMAIN
470 C> MISSING FOR ALL OTHER
471 C> MESSAGE TYPES)
472 C> (NOTE1: THE CASE DOFCST=FALSE & SOME_FCST=TRUE WRITES THE
473 C> FORECAST VALUES FOR THE TYPES MENTIONED ABOVE BECAUSE
474 C> THEY ARE NEEDED BY SUBSEQUENT QUALITY CONTROL PROGRAMS.)
475 C> (NOTE2: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE
476 C> NON-GLOBAL VERSIONS SINCE ONLY THE GLOBAL REQUIRES A
477 C> FIRST GUESS TO BE PRESENT FOR ALL CONVENTIONAL MESSAGE
478 C> TYPES IN THE PREPBUFR FILE.)
479 C> (NOTE3: IF DOANLS=TRUE, THEN DOFCST & SOME_FCST ARE FORCED TO BE
480 C> FALSE, MEANING A GUESS WILL NOT BE ENCODED FOR ANY
481 C> MESSAGE TYPE.)
482 C> (NOTE4: IF DOFCST=TRUE, THEN SOME_FCST IS MEANINGLESS.)
483 C> (NOTE5: SOME_FCST DEFAULTS TO TRUE.)
484 C>
485 C> - DOANLS - ENCODE ANALYZED VALUES, INTERPOLATED FROM EITHER A SIGIO
486 C> (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT GLOBAL FILE,
487 C> INTO THE PREPBUFR FILE -
488 C> POSTEVENTS MODE - ?
489 C> DOANLS = .TRUE. ---> YES, FOR ALL MESSAGE TYPES
490 C> DOANLS = .FALSE. ---> NO, FOR ALL MESSAGE TYPES
491 C> - PREVENTS MODE - (DEFAULT)
492 C> (NOTE: DOANLS=TRUE WILL OVERRIDE AND FORCE TO FALSE ALL OTHER
493 C> SWITCHES. IN ADDITION, THE FORECAST VALUES WILL NOT
494 C> BE ENCODED FOR ANY MESSAGE TYPE AND SPECIFIC HUMIDITY
495 C> WILL NOT BE RE-CALCULATED.)
496 C>
497 C> - DOBERR - ENCODE OBSERVATIONAL ERROR VALUES, AS READ FROM OBS.
498 C> ERROR FILE, INTO THE PREPBUFR FILE?
499 C> DOBERR = .TRUE. ---> YES (DEFAULT)
500 C> DOBERR = .FALSE. ---> NO (VALUES REMAIN MISSING)
501 C> (NOTE1: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE
502 C> RAP -AND PREVIOUS RUC- VERSION SINCE IT DOES NOT REQUIRE
503 C> OBSERVATIONAL ERRORS TO BE PRESENT IN THE PREPBUFR FILE.)
504 C> (NOTE2: IF DOANLS=TRUE, THEN DOBERR IS FORCED TO BE FALSE.)
505 C>
506 C> - QTOP_REJ - THE PRESSURE LEVEL (IN MB) ABOVE WHICH ALL SPECIFIC
507 C> HUMIDITY QUALITY MARKERS ARE "REJECTED" (THE QUALITY
508 C> MARKER IS SET TO 9 ON ALL PRESSURE LEVELS LESS THAN
509 C> THIS LEVEL) (DEFAULT=300.)
510 C>
511 C> - SATMQC - PERFORM SPECIAL QUALITY CONTROL ON SATELLITE TEMPERATURE
512 C> SOUNDINGS IN REPORT TYPES 160-179?
513 C> SATMQC = .TRUE. ---> YES
514 C> SATMQC = .FALSE. ---> NO (DEFAULT)
515 C> (NOTE: THIS APPLIES ONLY TO THE CDAS OR HISTORICAL RE-RUNS
516 C> WITH TEMPERATURE SOUNDINGS IN THESE REPORT TYPES)
517 C>
518 C> - DOPMSL - ENCODE DERIVED PMSL ("PMO") FOR ALL SURFACE REPORTS WHEN
519 C> REPORTED PMSL IS MISSING - ?
520 C> DOPMSL = .TRUE. ---> YES
521 C> DOPMSL = .FALSE. ---> NO ("PMO" REMAINS MISSING)(DEFAULT)
522 C> {NOTE: THIS APPLIES ONLY WHEN DOVTMP=TRUE AND DOANLS=FALSE
523 C> ("PREVENTS" MODE), VIRTUAL TEMPERATURE CAN BE CALCULATED,
524 C> AND STATION PRESSURE AND SURFACE HEIGHT/ELEVATION ARE
525 C> BOTH PRESENT. THE DERIVED PMSL EITHER GETS A QUALITY
526 C> MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE
527 C> QUALITY MARK ("PQM") PQM, WHICHEVER IS GREATER. THE VALUE
528 C> OF THE PMSL INDICATOR ("PMIN") IS SET TO 1 TO DENOTE PMSL
529 C> WAS DERIVED RATHER THAN OBSERVED.}
530 C>
531  SUBROUTINE gblevents(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET,
532  $ NEWTYP)
533 
534  USE sigio_module
535  USE sigio_r_module
536  USE nemsio_module
537  USE nemsio_openclose
538  USE nemsio_read
539  USE gblevn_module
540  use netcdf
541 
542  INTEGER, PARAMETER :: IM=384, jm=im/2+1
543  integer :: idrt = 0
544 
545  CHARACTER*80 HEADR,OBSTR,QMSTR,FCSTR,OESTR,ANSTR
546  CHARACTER*8 SUBSET
547  REAL(8) OBS_8,QMS_8,BAK_8,SID_8,HDR_8(10)
548  REAL(8) BMISS,GETBMISS
549  LOGICAL DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT,DOANLS,
550  $ satmqc,adpupa_virt,recalc_q,doprev,dopmsl
551  INTEGER*4 IRET,IRET1
552  INTEGER*4 IUNITF(2),ncid
553 
554  CHARACTER*20 CFILE1
555  TYPE(SIGIO_HEAD) :: HEAD1
556  TYPE(NEMSIO_GFILE) :: GFILE1
557 
558  COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
559  $ xob,yob,dhr,typ,nlev
560  COMMON /gbevbb/ pvcd,vtcd
561  COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
562  $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
563  COMMON /gbevdd/ errs(300,33,6)
564  COMMON /gbevff/ bmiss
565 
566  SAVE
567 
568  DATA ifirst/0/
569 
570  DATA headr /
571  $ 'SID XOB YOB DHR TYP '/
572  DATA obstr /
573  $ 'POB QOB TOB ZOB UOB VOB PWO PW1O PW2O PW3O PW4O CAT PRSS '/
574  DATA qmstr /
575  $ 'PQM QQM TQM ZQM WQM PWQ PW1Q PW2Q PW3Q PW4Q NUL NUL '/
576  DATA fcstr /
577  $ 'PFC QFC TFC ZFC UFC VFC PWF PW1F PW2F PW3F PW4F NUL '/
578  DATA anstr /
579  $ 'PAN QAN TAN ZAN UAN VAN PWA PW1A PW2A PW3A PW4A NUL '/
580  DATA oestr /
581  $ 'POE QOE TOE ZOE WOE PWE PW1E PW2E PW3E PW4E NUL NUL '/
582 
583  namelist /prevdata/dovtmp,dofcst,some_fcst,doberr,doanls,
584  $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
585 
586 C----------------------------------------------------------------------
587 C----------------------------------------------------------------------
588 
589  IF(ifirst.EQ.0) THEN
590 
591 C -------------------------------
592 C FIRST TIME IN DO A FEW THINGS...
593 C -------------------------------
594 
595  ifirst = 1
596  print 700
597  700 FORMAT(/1x,100('#')/' =====> SUBROUTINE GBLEVENTS INVOKED FOR ',
598  $ 'THE FIRST TIME - VERSION LAST UPDATED 2017-02-22'/)
599 
600  bmiss = getbmiss()
601  print *
602  print *, 'BUFRLIB value for missing passed into GBLEVENTS ',
603  $ 'is: ',bmiss
604  print *
605 
606 C INITIALIZE NAMELIST SWITCHES TO DEFAULT VALUES
607 C ----------------------------------------------
608 
609  dovtmp = .true.
610  doprev = .true.
611  recalc_q = .true.
612  dofcst = .true.
613  some_fcst = .true.
614  doberr = .true.
615  doanls = .false.
616  qtop_rej = 300.
617  satmqc = .false.
618  adpupa_virt = .false.
619  dopmsl = .false.
620  READ(5,prevdata,err=101,end=102)
621  GO TO 103
622 C-----------------------------------------------------------------------
623  101 CONTINUE
624 
625 C ERROR READING STANDARD INPUT - THIS DEFAULTS TO POSTEVENTS MODE
626 C ---------------------------------------------------------------
627 
628  print 7013
629  7013 FORMAT(/' ##> GBLEVENTS: ERROR READING STANDARD INPUT DATA CARDS',
630  $ ' -- DEFAULTS TO "POSTEVENTS" MODE'/)
631  doanls = .true.
632  GO TO 103
633 
634 C-----------------------------------------------------------------------
635  102 CONTINUE
636 
637 C STANDARD INPUT IS EMPTY - THIS DEFAULTS TO POSTEVENTS MODE
638 C ----------------------------------------------------------
639 
640  print 7014
641  7014 FORMAT(/' ##> GBLEVENTS: STANDARD INPUT DATA CARDS DO NOT ',
642  $ 'EXIST -- DEFAULTS TO "POSTEVENTS" MODE'/)
643  doanls = .true.
644 
645 C-----------------------------------------------------------------------
646  103 CONTINUE
647  IF(doanls) THEN
648  dovtmp = .false.
649  doprev = .false.
650  dofcst = .false.
651  some_fcst = .false.
652  doberr = .false.
653  adpupa_virt = .false.
654  dopmsl = .false.
655  ENDIF
656  IF(dovtmp) recalc_q=.true. ! RECALC_Q must be T if DOVTMP is T
657  WRITE (6,prevdata)
658 
659  fcst = dofcst
660  virt = .false.
661 
662 C CHECK VALID-TIME DATE OF GUESS/ANALYSIS FILE(S) AGAINST THE CENTER
663 C DATE FOR THE PREPBUFR FILE AND OBTAIN THE FIRST GUESS/ANALYSIS
664 C UNLESS ALL OF DOFCST, SOME_FCST, DOANLS ARE FALSE
665 C ------------------------------------------------------------------
666 
667  IF(.NOT.doanls) THEN
668  IF(.NOT.dofcst.AND..NOT.some_fcst) THEN
669  print 901
670  901 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - FIRST GUESS NOT READ ',
671  $ 'IN'/)
672  ELSE
673  print 701
674  701 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - DATE CHECK AND ',
675  $ 'TRANSFORM THE FIRST GUESS'/)
676  ENDIF
677  ELSE
678  print 7701
679  7701 FORMAT(/' --> GBLEVENTS: POSTEVENTS MODE - DATE CHECK AND ',
680  $ 'TRANSFORM THE ANALYSIS'/)
681  ENDIF
682 
683  IF(dofcst .OR. some_fcst .OR. doanls) THEN
684  WRITE(cfile1,'("fort.",I2.2)') iunitf(1)
685 
686  iret = nf90_open(trim(cfile1),nf90_nowrite,ncid)
687  if (iret == 0) then
688  print *,' ===> GFS FCST/ANAL INPUT IS NETCDF'
689  CALL gblevn13(iunitf,idatep,im,jm,idrt)
690  else
691  CALL sigio_rropen(iunitf(1),cfile1,iret)
692  CALL sigio_srhead(iunitf(1),head1,iret1)
693  IF(iret == 0 .AND. iret1 == 0) THEN
694  print *,' ===> GLOBAL FCST/ANAL INPUT IS SIGIO'
695  CALL sigio_sclose(iunitf(1),iret)
696  CALL gblevn10(iunitf,idatep,im,jm,idrt)
697  ELSE
698  CALL nemsio_open(gfile1,trim(cfile1),'read',iret=iret)
699  IF(iret == 0) THEN
700  CALL nemsio_close(gfile1,iret=iret)
701  print *,' ===> GFS FCST/ANAL INPUT IS NEMSIO'
702  CALL gblevn12(iunitf,idatep,im,jm,idrt)
703  ENDIF
704  ENDIF
705  ENDIF
706  ENDIF
707 
708  print*,'after returning from GBLEVN10, GBLEVN12,',
709  $ ' or GBLEVN13. idrt=',idrt
710 
711  IF(doberr) THEN
712 
713 C IF REQUESTED, READ ERROR FILES (ONLY POSSIBLE IN PREVENTS MODE)
714 C ---------------------------------------------------------------
715 
716  print 702
717  702 FORMAT(/' --> GBLEVENTS: READ ERROR FILES'/)
718 
719  CALL gblevn01(iunite)
720 
721  ELSE
722 
723  errs = 0
724  IF(.NOT.doanls) print 3702
725  3702 FORMAT(/' --> GBLEVENTS: OBS. ERROR NOT ENCODED IN PREPBUFR ',
726  $ '(BY CHOICE)'/)
727 
728  ENDIF
729 
730 C OBTAIN NECESSARY PROGRAM CODES (ONLY USED IN PREVENTS MODE)
731 C -----------------------------------------------------------
732 
733  CALL ufbqcd(iunitp,'PREVENT',pvcd)
734  CALL ufbqcd(iunitp,'VIRTMP ',vtcd)
735 
736  print 703
737  703 FORMAT(/1x,100('#')/)
738 
739 C SET-UP OUTPUT "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE
740 C (ONLY USED IN PREVENTS MODE)
741 C ----------------------------------------------------------------
742 
743  IF(.NOT.doanls) WRITE(iunits,1701) idatep
744  1701 FORMAT(//130('#')//38x,'*** "PREVENT" EVENTS DATA FILTERING ',
745  $ 'SUMMARY ***'/35x,'--> CENTER DATE FOR PREPBUFR FILE IS: ',i10,
746  $ ' <--'//)
747 
748 C----------------------------------------------------------------------
749 C----------------------------------------------------------------------
750 
751  ENDIF
752 
753  IF(.NOT.doanls) THEN
754 
755  IF(newtyp.EQ.1) WRITE(iunits,1702) subset
756  1702 FORMAT(130('-')/39x,'--> SUMMARY FOR TABLE A ENTRY "',a8,'" <--'/)
757 
758  IF(.NOT.dofcst .AND. some_fcst) fcst = (subset.EQ.'ADPUPA '
759  $ .OR.subset.EQ.'PROFLR '.OR.subset .EQ.'AIRCFT '.OR.subset
760  $ .EQ.'AIRCAR '.OR.subset .EQ.'VADWND ')
761 
762 C Will not subject ACARS reports to virtual temp. processing until
763 C spec. humidity is used in production
764 
765 ccccc VIRT = (SUBSET.EQ.'ADPSFC '.OR.SUBSET.EQ.'SFCSHP '.OR.
766 ccccc$ SUBSET.EQ.'MSONET '.OR.SUBSET.EQ.'AIRCAR '.OR.
767 ccccc$ SUBSET.EQ.'RASSDA '.OR.SUBSET.EQ.'SATEMP '.OR.
768 ccccc$ (SUBSET.EQ.'ADPUPA '.AND.ADPUPA_VIRT))
769  virt = (recalc_q.AND.(subset.EQ.'ADPSFC '.OR.
770  $ subset.EQ.'SFCSHP '.OR.
771  $ subset.EQ.'MSONET '.OR.
772  $ subset.EQ.'RASSDA '.OR.
773  $ subset.EQ.'SATEMP '.OR.
774  $ (subset.EQ.'ADPUPA '.AND.adpupa_virt)))
775 
776 
777  IF(.NOT.(fcst.OR.doberr.OR.virt.OR.doprev)) THEN
778  IF(newtyp.EQ.1) WRITE(iunits,1703)
779  1703 FORMAT(/' ==> DATA FILTERING NOT PERFORMED FOR THIS TABLE A ',
780  $ 'ENTRY -- FORECAST, OBS ERROR, "VIRTMP", "PREVENT" PROCESSING ',
781  $ 'NOT DONE'/)
782  RETURN
783  ENDIF
784 
785  ENDIF
786 
787 C READY TO RETRIEVE NECESSARY INFORMATION OUT OF THE NEXT REPORT WHICH
788 C HAS BEEN "UFB" ENCODED INTO THE PREPBUFR FILE BY THE CALLING PROGRAM
789 C (USE NEGATIVE UNIT NUMBER HERE SINCE FILE OPEN FOR OUTPUT)
790 C (NOTE: THE CALLING PROGRAM HAS NOT YET WRIITEN THE REPORT INTO
791 C THE PREPBUFR FILE VIA WRITSB!)
792 C ----------------------------------------------------------------
793 
794  CALL ufbint(-iunitp,obs_8,13,255,nlev,obstr)
795  CALL ufbint(-iunitp,qms_8,12,255,nlev,qmstr)
796  CALL ufbint(-iunitp,hdr_8,10, 1,iret,headr)
797  sid_8 = hdr_8(1)
798  xob = hdr_8(2)
799  yob = hdr_8(3)
800  dhr = hdr_8(4)
801  typ = hdr_8(5)
802 
803  IF(fcst.OR.doanls) THEN
804 
805 C PREVENTS MODE: ENCODE FIRST GUESS VALUES INTO PREPBUFR FILE
806 C ------------------------------------------------------------
807 
808 C POSTEVENTS MODE: ENCODE ANALYSIS VALUES INTO REPORT AND RETURN TO
809 C CALLING PROGRAM TO WRITE GBL-EVENTED REPORT
810 C (SUBSET) INTO PREPBUFR FILE
811 C -----------------------------------------------------------------
812 
813  CALL gblevn03(subset)
814  IF(nlev.GT.0) THEN
815  IF(fcst) THEN
816  CALL ufbint(iunitp,bak_8,12,nlev,iret,fcstr)
817  ELSE
818  CALL ufbint(iunitp,bak_8,12,nlev,iret,anstr)
819  RETURN
820  ENDIF
821  ENDIF
822  ENDIF
823 
824 C --------------------------------------------------------------------
825 C LOGIC FROM HERE ON PERTAINS ONLY TO PREVENTS MODE OF THIS SUBROUTINE
826 C --------------------------------------------------------------------
827 
828 C ENCODE OBSERVATION ERRORS INTO REPORT
829 C -------------------------------------
830 
831  IF(doberr) THEN
832  IF(newtyp.EQ.1) WRITE(iunits,1710)
833  1710 FORMAT(/' ==> OBS ERROR VALUES ARE ENCODED FOR THIS TABLE A ',
834  $ 'ENTRY'//' ==> FILTERING VIA MISSING OBS ERROR TEST IS ',
835  $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES ARE ',
836  $ 'PROCESSED/STORED'/)
837  CALL gblevn04
838  IF(nlev.GT.0) CALL ufbint(iunitp,bak_8,12,nlev,iret,oestr)
839  ELSE
840  IF(newtyp.EQ.1) WRITE(iunits,1705)
841  1705 FORMAT(/' ==> OBS ERROR VALUES NOT ENCODED FOR THIS TABLE A ',
842  $ 'ENTRY'//' ==> FILTERING VIA MISSING OBS ERROR TEST NOT ',
843  $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES NOT ',
844  $ 'PROCESSED/STORED'/)
845  ENDIF
846 
847 C MAKE THE GBLEVENTS EVENTS AND ENCODE INTO REPORT
848 C ------------------------------------------------
849 
850  IF(.NOT.fcst) THEN
851  IF(newtyp.EQ.1) WRITE(iunits,1704)
852  1704 FORMAT(/' ==> FORECAST VALUES NOT ENCODED FOR THIS TABLE A ',
853  $ 'ENTRY'//' ==> FILTERING VIA POB VS. GESS PSFC TEST NOT ',
854  $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES NOT ',
855  $ 'PROCESSED/STORED'/)
856  ELSE
857  IF(newtyp.EQ.1) WRITE(iunits,1708)
858  1708 FORMAT(/' ==> FORECAST VALUES ARE ENCODED FOR THIS TABLE A ',
859  $ 'ENTRY'//' ==> FILTERING VIA POB VS. GESS PSFC TEST IS ',
860  $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES ARE ',
861  $ 'PROCESSED/STORED'/)
862  ENDIF
863 
864  IF(doprev) THEN
865  IF(newtyp.EQ.1) WRITE(iunits,1807)
866  1807 FORMAT(/' ==> "PREVENT" EVENT PROCESSING IS PERFORMED FOR THIS',
867  $ ' TABLE A ENTRY'/)
868  CALL gblevn02(iunitp,iunits,newtyp,subset)
869  ELSE
870  IF(newtyp.EQ.1) WRITE(iunits,1806)
871  1806 FORMAT(/' ==> "PREVENT" EVENT PROCESSING NOT PERFORMED FOR THIS',
872  $ ' TABLE A ENTRY'/)
873  ENDIF
874 
875 C MAKE THE VIRTUAL TEMPERATURE EVENTS AND ENCODE INTO REPORT
876 C ----------------------------------------------------------
877 
878  IF(.NOT.virt) THEN
879  IF(newtyp.EQ.1) WRITE(iunits,1706)
880  1706 FORMAT(/' ==> "VIRTMP" EVENT PROCESSING NOT PERFORMED FOR THIS ',
881  $ 'TABLE A ENTRY'/)
882  ELSE
883  IF(newtyp.EQ.1) WRITE(iunits,1707)
884  1707 FORMAT(/' ==> "VIRTMP" EVENT PROCESSING IS PERFORMED FOR THIS ',
885  $ 'TABLE A ENTRY'/)
886  CALL gblevn08(iunitp,subset)
887  ENDIF
888 
889 C RETURN TO CALLING PROGRAM TO WRITE GBL-EVENTED REPORT (SUBSET) INTO
890 C PREPBUFR FILE
891 C -------------------------------------------------------------------
892 
893  RETURN
894 
895  END
896 C***********************************************************************
897 C***********************************************************************
898 C> Read observation error table
899 C> @param IUNITE - UNIT NUMBER OF INPUT OBSERVATION ERROR FILE (USED ONLY IN PREVENTS MODE)
900  SUBROUTINE gblevn01(IUNITE) ! FORMERLY SUBROUTINE ETABLE
901 
902  COMMON /gbevdd/ errs(300,33,6)
903 
904 C READ THE OBSERVATION ERROR TABLES
905 C ---------------------------------
906 
907  rewind iunite
908 
909  irec = 0
910 
911  10 CONTINUE
912  READ(iunite,'(1X,I3)',end=100) kx
913  irec = irec + 1
914  DO k=1,33
915  READ(iunite,'(1X,6E12.5)') (errs(kx,k,m),m=1,6)
916  ENDDO
917  GO TO 10
918 
919  100 CONTINUE
920  IF(irec.LE.0) THEN
921  print'(" ##GBLEVENTS/GBLEVN01 - OBS. ERROR TABLE EMPTY OR ",
922  $ "DOES NOT EXIST - STOP 60")'
923  CALL errexit(60)
924  ENDIF
925 
926  RETURN
927 
928  END
929 C***********************************************************************
930 C***********************************************************************
931 C> Filter data
932 C> @param IUNITP - UNIT NUMBER OF OUTPUT PREPBUFR DATA SET
933 C> @param IUNITS - UNIT NUMBER OF "PREVENT" EVENTS DATA FILTERING
934 C> - SUMMARY PRINT FILE
935 C> - (USED ONLY IN PREVENTS MODE)
936 C> @param SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
937 C> - REPORT BEING PROCESSED
938 C> @param NEWTYP - INDICATOR IF THE BUFR MESSAGE TABLE A ENTRY HAS
939 C> - CHANGED FROM THAT OF THE PREVIOUS REPORT (=0 - NO,
940 C> - =1 - YES)
941 C>
942  SUBROUTINE gblevn02(IUNITP,IUNITS,NEWTYP,subset)
943  ! FORMERLY SUBROUTINE FILTAN
944 
945  dimension nflgrt(100:299,12),oemin(2:6)
946  CHARACTER*8 STNID,subset
947  CHARACTER*40 PEVN,QEVN,TEVN,WEVN,PWVN,PW1VN,PW2VN,PW3VN,PW4VN
948  REAL(8) PEV_8(4,255),QEV_8(4,255),TEV_8(4,255),WEV_8(5,255),
949  $ pwv_8(4,255),pw1v_8(4,255),pw2v_8(4,255),
950  $ pw3v_8(4,255),pw4v_8(4,255),obs_8,qms_8,bak_8,sid_8,
951  $ ufc_8,vfc_8
952  LOGICAL FCST,REJP_PS,REJPS,REJT,REJQ,REJW,REJPW,REJPW1,
953  $ rejpw2,rejpw3,rejpw4,satmqc,satemp,soln60,sols60,
954  $ moerr_p,moerr_t,adpupa_virt,doberr,dofcst,some_fcst,
955  $ dovtmp,virt,recalc_q,doprev
956  REAL(8) BMISS
957 
958  COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
959  $ xob,yob,dhr,typ,nlev
960  COMMON /gbevbb/ pvcd,vtcd
961  COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
962  $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
963  COMMON /gbevee/psg01,zsg01,tg01(500),ug01(500),vg01(500),
964  x qg01(500),zint(500),pint(500),pintlog(500),plev(500),
965  x plevlog(500)
966  COMMON /gbevff/ bmiss
967 
968  equivalence(sid_8,stnid)
969 
970  DATA pevn /'POB PQM PPC PRC '/
971  DATA qevn /'QOB QQM QPC QRC '/
972  DATA tevn /'TOB TQM TPC TRC '/
973  DATA wevn /'UOB VOB WQM WPC WRC '/
974  DATA pwvn /'PWO PWQ PWP PWR '/
975  DATA pw1vn /'PW1O PW1Q PW1P PW1R '/
976  DATA pw2vn /'PW2O PW2Q PW2P PW2R '/
977  DATA pw3vn /'PW3O PW3Q PW3P PW3R '/
978  DATA pw4vn /'PW4O PW4Q PW4P PW4R '/
979 
980  DATA nflgrt/2400*0/
981 
982  DATA oemin /0.5,0.1,1.0,0.5,1.0/
983 
984  ni = mod((nint(typ)/10),10)
985 
986  IF(newtyp.EQ.1) nflgrt = 0
987 
988 C LOGICAL SWITCHES FOR OBSERVATION LOCATION FILTERING
989 C ---------------------------------------------------
990 
991  satemp = ((typ.GE.160.AND.typ.LE.179).AND.satmqc)
992  soln60 = ((typ.GE.160.AND.typ.LE.163).AND.yob.GE.-60.AND.satmqc)
993  sols60 = ((typ.EQ.160.OR.typ.EQ.162.OR.typ.EQ.163).AND.yob.LT.-60
994  $ .AND.satmqc)
995 
996 C CLEAR THE EVENT ARRAYS
997 C ----------------------
998 
999  pev_8 = bmiss
1000  qev_8 = bmiss
1001  tev_8 = bmiss
1002  wev_8 = bmiss
1003  pwv_8 = bmiss
1004  pw1v_8 = bmiss
1005  pw2v_8 = bmiss
1006  pw3v_8 = bmiss
1007  pw4v_8 = bmiss
1008 
1009  maxpev = 0
1010  maxqev = 0
1011  maxtev = 0
1012  maxwev = 0
1013  maxpwv = 0
1014  maxpw1v = 0
1015  maxpw2v = 0
1016  maxpw3v = 0
1017  maxpw4v = 0
1018 
1019 C LOOP OVER LEVELS APPLYING UNDERGROUND FILTERING AND SPECIAL RULES
1020 C -----------------------------------------------------------------
1021 
1022  IF(nlev.GT.0) THEN
1023  DO l=1,nlev
1024 
1025  pob = obs_8( 1,l)
1026  qob = obs_8( 2,l)
1027  tob = obs_8( 3,l)
1028  uob = obs_8( 5,l)
1029  vob = obs_8( 6,l)
1030  pwo = obs_8( 7,l)
1031  pw1o = obs_8( 8,l)
1032  pw2o = obs_8( 9,l)
1033  pw3o = obs_8(10,l)
1034  pw4o = obs_8(11,l)
1035  cat = obs_8(12,l)
1036  prss = obs_8(13,l)
1037 
1038  pqm = qms_8( 1,l)
1039  qqm = qms_8( 2,l)
1040  tqm = qms_8( 3,l)
1041  zqm = qms_8( 4,l)
1042  wqm = qms_8( 5,l)
1043  pwq = qms_8( 6,l)
1044  pw1q = qms_8( 7,l)
1045  pw2q = qms_8( 8,l)
1046  pw3q = qms_8( 9,l)
1047  pw4q = qms_8(10,l)
1048 
1049  rejp_ps = .false.
1050  moerr_p = .false.
1051  moerr_t = .false.
1052  rcd = 99999
1053 
1054 C -------------------------------------------------------------------
1055 C RULES FOR PRESSURE (ON ANY LEVEL) -- ALL DATA (MASS AND WIND) ON
1056 C LEVEL REJECTED IF:
1057 C - PRESSURE MORE THAN 100 MB BELOW MODEL (GUESS) SURFACE PRESSURE
1058 C (AND SWITCH FCST=TRUE) -- "PREVENT" PGM REASON CODE 1
1059 C - PRESSURE IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON CODE 2
1060 C - SURFACE (MASS OR WIND) REPORT PRESSURE IS REPORTED ABOVE 450 MB
1061 C OR BELOW 1100 MB -- "PREVENT" PGM REASON CODE 2
1062 C REJECTION MEANS Q.M. SET TO 8
1063 C -------------------------------------------------------------------
1064 
1065  IF(pob.LT.bmiss) THEN
1066  IF(.NOT.fcst) psg01 = pob
1067  IF(pob-psg01.GE.100. .OR. pob.LE.0. .OR.
1068  $ ((pob.LE.450..OR.pob.GE.1100.) .AND. ni.EQ.8)) THEN
1069  IF(pob.LE.0..OR.pob.LE.450..OR.pob.GE.1100.) THEN
1070  IF(ni.EQ.8) THEN
1071  WRITE(iunits,302) stnid,nint(typ),yob,xob,pob
1072  302 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1073  $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, FAILS SANITY ',
1074  $ 'CHECK')
1075  ELSE
1076  WRITE(iunits,101) stnid,nint(typ),yob,xob,pob
1077  101 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1078  $'E, REJECT ALL DATA ON LVL - POB=',f6.1,' MB, FAILS SANITY CHECK')
1079  ENDIF
1080  rcd = 2
1081  ELSE
1082  IF(ni.EQ.8) THEN
1083  WRITE(iunits,303) stnid,nint(typ),yob,xob,pob,psg01
1084  303 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1085  $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, > 100 MB ',
1086  $ 'BELOW GES PSFC(=',f6.1,'MB)')
1087  ELSE
1088  WRITE(iunits,102) stnid,nint(typ),yob,xob,pob,psg01
1089  102 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1090  $ 'E, REJECT ALL DATA ON LVL - POB=',f6.1,' MB, > 100 MB BELOW ',
1091  $ 'GES PSFC(=',f6.1,' MB)')
1092  ENDIF
1093  rcd = 1
1094  ENDIF
1095  rej = 8
1096  rejp_ps = .true.
1097  pev_8(1,l) = pob
1098  pev_8(2,l) = rej
1099  pev_8(3,l) = pvcd
1100  pev_8(4,l) = rcd
1101  maxpev = l
1102  ENDIF
1103  ENDIF
1104 
1105 C -------------------------------------------------------------------
1106 C RULES FOR SURFACE PRESSURE -- ALL MASS DATA ON SURFACE LEVEL
1107 C REJECTED IF:
1108 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1109 C "PREVENT" PGM REASON CODE 3
1110 C - PRESSURE IS MORE THAN 100 MB ABOVE OR BELOW MODEL (GUESS)
1111 C SURFACE PRESSURE (AND SWITCH FCST=TRUE) --
1112 C "PREVENT" PGM REASON CODE 4
1113 C - PRESSURE IS REPORTED ABOVE 450 MB OR BELOW 1100 MB -- "PREVENT"
1114 C PGM REASON CODE 2 (NOTE: DOES NOT APPLY TO SURFACE REPORTS,
1115 C THESE WERE TESTED FOR THIS CRITERION IN ABOVE PRESSURE TEST)
1116 C - PRESSURE VIOLATES RULES FOR PRESSURE ON ANY LEVEL (SEE ABOVE)
1117 C REJECTION FOR FIRST RULE MEANS Q.M. SET TO 9 UNLESS:
1118 C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
1119 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1120 C AS IS
1121 C -------------------------------------------------------------------
1122 
1123  IF(pob.LT.bmiss .AND. cat.EQ.0) THEN
1124  IF(.NOT.fcst) psg01 = pob
1125  rejps = oefg01(pob,typ,5,oemin(5)).GE.bmiss .OR.
1126  $ abs(pob-psg01).GE.100. .OR.
1127  $ pob.LE.450. .OR.
1128  $ pob.GE.1100.
1129  IF(rejps.OR.rejp_ps) THEN
1130  rej = 8
1131  IF(.NOT.rejp_ps) THEN
1132  IF(abs(pob-psg01).GE.100.) THEN
1133  WRITE(iunits,104) stnid,nint(typ),yob,xob,pob,psg01
1134  104 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1135  $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, > 100 MB ',
1136  $ 'ABOVE GES PSFC(=',f6.1,'MB)')
1137  rcd = 4
1138  ELSE IF(pob.LE.450..OR.pob.GE.1100.) THEN
1139  WRITE(iunits,105) stnid,nint(typ),yob,xob,pob
1140  105 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1141  $ 'E, REJECT ALL DATA ON SFC LVL - POB=',f6.1,' MB, FAILS SANITY ',
1142  $ 'CHECK - this should never be printed since test now made in ',
1143  $ 'section above ')
1144  rcd = 2
1145  ELSE
1146  IF(nflgrt(nint(typ),1).EQ.0) THEN
1147  WRITE(iunits,201) nint(typ)
1148  201 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1149  $ 'ALL DATA ON SURFACE LEVEL DUE TO MISSING SFC-P OBS ERROR'/)
1150  nflgrt(nint(typ),1) = 1
1151  ENDIF
1152 CDAK CDAK CDAK CDAK WRITE(IUNITS,103) STNID,NINT(TYP),YOB,XOB,POB
1153 CD103 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1154 CDAK $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,'MB, MISSING OBS.',
1155 CDAK $ ' ERROR')
1156  rcd = 3
1157  rej = 9
1158  ENDIF
1159  ENDIF
1160  rejp_ps = .true.
1161  IF(rcd.EQ.3) moerr_p = .true.
1162  IF(rej.EQ.9.AND.(pqm.GT.3.AND.pqm.LT.15)) THEN
1163 CDAKCDAKCDAKCDAK WRITE(IUNITS,1401) STNID,NINT(TYP),YOB,XOB,PQM
1164  1401 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1165  $ 'E, INPUT PQM =',f4.0,' -- DO NOT APPLY PSFC QM=9 EVENT')
1166  ELSE
1167  pev_8(1,l) = pob
1168  pev_8(2,l) = rej
1169  pev_8(3,l) = pvcd
1170  pev_8(4,l) = rcd
1171  maxpev = l
1172  ENDIF
1173  ENDIF
1174  ENDIF
1175 
1176 C -------------------------------------------------------------------
1177 C RULES FOR TEMPERATURE -- TOB AND QOB ON LEVEL REJECTED IF:
1178 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1179 C "PREVENT" PGM REASON CODE 3
1180 C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
1181 C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
1182 C - REPORT IS TYPE 160-163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE
1183 C SOUNDINGS, ALL PATHS), AND IS AT OR NORTH OF 60 DEGREES SOUTH
1184 C LATITUDE, AND PRESSURE ON LEVEL IS AT OR BELOW 100 MB (AND
1185 C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 6
1186 C - REPORT IS TYPE 160,162,163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE
1187 C SOUNDINGS, ALL PATHS BUT CLEAR), AND IS SOUTH OF 60 DEGREES
1188 C SOUTH LATITUDE, AND PRESSURE ON LEVEL IS BELOW 100 MB (AND
1189 C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 6
1190 C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
1191 C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
1192 C 2 ABOVE) (SEE ABOVE)
1193 C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
1194 C REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS:
1195 C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
1196 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1197 C AS IS
1198 C -------------------------------------------------------------------
1199 
1200  IF(tob.LT.bmiss) THEN
1201  rejt = oefg01(pob,typ,2,oemin(2)).GE.bmiss .OR.
1202  $ (soln60.AND.nint(pob*10.).GE.1000) .OR.
1203  $ (sols60.AND.nint(pob*10.).GT.1000)
1204  IF(rejt.OR.rejp_ps) THEN
1205  rej = 8
1206  IF(.NOT.rejp_ps) THEN
1207  IF(soln60.AND.nint(pob*10.).GE.1000) THEN
1208  IF(nflgrt(nint(typ),6).EQ.0) THEN
1209  WRITE(iunits,7304) nint(typ)
1210  7304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1211  $'TOB/QOB AT AND BELOW 100 MB IF REPORT IS NORTH OF 60S LATITUDE'/)
1212  nflgrt(nint(typ),6) = 1
1213  ENDIF
1214  rcd = 6
1215  ELSE IF(sols60.AND.nint(pob*10.).GT.1000) THEN
1216  IF(nflgrt(nint(typ),7).EQ.0) THEN
1217  WRITE(iunits,7305) nint(typ)
1218  7305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1219  $ 'TOB/QOB BELOW 100 MB IF REPORT IS SOUTH OF 60S LATITUDE'/)
1220  nflgrt(nint(typ),7) = 1
1221  ENDIF
1222  rcd = 6
1223  ELSE
1224  IF(nflgrt(nint(typ),2).EQ.0) THEN
1225  IF(ni.EQ.8) THEN
1226  WRITE(iunits,304) nint(typ)
1227  304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1228  $ 'TOB/QOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
1229  ELSE
1230  WRITE(iunits,202) nint(typ)
1231  202 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1232  $ 'TOB/QOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ',
1233  $ 'MISSING OBS ERROR'/)
1234  ENDIF
1235  nflgrt(nint(typ),2) = 1
1236 cdak cdak cdak cdak cdakWRITE(IUNITS,106) STNID,NINT(TYP),YOB,XOB,TOB
1237 cd106 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1238 cdak $ 'E, REJECT TOB/QOB ON LVL - TOB=',F5.1,'C, MISSING OBS. ERROR')
1239  ENDIF
1240  rcd = 3
1241  rej = 9
1242  ENDIF
1243  ELSE
1244  IF(moerr_p) THEN
1245  rcd = 3
1246  rej = 9
1247  ENDIF
1248  ENDIF
1249  IF(rcd.EQ.3) moerr_t = .true.
1250  IF(rej.EQ.9.AND.(tqm.GT.3.AND.tqm.LT.15)) THEN
1251 CDAKCDAKCDAKCDAK WRITE(IUNITS,1402) STNID,NINT(TYP),YOB,XOB,TQM
1252  1402 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1253  $ 'E, INPUT TQM =',f4.0,' -- DO NOT APPLY TEMP QM=9 EVENT')
1254  ELSE
1255  tev_8(1,l) = tob
1256  tev_8(2,l) = rej
1257  tev_8(3,l) = pvcd
1258  tev_8(4,l) = rcd
1259  maxtev = l
1260  ENDIF
1261  ENDIF
1262  ENDIF
1263 
1264 C -------------------------------------------------------------------
1265 C RULES FOR SPECIFIC HUMIDITY -- QOB ON LEVEL REJECTED IF:
1266 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1267 C "PREVENT" PGM REASON CODE 3
1268 C - PRESSURE ON LEVEL IS ABOVE "QTOP_REJ" MB {WHERE QTOP_REJ IS
1269 C READ IN FROM NAMELIST "PREVDATA" (SEE DOCBLOCK)} -- "PREVENT"
1270 C PGM REASON CODE 5
1271 C - OBSERVATION ERROR FOR TEMPERATURE ON LEVEL IS MISSING (AND
1272 C SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
1273 C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
1274 C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
1275 C - TEMPERATURE ON LEVEL IS MISSING OR IS LESS THAN -150 DEG. C --
1276 C "PREVENT" PGM REASON CODE 2
1277 C - SPECIFIC HUMIDITY IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON
1278 C CODE 2
1279 C - REPORT IS TYPE 160-179 (SATELLITE TEMPERATURE SOUNDINGS, ALL
1280 C TYPES, ALL PATHS, LAND AND SEA), ALL PRESSURE LEVELS (AND
1281 C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 7
1282 C - TEMPERATURE ON LEVEL VIOLATES RULES FOR TEMPERATURE (EXCEPT FOR
1283 C MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE 2 ABOVE)
1284 C (SEE ABOVE)
1285 C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
1286 C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
1287 C 4 ABOVE) (SEE ABOVE)
1288 C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
1289 C REJECTION FOR FIRST FOUR RULES MEANS Q.M. SET TO 9 UNLESS:
1290 C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
1291 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1292 C AS IS
1293 C -------------------------------------------------------------------
1294 
1295  IF(qob.LT.bmiss) THEN
1296  rejq = oefg01(pob,typ,3,oemin(3)).GE.bmiss .OR.
1297  $ tob.GE.bmiss .OR.
1298  $ tob.LE.-150. .OR.
1299  $ nint(pob * 10.).LT.nint(qtop_rej * 10.) .OR.
1300  $ qob.LE.0. .OR.
1301  $ satemp .OR.
1302  $ rejt
1303  IF(rejq.OR.rejp_ps) THEN
1304  rej = 8
1305  IF(.NOT.rejp_ps.AND..NOT.rejt) THEN
1306  IF(satemp) THEN
1307  IF(nflgrt(nint(typ),8).EQ.0) THEN
1308  WRITE(iunits,7306) nint(typ)
1309  7306 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1310  $ 'QOB ON ALL LEVELS'/)
1311  nflgrt(nint(typ),8) = 1
1312  ENDIF
1313  rcd = 7
1314  ELSE IF(qob.LE.0..OR.tob.GE.bmiss.OR.tob.LE.-150.)THEN
1315  WRITE(iunits,111) stnid,nint(typ),yob,xob,
1316  $ qob/1000.,tob
1317  111 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1318  $ 'E, REJECT QOB ON LVL - QOB=',f6.3,' G/KG, FAILS SANITY CHECK ',
1319  $ '(TOB=',f5.1,' C)')
1320  rcd = 2
1321  ELSE IF(oefg01(pob,typ,3,oemin(3)).GE.bmiss) THEN
1322  IF(nflgrt(nint(typ),3).EQ.0) THEN
1323  IF(ni.EQ.8) THEN
1324  WRITE(iunits,305) nint(typ)
1325  305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1326  $ 'QOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
1327  ELSE
1328  WRITE(iunits,203) nint(typ)
1329  203 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1330  $ 'QOB ON AT LEAST ONE LEVEL (IF AVAILABLE ON THAT LEVEL) DUE TO ',
1331  $ 'MISSING OBS ERROR'/)
1332  ENDIF
1333  nflgrt(nint(typ),3) = 1
1334  ENDIF
1335  rcd = 3
1336  rej = 9
1337 cdak cdak cdak cdak WRITE(IUNITS,108) STNID,NINT(TYP),YOB,XOB,QOB/1000.
1338 cd108 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1339 cdak $'E, REJECT QOB ON LVL - QOB=',F6.3,'G/KG, MISSING OBS. ERROR')
1340  ELSE
1341  WRITE(iunits,109) stnid,nint(typ),yob,xob,
1342  $ qob/1000.,qtop_rej,pob
1343  109 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1344  $ 'E, REJECT QOB ON LVL - QOB=',f6.3,' G/KG, ABOVE ',f6.1,
1345  $ 'MB (POB=',f6.1,' MB)')
1346  rcd = 5
1347  rej = 9
1348  ENDIF
1349  ELSE
1350  IF(moerr_p.OR.moerr_t) THEN
1351  rcd = 3
1352  rej = 9
1353  ENDIF
1354  ENDIF
1355  IF(rej.EQ.9.AND.(qqm.GT.3.AND.qqm.LT.15)) THEN
1356 CDAKCDAKCDAKCDAK WRITE(IUNITS,1403) STNID,NINT(TYP),YOB,XOB,QQM
1357  1403 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1358  $ 'E, INPUT QQM =',f4.0,' -- DO NOT APPLY MSTR QM=9 EVENT')
1359  ELSE
1360  qev_8(1,l) = qob
1361  qev_8(2,l) = rej
1362  qev_8(3,l) = pvcd
1363  qev_8(4,l) = rcd
1364  maxqev = l
1365  ENDIF
1366  ENDIF
1367  ENDIF
1368 
1369 C -------------------------------------------------------------------
1370 C RULES FOR WIND -- UOB AND VOB ON LEVEL REJECTED IF:
1371 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1372 C "PREVENT" PGM REASON CODE 3
1373 C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS
1374 C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3
1375 C (Note: This can currently never happen because earlier check
1376 C for missing obs error for sfc pressure is only done if
1377 C "surface" level is category 0 and this is not possible
1378 C for wind reports.)
1379 C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE
1380 C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE
1381 C 2 ABOVE) (SEE ABOVE)
1382 C - This is a SFCSHP report with calm winds and non-missing
1383 C background u- or v-component wind is .GE. 5 m/sec --
1384 C "PREVENT" PGM REASON CODE 8, Q.M. set to 8
1385 C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE)
1386 C REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS:
1387 C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8
1388 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1389 C AS IS
1390 C -------------------------------------------------------------------
1391 
1392  IF(min(uob,vob).LT.bmiss) THEN
1393  rejw = oefg01(pob,typ,4,oemin(4)).GE.bmiss
1394  IF(rejw.OR.rejp_ps) THEN
1395  rej = 8
1396  IF(.NOT.rejp_ps) THEN
1397  IF(nflgrt(nint(typ),4).EQ.0) THEN
1398  IF(ni.EQ.8) THEN
1399  WRITE(iunits,1304) nint(typ)
1400  1304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1401  $ 'UOB/VOB ON SFC LVL DUE TO MISSING OBS ERROR'/)
1402  ELSE
1403  WRITE(iunits,204) nint(typ)
1404  204 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1405  $ 'UOB/VOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ',
1406  $ 'MISSING OBS ERROR'/)
1407  ENDIF
1408  nflgrt(nint(typ),4) = 1
1409  ENDIF
1410 cdak cdak cdak WRITE(IUNITS,112) STNID,NINT(TYP),YOB,XOB
1411 cd112 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1412 cdak $'E, REJECT UOB/VOB ON LVL - MISSING OBS. ERROR')
1413  rcd = 3
1414  rej = 9
1415  ELSE
1416  IF(moerr_p) THEN ! This currently can never be TRUE
1417  ! since CAT is never 0 for "sfc"
1418  ! level of wind reports
1419  rcd = 3
1420  rej = 9
1421  ENDIF
1422  ENDIF
1423  IF(rej.EQ.9.AND.(wqm.GT.3.AND.wqm.LT.15)) THEN
1424 CDAKCDAKCDAKCDAK WRITE(IUNITS,1404) STNID,NINT(TYP),YOB,XOB,WQM
1425  1404 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1426  $ 'E, INPUT WQM =',f4.0,' -- DO NOT APPLY WIND QM=9 EVENT')
1427  ELSE
1428  wev_8(1,l) = uob
1429  wev_8(2,l) = vob
1430  wev_8(3,l) = rej
1431  wev_8(4,l) = pvcd
1432  wev_8(5,l) = rcd
1433  maxwev = l
1434  ENDIF
1435  ENDIF
1436 c Reject calm sfc marine winds if background u- or v- wind .ge. 5 m/s
1437  if(subset.eq."SFCSHP".and.uob.eq.0..and.vob.eq.0.) then
1438  call ufbint(-iunitp,ufc_8,1,1,iret,'UFC')
1439  call ufbint(-iunitp,vfc_8,1,1,iret,'VFC')
1440 cccccccc if(ibfms(ufc_8).eq.0.or.ibfms(vfc_8).eq.0) then
1441  ! DAK: changed to only allow this test if
1442  ! both UFC and VFC are non-missing
1443  if(ibfms(ufc_8).eq.0.and.ibfms(vfc_8).eq.0) then
1444  if(abs(ufc_8).ge.5..or.abs(vfc_8).ge.5.) then
1445  if(wqm.le.3) then
1446  write(iunits,1504) stnid,nint(typ),ufc_8,vfc_8
1447  1504 FORMAT(/' --> ID ',a8,', (RTP ',i3,' UFC=',f5.2,' VFC=',f5.2,') ',
1448  $ 'NEW WIND EVENT, UOB/VOB CALM (0 M/S) WHILE |UFC| AND/OR |VFC| ',
1449  $ '>= 5 M/S, QM SET TO 8'/)
1450  wev_8(1,1) = uob
1451  wev_8(2,1) = vob
1452  wev_8(3,1) = 8
1453  wev_8(4,1) = pvcd
1454  wev_8(5,1) = 8
1455  maxwev = 1
1456  end if
1457  end if
1458  end if
1459  end if
1460  ENDIF
1461 
1462 C -------------------------------------------------------------------
1463 C RULES FOR TOTAL COLUMN PRECIPITABLE WATER -- PWO REJECTED IF:
1464 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1465 C "PREVENT" PGM REASON CODE 3
1466 C REJECTION MEANS Q.M. SET TO 9 UNLESS:
1467 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1468 C AS IS
1469 C -------------------------------------------------------------------
1470 
1471  IF(pwo.LT.bmiss) THEN
1472  rejpw = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
1473  IF(rejpw) THEN
1474  IF(nflgrt(nint(typ),5).EQ.0) THEN
1475  WRITE(iunits,205) nint(typ)
1476  205 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1477  $ 'PWO DUE TO MISSING OBS ERROR'/)
1478  nflgrt(nint(typ),5) = 1
1479  ENDIF
1480 cdakcdakcdak WRITE(IUNITS,113) STNID,NINT(TYP),YOB,XOB,PWO
1481 cd113 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1482 cdak $'E, REJECT PWO ON LVL - PWO=',F5.1,'MM, MISSING OBS. ERROR')
1483  IF(pwq.GT.3.AND.pwq.LT.15) THEN
1484 CDAKCDAKCDAKCDAK WRITE(IUNITS,1405) STNID,NINT(TYP),YOB,XOB,PWQ
1485  1405 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1486  $ 'E, INPUT PWQ =',f4.0,' -- DO NOT APPLY PWtO QM=9 EVENT')
1487  ELSE
1488  pwv_8(1,l) = pwo
1489  pwv_8(2,l) = 9
1490  pwv_8(3,l) = pvcd
1491  pwv_8(4,l) = 3
1492  maxpwv = l
1493  ENDIF
1494  ENDIF
1495  ENDIF
1496 
1497 C -------------------------------------------------------------------
1498 C RULES FOR LAYER 1 PRECIPITABLE WATER -- PW1O REJECTED IF:
1499 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1500 C "PREVENT" PGM REASON CODE 3
1501 C REJECTION MEANS Q.M. SET TO 9 UNLESS:
1502 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1503 C AS IS
1504 C -------------------------------------------------------------------
1505 
1506  IF(pw1o.LT.bmiss) THEN
1507  rejpw1 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
1508  IF(rejpw1) THEN
1509  IF(nflgrt(nint(typ),9).EQ.0) THEN
1510  WRITE(iunits,206) nint(typ)
1511  206 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1512  $ 'PW1O DUE TO MISSING OBS ERROR'/)
1513  nflgrt(nint(typ),9) = 1
1514  ENDIF
1515 cdakcdakcdak WRITE(IUNITS,114) STNID,NINT(TYP),YOB,XOB,PW1O
1516 cd114 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1517 cdak $'E, REJECT PW1O ON LVL - PW1O=',F5.1,'MM, MISSING OBS. ERROR')
1518  IF(pw1q.GT.3.AND.pw1q.LT.15) THEN
1519 CDAKCDAKCDAKCDAK WRITE(IUNITS,1406) STNID,NINT(TYP),YOB,XOB,PW1Q
1520  1406 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1521  $ 'E, INPUT PW1Q =',f4.0,' -- DO NOT APPLY PW1O QM=9 EVENT')
1522  ELSE
1523  pw1v_8(1,l) = pw1o
1524  pw1v_8(2,l) = 9
1525  pw1v_8(3,l) = pvcd
1526  pw1v_8(4,l) = 3
1527  maxpw1v = l
1528  ENDIF
1529  ENDIF
1530  ENDIF
1531 
1532 C -------------------------------------------------------------------
1533 C RULES FOR LAYER 2 PRECIPITABLE WATER -- PW2O REJECTED IF:
1534 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1535 C "PREVENT" PGM REASON CODE 3
1536 C REJECTION MEANS Q.M. SET TO 9 UNLESS:
1537 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1538 C AS IS
1539 C -------------------------------------------------------------------
1540 
1541  IF(pw2o.LT.bmiss) THEN
1542  rejpw2 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
1543  IF(rejpw2) THEN
1544  IF(nflgrt(nint(typ),10).EQ.0) THEN
1545  WRITE(iunits,207) nint(typ)
1546  207 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1547  $ 'PW2O DUE TO MISSING OBS ERROR'/)
1548  nflgrt(nint(typ),10) = 1
1549  ENDIF
1550 cdakcdakcdak WRITE(IUNITS,115) STNID,NINT(TYP),YOB,XOB,PW2O
1551 cd115 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1552 cdak $'E, REJECT PW2O ON LVL - PW2O=',F5.1,'MM, MISSING OBS. ERROR')
1553  IF(pw2q.GT.3.AND.pw2q.LT.15) THEN
1554 CDAKCDAKCDAKCDAK WRITE(IUNITS,1407) STNID,NINT(TYP),YOB,XOB,PW2Q
1555  1407 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1556  $ 'E, INPUT PW2Q =',f4.0,' -- DO NOT APPLY PW2O QM=9 EVENT')
1557  ELSE
1558  pw2v_8(1,l) = pw2o
1559  pw2v_8(2,l) = 9
1560  pw2v_8(3,l) = pvcd
1561  pw2v_8(4,l) = 3
1562  maxpw2v = l
1563  ENDIF
1564  ENDIF
1565  ENDIF
1566 
1567 C -------------------------------------------------------------------
1568 C RULES FOR LAYER 3 PRECIPITABLE WATER -- PW3O REJECTED IF:
1569 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1570 C "PREVENT" PGM REASON CODE 3
1571 C REJECTION MEANS Q.M. SET TO 9 UNLESS:
1572 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1573 C AS IS
1574 C -------------------------------------------------------------------
1575 
1576  IF(pw3o.LT.bmiss) THEN
1577  rejpw3 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
1578  IF(rejpw3) THEN
1579  IF(nflgrt(nint(typ),11).EQ.0) THEN
1580  WRITE(iunits,208) nint(typ)
1581  208 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1582  $ 'PW3O DUE TO MISSING OBS ERROR'/)
1583  nflgrt(nint(typ),11) = 1
1584  ENDIF
1585 cdakcdakcdak WRITE(IUNITS,116) STNID,NINT(TYP),YOB,XOB,PW3O
1586 cd116 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1587 cdak $'E, REJECT PW3O ON LVL - PW3O=',F5.1,'MM, MISSING OBS. ERROR')
1588  IF(pw3q.GT.3.AND.pw3q.LT.15) THEN
1589 CDAKCDAKCDAKCDAK WRITE(IUNITS,1408) STNID,NINT(TYP),YOB,XOB,PW3Q
1590  1408 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1591  $ 'E, INPUT PW3Q =',f4.0,' -- DO NOT APPLY PW3O QM=9 EVENT')
1592  ELSE
1593  pw3v_8(1,l) = pw3o
1594  pw3v_8(2,l) = 9
1595  pw3v_8(3,l) = pvcd
1596  pw3v_8(4,l) = 3
1597  maxpw3v = l
1598  ENDIF
1599  ENDIF
1600  ENDIF
1601 
1602 C -------------------------------------------------------------------
1603 C RULES FOR LAYER 4 PRECIPITABLE WATER -- PW4O REJECTED IF:
1604 C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) --
1605 C "PREVENT" PGM REASON CODE 3
1606 C REJECTION MEANS Q.M. SET TO 9 UNLESS:
1607 C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M.
1608 C AS IS
1609 C -------------------------------------------------------------------
1610 
1611  IF(pw4o.LT.bmiss) THEN
1612  rejpw4 = oefg01(prss*0.01,typ,6,oemin(6)).GE.bmiss
1613  IF(rejpw4) THEN
1614  IF(nflgrt(nint(typ),12).EQ.0) THEN
1615  WRITE(iunits,209) nint(typ)
1616  209 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',i3,', REJECT ',
1617  $ 'PW4O DUE TO MISSING OBS ERROR'/)
1618  nflgrt(nint(typ),12) = 1
1619  ENDIF
1620 cdakcdakcdak WRITE(IUNITS,117) STNID,NINT(TYP),YOB,XOB,PW4O
1621 cd117 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2,
1622 cdak $'E, REJECT PW4O ON LVL - PW4O=',F5.1,'MM, MISSING OBS. ERROR')
1623  IF(pw4q.GT.3.AND.pw4q.LT.15) THEN
1624 CDAKCDAKCDAKCDAK WRITE(IUNITS,1409) STNID,NINT(TYP),YOB,XOB,PW4Q
1625  1409 FORMAT(' ~~> ID ',a8,' (RTP ',i3,'), LAT=',f6.2,'N, LON=',f6.2,
1626  $ 'E, INPUT PW4Q =',f4.0,' -- DO NOT APPLY PW4O QM=9 EVENT')
1627  ELSE
1628  pw4v_8(1,l) = pw4o
1629  pw4v_8(2,l) = 9
1630  pw4v_8(3,l) = pvcd
1631  pw4v_8(4,l) = 3
1632  maxpw4v = l
1633  ENDIF
1634  ENDIF
1635  ENDIF
1636 
1637  ENDDO
1638  ENDIF
1639 
1640 C APPLY THE PROPER EVENTS
1641 C -----------------------
1642 
1643  IF(maxpev .GT.0) CALL ufbint(iunitp,pev_8, 4,maxpev, iret,pevn)
1644  IF(maxqev .GT.0) CALL ufbint(iunitp,qev_8, 4,maxqev, iret,qevn)
1645  IF(maxtev .GT.0) CALL ufbint(iunitp,tev_8, 4,maxtev, iret,tevn)
1646  IF(maxwev .GT.0) CALL ufbint(iunitp,wev_8, 5,maxwev, iret,wevn)
1647  IF(maxpwv .GT.0) CALL ufbint(iunitp,pwv_8, 4,maxpwv, iret,pwvn)
1648  IF(maxpw1v.GT.0) CALL ufbint(iunitp,pw1v_8,4,maxpw1v,iret,pw1vn)
1649  IF(maxpw2v.GT.0) CALL ufbint(iunitp,pw2v_8,4,maxpw2v,iret,pw2vn)
1650  IF(maxpw3v.GT.0) CALL ufbint(iunitp,pw3v_8,4,maxpw3v,iret,pw3vn)
1651  IF(maxpw4v.GT.0) CALL ufbint(iunitp,pw4v_8,4,maxpw4v,iret,pw4vn)
1652 
1653  RETURN
1654  END
1655 C***********************************************************************
1656 C***********************************************************************
1657 C> GBLEVN03 - INTERPOLATE MODEL DATA (FIRST GUESS OR ANALYSIS) TO OB LOCATIONS
1658 C> @param SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR REPORT BEING PROCESSED
1659 C-----------------------------------------------------------------------
1660  SUBROUTINE gblevn03(SUBSET) ! FORMERLY SUBROUTINE GETFC
1661 
1662  USE gblevn_module
1663 
1664  REAL(8) OBS_8,QMS_8,BAK_8,SID_8
1665  REAL(8) BMISS
1666  CHARACTER*8 SUBSET
1667 
1668  COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
1669  $ xob,yob,dhr,typ,nlev
1670  COMMON /gbevee/psg01,zsg01,tg01(500),ug01(500),vg01(500),
1671  x qg01(500),zint(500),pint(500),pintlog(500),plev(500),
1672  x plevlog(500)
1673  COMMON /gbevff/ bmiss
1674 
1675 
1676  DATA tzero / 273.15 /
1677  DATA betap / .0552 /
1678  DATA beta / .00650 /
1679  DATA rog / 29.261 /
1680 
1681 C CLEAR THE BACKGROUND EVENT ARRAY
1682 C --------------------------------
1683 
1684  bak_8 = bmiss
1685 
1686 C GET GUESS PROFILE AT OB LOCATION
1687 C --------------------------------
1688  CALL gblevn06(xob,yob)
1689 
1690 
1691 C INTERPOLATE GUESS PROFILES TO OB PRESSURES
1692 C ------------------------------------------
1693 
1694  IF(nlev.GT.0) THEN
1695  DO 10 l=1,nlev
1696 
1697  pob = obs_8( 1,l)
1698  qob = obs_8( 2,l)
1699  tob = obs_8( 3,l)
1700  zob = obs_8( 4,l)
1701  uob = obs_8( 5,l)
1702  vob = obs_8( 6,l)
1703  pwo = obs_8( 7,l)
1704  pw1o = obs_8( 8,l)
1705  pw2o = obs_8( 9,l)
1706  pw3o = obs_8(10,l)
1707  pw4o = obs_8(11,l)
1708  cat = obs_8(12,l)
1709 
1710  IF(pob.LE.0. .OR. pob.GE.bmiss) GOTO 10
1711 
1712  poblog = log(pob)
1713 
1714  la = -999
1715  lb = -999
1716  do k=1,kmax-1
1717  if (poblog<=plevlog(k) .and. poblog>plevlog(k+1)) then
1718  la = k
1719  lb = k+1
1720  exit
1721  endif
1722  enddo
1723  if (la > 0) then
1724  wt = (poblog-plevlog(lb)) / (plevlog(la)-plevlog(lb))
1725  else
1726  la = 1
1727  lb = la+1
1728  wt = 0.0
1729  endif
1730 
1731  li=0
1732  do k=1,kmax-1
1733  if (poblog<=pintlog(k) .and. poblog>pintlog(k+1)) then
1734  li = k
1735  exit
1736  endif
1737  enddo
1738 
1739 C SURFACE PRESSURE
1740 C ----------------
1741 
1742  IF(cat.EQ.0 .AND. zob.LT.bmiss) THEN
1743  ts = tg01(1) + (psg01-plev(1))*betap
1744  dz = zob-zsg01
1745  tm = ts - dz*beta*.5
1746  pfc = psg01*exp(-dz/(tm*rog))
1747  ELSE
1748  pfc = bmiss
1749  ENDIF
1750 
1751 C SPECIFIC HUMIDITY
1752 C -----------------
1753 
1754  IF(qob.LT.bmiss.OR.tob.LT.bmiss.OR.typ.EQ.111) THEN
1755 
1756 C (QFC NEEDED BY SYNDATA PROGRAM BUT ONLY FOR REPORT TYPE 111)
1757 
1758  qob = qg01(lb) + (qg01(la)-qg01(lb))*wt
1759  ENDIF
1760 
1761 C TEMPERATURE
1762 C -----------
1763 
1764  IF(tob.LT.bmiss.OR.subset.EQ.'VADWND '.OR.typ.EQ.111) THEN
1765 
1766 C (TFC NEEDED BY CQCVAD AND SYNDATA PROGRAMS, LATTER ONLY FOR REPORT
1767 C TYPE 111)
1768 
1769  IF(pob.GT.plev(1)) THEN
1770  tob = tg01(1) + (pob-plev(1))*betap
1771  ELSE
1772  tob = tg01(lb) + (tg01(la)-tg01(lb))*wt
1773  ENDIF
1774  tob = tob - tzero
1775  ENDIF
1776 
1777 C HEIGHT
1778 C ------
1779 
1780  IF(zob.LT.bmiss) THEN
1781  IF(pob.GT.plev(1)) THEN
1782  tm = tg01(1) + (.5*(pint(1)+pob)-plev(1))*betap
1783  zob = zint(1) - rog*tm*log(pob/pint(1))
1784  ELSE
1785  tm = tg01(lb) + (tg01(la)-tg01(lb))*wt
1786  zob = zint(li) - rog*tm*log(pob/pint(li))
1787  ENDIF
1788  ENDIF
1789 
1790 C U AND V COMPONENTS
1791 C ------------------
1792 
1793  IF(uob.LT.bmiss .OR. vob.LT.bmiss) THEN
1794  uob = ug01(lb) + (ug01(la)-ug01(lb))*wt
1795  vob = vg01(lb) + (vg01(la)-vg01(lb))*wt
1796  ENDIF
1797 
1798 
1799 C PRECIPITABLE WATER
1800 C ------------------
1801 
1802  pwo = bmiss
1803  pw1o = bmiss
1804  pw2o = bmiss
1805  pw3o = bmiss
1806  pw4o = bmiss
1807 
1808 C RELATIVE HUMIDITY
1809 C -----------------
1810 
1811  rho = bmiss
1812 
1813 C SCATTER THE PROPER FIRST GUESS/ANALYSIS VALUES
1814 C ----------------------------------------------
1815 
1816  bak_8(1,l) = pfc
1817  bak_8(2,l) = qob
1818  bak_8(3,l) = tob
1819  bak_8(4,l) = zob
1820  bak_8(5,l) = uob
1821  bak_8(6,l) = vob
1822  bak_8(7,l) = pwo
1823  bak_8(8,l) = pw1o
1824  bak_8(9,l) = pw2o
1825  bak_8(10,l) = pw3o
1826  bak_8(11,l) = pw4o
1827  bak_8(12,l) = rho
1828 
1829  10 ENDDO
1830  ENDIF
1831 
1832  RETURN
1833  END
1834 C***********************************************************************
1835 C***********************************************************************
1836 C> Get observation error
1837  SUBROUTINE gblevn04 ! FORMERLY SUBROUTINE GETOE
1838 
1839  dimension oemin(2:6)
1840  REAL(8) OBS_8,QMS_8,BAK_8,SID_8
1841  REAL(8) BMISS
1842 
1843  COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
1844  $ xob,yob,dhr,typ,nlev
1845  COMMON /gbevff/ bmiss
1846 
1847  DATA oemin /0.5,0.1,1.0,0.5,1.0/
1848 
1849 C CLEAR THE EVENT ARRAY
1850 C ---------------------
1851 
1852  bak_8 = bmiss
1853 
1854 C LOOP OVER LEVELS LOOKING UP THE OBSERVATION ERROR
1855 C -------------------------------------------------
1856 
1857  IF(nlev.GT.0) THEN
1858  DO l=1,nlev
1859 
1860  pob = obs_8( 1,l)
1861  qob = obs_8( 2,l)
1862  tob = obs_8( 3,l)
1863  wob = max(obs_8(5,l),obs_8(6,l))
1864  pwo = obs_8( 7,l)
1865  pw1o = obs_8( 8,l)
1866  pw2o = obs_8( 9,l)
1867  pw3o = obs_8(10,l)
1868  pw4o = obs_8(11,l)
1869  cat = obs_8(12,l)
1870 
1871  IF(cat .EQ.0 ) bak_8( 1,l) = oefg01(pob,typ,5,oemin(5))
1872  IF(qob .LT.bmiss) bak_8( 2,l) = oefg01(pob,typ,3,oemin(3))
1873  IF(tob .LT.bmiss) bak_8( 3,l) = oefg01(pob,typ,2,oemin(2))
1874  IF(wob .LT.bmiss) bak_8( 5,l) = oefg01(pob,typ,4,oemin(4))
1875  IF(pwo .LT.bmiss) bak_8( 6,l) = oefg01(pob,typ,6,oemin(6))
1876  IF(pw1o.LT.bmiss) bak_8( 7,l) = oefg01(pob,typ,6,oemin(6))
1877  IF(pw2o.LT.bmiss) bak_8( 8,l) = oefg01(pob,typ,6,oemin(6))
1878  IF(pw3o.LT.bmiss) bak_8( 9,l) = oefg01(pob,typ,6,oemin(6))
1879  IF(pw4o.LT.bmiss) bak_8(10,l) = oefg01(pob,typ,6,oemin(6))
1880 
1881  ENDDO
1882  ENDIF
1883 
1884  RETURN
1885  END
1886 C***********************************************************************
1887 C***********************************************************************
1888 C> SUBROUTINE GBLEVN06 - 2D LINEAR HORIZONTAL INTERPOLATION
1889 C-----------------------------------------------------------------------
1890  SUBROUTINE gblevn06(XOB,YOB) ! FORMERLY SUBROUTINE HTERP
1891 
1892  USE gblevn_module
1893 
1894  COMMON /gbevee/ psi,zsi,ti(500),ui(500),vi(500),qi(500),
1895  x zint(500),pint(500),pintlog(500),plev(500),plevlog(500)
1896 
1897  DATA rog / 29.261 /
1898 
1899 
1900 C CALCULATE HORIZONTAL WEIGHTS AND INTERPOLATE
1901 C --------------------------------------------
1902 
1903  wx = xob/dlon + 1.0
1904  i0 = wx
1905  i1 = mod(i0,imax) + 1
1906  wx = wx-i0
1907 
1908  wy = (yob+90.)/dlat + 1.0
1909  j0 = wy
1910  j1 = min(j0+1,jmax)
1911  wy = wy-j0
1912 
1913 C HTERP FOR SURFACE HEIGHT
1914 C ------------------------
1915 
1916  p1 = iar12z(i0,j0)
1917  p2 = iar12z(i0,j1)
1918  p3 = iar12z(i1,j0)
1919  p4 = iar12z(i1,j1)
1920  p5 = p1+(p2-p1)*wy
1921  p6 = p3+(p4-p3)*wy
1922  zsi = p5+(p6-p5)*wx
1923 
1924 C HTERP FOR SURFACE PRESSURE
1925 C --------------------------
1926 
1927  p1 = iar13p(i0,j0)
1928  p2 = iar13p(i0,j1)
1929  p3 = iar13p(i1,j0)
1930  p4 = iar13p(i1,j1)
1931  p5 = p1+(p2-p1)*wy
1932  p6 = p3+(p4-p3)*wy
1933  psi = p5+(p6-p5)*wx
1934 
1935 C HTERP FOR UPA T,U,V,Q
1936 C ---------------------
1937 
1938  DO k=1,kmax
1939 
1940  p1 = iar14t(i0,j0,k)
1941  p2 = iar14t(i0,j1,k)
1942  p3 = iar14t(i1,j0,k)
1943  p4 = iar14t(i1,j1,k)
1944  p5 = p1+(p2-p1)*wy
1945  p6 = p3+(p4-p3)*wy
1946  ti(k) = p5+(p6-p5)*wx
1947 
1948  p1 = iar15u(i0,j0,k)
1949  p2 = iar15u(i0,j1,k)
1950  p3 = iar15u(i1,j0,k)
1951  p4 = iar15u(i1,j1,k)
1952  p5 = p1+(p2-p1)*wy
1953  p6 = p3+(p4-p3)*wy
1954  ui(k) = p5+(p6-p5)*wx
1955 
1956  p1 = iar16v(i0,j0,k)
1957  p2 = iar16v(i0,j1,k)
1958  p3 = iar16v(i1,j0,k)
1959  p4 = iar16v(i1,j1,k)
1960  p5 = p1+(p2-p1)*wy
1961  p6 = p3+(p4-p3)*wy
1962  vi(k) = p5+(p6-p5)*wx
1963 
1964  p1 = iar17q(i0,j0,k)
1965  p2 = iar17q(i0,j1,k)
1966  p3 = iar17q(i1,j0,k)
1967  p4 = iar17q(i1,j1,k)
1968  p5 = p1+(p2-p1)*wy
1969  p6 = p3+(p4-p3)*wy
1970  qi(k) = p5+(p6-p5)*wx
1971 
1972 C Layer Pressure
1973 C --------------
1974 
1975  p1 = iarpsl(i0,j0,k)
1976  p2 = iarpsl(i0,j1,k)
1977  p3 = iarpsl(i1,j0,k)
1978  p4 = iarpsl(i1,j1,k)
1979  p5 = p1+(p2-p1)*wy
1980  p6 = p3+(p4-p3)*wy
1981  plev(k) = p5+(p6-p5)*wx
1982 
1983 C Interface Pressure
1984 C ------------------
1985 
1986  p1 = iarpsi(i0,j0,k+1)
1987  p2 = iarpsi(i0,j1,k+1)
1988  p3 = iarpsi(i1,j0,k+1)
1989  p4 = iarpsi(i1,j1,k+1)
1990  p5 = p1+(p2-p1)*wy
1991  p6 = p3+(p4-p3)*wy
1992  pint(k+1) = p5+(p6-p5)*wx
1993 
1994  ENDDO
1995 
1996 C Compute interface heights
1997 C -------------------------
1998 
1999  zint(1) = zsi
2000  pint(1) = psi
2001  pintlog(1) = log(pint(1))
2002  do k=2,kmax
2003  k0 = k-1
2004  zint(k) = zint(k0) - rog*ti(k0)*log(pint(k)/pint(k0))
2005  pintlog(k) = log(pint(k))
2006  enddo
2007  pint(kmax+1) = 0.0
2008 
2009 C Compute log(pressure) at layer midpoints
2010 C ----------------------------------------
2011 
2012  do k=1,kmax
2013  plevlog(k) = log(plev(k))
2014  enddo
2015 
2016 ccccc print *,' pint=',pint(1:kmax)
2017 ccccc print *,' zint=',zint(1:kmax)
2018 
2019  RETURN
2020  END
2021 
2022 C$$$ SUBPROGRAM DOCUMENTATION BLOCK
2023 C>
2024 C> SUBPROGRAM: OEFG01
2025 C> PRGMMR: D. A. KEYSER ORG: NP22 DATE: 2007-09-14
2026 C>
2027 C> ABSTRACT: FUNCTION WHICH RETURNS THE OBSERVATION ERROR FOR A
2028 C> REQUESTED VARIABLE INTERPOLATED TO A DEFINED PRESSURE LEVEL FOR A
2029 C> DEFINED REPORT TYPE. IT IS OBTAINED FROM AN INPUT ARRAY CONTAINING
2030 C> OBSERVATION ERRORS ON FIXED PRESSURE LEVELS BY VARIABLE AND REPORT
2031 C> TYPE (READ EARLIER FROM THE EXTERNAL OBSERVATION ERROR TABLE)
2032 C>
2033 C> PROGRAM HISTORY LOG:
2034 C> 1995-05-17 J. WOOLLEN (NP20) - ORIGINAL AUTHOR (FUNCTION OEF)
2035 C> 2007-09-14 D. A. KEYSER -- MODIFIED TO USE EXACT LOGIC AS IN GSI
2036 C> (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE TYPE,
2037 C> LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO MB, A
2038 C> FEW OTHER MINOR CHANGES)
2039 C>
2040 C> USAGE: XX = OEFG01(P,TYP,IE,OEMIN)
2041 C> INPUT ARGUMENT LIST:
2042 C> @param P - REAL PRESSURE LEVEL (MB) TO INTERPOLATE OBS ERROR TO
2043 C> @param TYP - REAL PREPBUFR REPORT TYPE
2044 C> @param IE - VARIABLE TYPE BEING INTERPOLATED (=2 - TEMPERATURE,
2045 C> - =3 - MOISTURE, =4 - WIND, =5 - SURFACE PRESSURE, =6 -
2046 C> - PRECIPITABLE WATER)
2047 C> - (USED ONLY IN PREVENTS MODE)
2048 C> @param OEMIN - REAL MINIMUM VALUE FOR OBS ERROR (FOR VARIABLE BEING
2049 C> - INTERPOLATED)
2050 C>
2051 C> REMARKS: 'OEFG01' RETURNED IS OBSERVATION ERROR FOR VARIABLE "IE" IN
2052 C> REPORT TYPE "TYP", INTERPOLATED TO PRESSURE LEVEL "P".
2053 C>
2054 C>
2055 C> ATTRIBUTES:
2056 C> LANGUAGE: FORTRAN 90
2057 C> MACHINE: NCEP WCOSS
2058 C>
2059  FUNCTION oefg01(P,TYP,IE,OEMIN)
2060 
2061  REAL(8) bmiss
2062 
2063  COMMON /gbevdd/errs(300,33,6)
2064  COMMON /gbevff/ bmiss
2065 
2066  oefg01 = bmiss
2067 
2068 C LOOK UP ERRORS FOR PARTICULAR OB TYPES
2069 C --------------------------------------
2070 
2071  kx = typ
2072 
2073  p = max(0.,min(p,2000.))
2074 
2075  IF(p.GE.errs(kx,1,1)) k1 = 1
2076 
2077  DO kl = 1,32
2078  IF(p.GE.errs(kx,kl+1,1).AND.p.LE.errs(kx,kl,1)) k1 = kl
2079  ENDDO
2080 
2081  IF(p.LE.errs(kx,33,1)) k1 = 5
2082 
2083  k2 = k1 + 1
2084 
2085  ediff = errs(kx,k2,1) - errs(kx,k1,1)
2086  IF(abs(ediff).GT.0.0) THEN
2087  del = (p - errs(kx,k1,1))/ediff
2088  ELSE
2089  del = bmiss
2090  ENDIF
2091 
2092  del = max(0.,min(del,1.0))
2093 
2094  oefg01 = ((1.0 - del) * errs(kx,k1,ie)) + (del * errs(kx,k2,ie))
2095  oefg01 = max(oefg01,oemin)
2096 
2097 C SET MISSING ERROR VALUE TO "BMISS"
2098 C ----------------------------------
2099 
2100  IF(oefg01.GE.5e5) oefg01 = bmiss
2101 
2102  RETURN
2103 
2104  END
2105 
2106 C$$$ SUBPROGRAM DOCUMENTATION BLOCK
2107 C>
2108 C> SUBPROGRAM: GBLEVN08
2109 C> PRGMMR: S. MELCHIOR ORG: NP22 DATE: 2014-03-25
2110 C>
2111 C> ABSTRACT: CREATE VIRTUAL TEMPERATURE EVENTS WITHIN GBLEVENTS
2112 C> SUBROUTINE. FOR ALL TYPES EXCEPT RASS, THIS CONSISTS OF FIRST RE-
2113 C> CALCULATING THE SPECIFIC HUMIDITY FROM THE REPORTED DEWPOINT
2114 C> TEMPERATURE AND PRESSURE, FOLLOWED BY THE CALCULATION OF VIRTUAL
2115 C> TEMPERATURE FROM THE JUST-CALCULATED SPECIFIC HUMIDITY AND THE
2116 C> REPORTED (SENSIBLE) TEMPERATURE. THE RE-CALCULATED SPECIFIC
2117 C> HUMIDITY IS THEN ENCODED AS A STACKED EVENT TO BE LATER WRITTEN
2118 C> INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP", REASON CODE 0).
2119 C> IF THE NAMELIST SWITCH DOVTMP IS TRUE, THEN THE JUST-CALCULATED
2120 C> VIRTUAL TEMPERATURE IS THEN ALSO ENCODED AS A STACKED EVENT TO BE
2121 C> LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP",
2122 C> REASON CODE 0, 2 OR 6). FOR RASS DATA, SPECIFIC HUMIDITY IS
2123 C> MISSING HOWEVER IF THE NAMELIST SWITCH DOVTMP IS TRUE, A SIMPLE
2124 C> COPY OF THE REPORTED (VIRTUAL) TEMPERATURE IS ENCODED AS A STACKED
2125 C> EVENT TO BE LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM
2126 C> "VIRTMP", REASON CODE 3). FOR SURFACE DATA WITH A MISSING PMSL, IF
2127 C> DOVTMP=T AND DOPMSL=T AND A VIRTUAL TEMPERATURE HAS BEEN COMPUTED,
2128 C> CALCULATE AN ESTIMATED PMSL AND ENCODE IT INTO PREPBUFR FILE ALONG
2129 C> WITH AN INDICATOR THAT IS WAS DERIVED HERE. THIS SUBROUTINE IS
2130 C> CURRENTLY ONLY CALLED FOR SURFACE LAND ("ADPSFC"), MARINE
2131 C> ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA") OR SATELLITE
2132 C> TEMPERATURE RETRIEVAL ("SATEMP") DATA TYPES WHEN SWITCH
2133 C> "ADPUPA_VIRT" IS FALSE AND ONLY FOR SURFACE LAND ("ADPSFC"), MARINE
2134 C> ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA"), SATELLITE
2135 C> TEMPERATURE RETRIEVAL ("SATEMP") OR RAOB/DROP/MULTI-LVL RECCO
2136 C> ("ADPUPA") DATA TYPES WHEN SWITCH "ADPUPA_VIRT" IS TRUE. IT IS
2137 C> ALSO ONLY CALLED IN THE PREVENTS MODE. THIS ROUTINE IS CALLED ONCE
2138 C> FOR EACH VALID REPORT IN THE PREPBUFR FILE.
2139 C>
2140 C> PROGRAM HISTORY LOG:
2141 C> 1995-05-17 J. WOOLLEN (NP20) - ORIGINAL AUTHOR
2142 C> 1997-06-01 D.A. KEYSER - STREAMLINED, ADDED SWITCH DOVTMP
2143 C> 1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE
2144 C> NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT
2145 C> TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD
2146 C> QUAL. MARKER (8)
2147 C> 2004-08-30 D. A. KEYSER -- FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY
2148 C> OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF
2149 C> DOVTMP IS TRUE, GETS NEW REASON CODE 3
2150 C> 2006-07-14 D. A. KEYSER -- PROCESSES REPORTS IN MESSAGE TYPE ADPUPA
2151 C> (I.E., RAOBS, DROPS, MULTI-LEVEL RECCOS) WITH SAME LOGIC AS IN
2152 C> SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR WHEN NEW NAMELIST
2153 C> SWITCH "ADPUPA_VIRT" IS TRUE {NORMALLY "ADPUPA_VIRT" IS FALSE
2154 C> (DEFAULT) BECAUSE SUBSEQUENT PROGRAM PREPOBS_CQCBUFR PERFORMS
2155 C> THIS FUNCTION}
2156 C> 2007-09-14 D. A. KEYSER -- FOR NON-"ADPUPA" TYPES, Q.M. 9 IS NOW
2157 C> ASSIGNED TO CALCULATED VIRT. TEMPS IF THE MOISTURE Q.M. IS 9 OR
2158 C> 15 AND ORIG. TEMP NOT "BAD", THESE "VIRTMP" EVENTS RECEIVE NEW
2159 C> REASON CODE 4, HAD RECEIVED Q.M. 8 WITH REASON CODE 2 LIKE VIRT.
2160 C> TEMPS CALCULATED FROM "BAD" MOISTURE - THIS MEANS ONLY TRULY
2161 C> "BAD" VIRT. TEMPS WILL NOW GET Q.M. 8 AND VIRT. TEMPS FLAGGED FOR
2162 C> NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M. 9
2163 C> (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES NOT
2164 C> EVEN CONSIDER OBS WITH Q.M. 8); FOR "ADPUPA" TYPES, Q.M. 3 IS NOW
2165 C> ASSIGNED TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS
2166 C> TRULY BAD (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM
2167 C> IS 1 OR 2 AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN
2168 C> QQM WAS 9 OR 15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES,
2169 C> ENCODES A SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A
2170 C> "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO
2171 C> WHAT IS ALREADY DONE FOR "RASSDA" TYPES)
2172 C> 2013-04-12 D. A. KEYSER -- DON'T ALLOW CALCULATED Q TO BE < 0 WHICH
2173 C> CAN OCCUR ON WCOSS FOR CASES OF HORRIBLY BAD MESONET DATA
2174 C> 2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH,
2175 C> WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION
2176 C> PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST-
2177 C> COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN
2178 C> REPORTED PMSL IS MISSING. DOVTMP MUST BE TRUE AND DOANLS MUST BE
2179 C> FALSE ("PREVENTS" MODE). THE DERIVED PMSL EITHER GETS A QUALITY
2180 C> MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE QUALITY MARK
2181 C> ("PQM"), WHICHEVER IS GREATER. THE VALUE OF THE PMSL INDICATOR
2182 C> (NEW MNEMONIC "PMIN") IS SET TO 1 TO DENOTE PMSL WAS DERIVED
2183 C> RATHER THAN OBSERVED. THE DEFAULT FOR "DOPMSL" IS FALSE
2184 C> (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS). IT IS FORCED TO BE
2185 C> FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE).
2186 C>
2187 C> USAGE: CALL GBLEVN08(IUNITP)
2188 C> INPUT ARGUMENT LIST:
2189 C> @param IUNITP - BUFR OUTPUT FILE UNIT
2190 C> @param SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR
2191 C> - REPORT BEING PROCESSED
2192 C>
2193 C> REMARKS: WILL IMMEDIATELY RETURN TO CALLING PROGRAM IF ANY OF THE
2194 C> FOLLOWING CONDITIONS EXIST: THERE ARE NO LEVELS OF VALID DEWPOINT,
2195 C> OBS, TEMPERATURE Q.M. OR SPEC. HUMIDITY Q.M. IN THE INPUT PREPBUFR
2196 C> FILE FOR THE REPORT. WILL NOT ATTEMPT EITHER SPEC. HUMIDITY NOR
2197 C> VIRT. TEMP CALC. ON A GIVEN LEVEL IF ANY OF THE FOLLOWING
2198 C> CONDITIONS EXIST: REPORTED PRESSURE OBS IS MISSING, REPORTED
2199 C> (SENSIBLE) TEMPERATURE OBS IS MISSING, OR REPORTED DEWPOINT OBS IS
2200 C> MISSING.
2201 C>
2202 C> ATTRIBUTES:
2203 C> LANGUAGE: FORTRAN 90
2204 C> MACHINE: NCEP WCOSS
2205 C>
2206  SUBROUTINE gblevn08(IUNITP,SUBSET) ! FORMERLY SUBROUTINE VTPEVN
2207 
2208  parameter(rd=287.04, g=9.81)
2209  CHARACTER*80 EVNSTQ,EVNSTV,evnstp
2210  CHARACTER*8 SUBSET,stnid
2211  REAL(8) TDP_8(255),TQM_8(255),QQM_8(255),BAKQ_8(4,255),
2212  $ bakv_8(4,255),bakp_8(3),obs_8,qms_8,bak_8,
2213  $ sid_8,pqm_8
2214  real(8) pmo_8,zob_8,pmsl_8
2215  REAL(8) BMISS
2216 
2217  LOGICAL EVNQ,EVNV,DOVTMP,TROP,ADPUPA_VIRT,DOBERR,DOFCST,
2218  $ some_fcst,fcst,virt,satmqc,recalc_q,doprev,
2219  $ evnp,dopmsl,surf
2220 
2221  COMMON /gbevaa/ sid_8,obs_8(13,255),qms_8(12,255),bak_8(12,255),
2222  $ xob,yob,dhr,typ,nlev
2223  COMMON /gbevbb/ pvcd,vtcd
2224  COMMON /gbevcc/ dovtmp,dofcst,some_fcst,doberr,fcst,virt,
2225  $ qtop_rej,satmqc,adpupa_virt,recalc_q,doprev,dopmsl
2226  COMMON /gbevff/ bmiss
2227 
2228  DATA evnstq /'QOB QQM QPC QRC'/
2229  DATA evnstv /'TOB TQM TPC TRC'/
2230  data evnstp /'PMO PMQ PMIN'/
2231 
2232  equivalence(sid_8,stnid)
2233 
2234 C-----------------------------------------------------------------------
2235 C FCNS BELOW CONVERT TEMP/TD (K) & PRESS (MB) INTO SAT./ SPEC. HUM.(G/G)
2236 C-----------------------------------------------------------------------
2237  es(t) = 6.1078*exp((17.269*(t - 273.16))/((t - 273.16)+237.3))
2238  qs(t,p) = (0.622*es(t))/(p-(0.378*es(t)))
2239 C-----------------------------------------------------------------------
2240 C FCN BELOW derives mean sea-level pressure (mb)from P (station
2241 C pressure, mb), Tv (virtual temperature, K), and Z (station height/
2242 C elevation, m) for surface reports
2243 C-----------------------------------------------------------------------
2244  pmsl_fcn(p,tv,z) = p*exp((g*z)/(rd*tv))
2245 C-----------------------------------------------------------------------
2246 
2247 C CLEAR TEMPERATURE AND SPECIFIC HUMIDITY EVENTS
2248 C ----------------------------------------------
2249 
2250  evnq = .false.
2251  evnv = .false.
2252  evnp = .false.
2253  bakq_8 = bmiss
2254  bakv_8 = bmiss
2255  bakp_8 = bmiss
2256  trop = .false.
2257 
2258  surf = (subset.eq.'ADPSFC'.or.subset.eq.'SFCSHP'.or.
2259  $ subset.eq.'MSONET')
2260 
2261 C GET DEWPOINT TEMPERATURE AND CURRENT T,Q QUALITY MARKERS
2262 C Also get mean sea level pressure, station pressure quality mark,
2263 C and station height (elevation) for surface reports if dopmsl=T
2264 C ----------------------------------------------------------------
2265 
2266  CALL ufbint(-iunitp,tdp_8,1,255,nltd,'TDO')
2267  CALL ufbint(-iunitp,qqm_8,1,255,nlqq,'QQM')
2268  CALL ufbint(-iunitp,tqm_8,1,255,nltq,'TQM')
2269  if(surf.and.dopmsl) then
2270  call ufbint(-iunitp,pmo_8,1,1,nlpm,'PMO')
2271  call ufbint(-iunitp,zob_8,1,1,nlzo,'ZOB')
2272  call ufbint(-iunitp,pqm_8,1,1,nlpq,'PQM')
2273  endif
2274  IF(subset.NE.'RASSDA '.AND.subset.NE.'SATEMP ') THEN
2275  IF(nltd.EQ.0) RETURN
2276  IF(nlqq.EQ.0) RETURN
2277  ENDIF
2278  IF(nltq.EQ.0) RETURN
2279  IF(subset.NE.'RASSDA '.AND.subset.NE.'SATEMP ') THEN
2280  IF(nltd.NE.nlev) THEN
2281  print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLTD NLEV - STOP 61")'
2282  CALL errexit(61)
2283  ENDIF
2284  IF(nlqq.NE.nlev) THEN
2285  print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLQQ NLEV - STOP 63")'
2286  CALL errexit(63)
2287  ENDIF
2288  ENDIF
2289  IF(nltq.NE.nlev) THEN
2290  print.NE.'(" ##GBLEVENTS/GBLEVN08 - NLTQ NLEV - STOP 62")'
2291  CALL errexit(62)
2292  ENDIF
2293 
2294 C COMPUTE SPECIFIC HUMIDITY AND VIRTUAL TEMPERATURE USING REPORTED DEWP
2295 C For surface reports, also calculate PMSL if it is missing, dopmsl=T,
2296 C and a virtual temperature has been computed
2297 C ---------------------------------------------------------------------
2298 
2299  IF(nlev.GT.0) THEN
2300  DO l=1,nlev
2301  pob = obs_8(1,l)
2302  tdo = tdp_8(l)
2303  tob = obs_8(3,l)
2304  cat = obs_8(12,l)
2305  IF(dovtmp) THEN
2306  IF(subset.EQ.'RASSDA '.OR.subset.EQ.'SATEMP ') THEN
2307  IF(tob.LT.bmiss) THEN
2308  bakv_8(1,l) = tob
2309  bakv_8(2,l) = tqm_8(l)
2310  bakv_8(3,l) = vtcd
2311  bakv_8(4,l) = 3
2312  evnv = .true.
2313  cycle
2314  ENDIF
2315  ENDIF
2316  ENDIF
2317  IF(pob.LT.bmiss .AND. tob.LT.bmiss
2318  $ .AND. tdo.LT.bmiss) THEN
2319  IF(qqm_8(l).GT.3) THEN
2320 C Don't update q or calculate Tv if bad moisture obs fails sanity check
2321 C (also don't calc. PMSL if it is missing and dopmsl=T for sfc rpts)
2322 cdak IF(TDO.LT.-103.15 .OR. TDO.GT.46.83 .OR. POB.LT.0.1 .OR.
2323 cdak $ POB.GT.1100.)
2324 cdak $ print *, '&&& bad QM fails sanity check'
2325  IF(tdo.LT.-103.15 .OR. tdo.GT.46.83 .OR. pob.LT.0.1 .OR.
2326  $ pob.GT.1100.) cycle
2327  ENDIF
2328  qob = qs(tdo+273.16,pob)
2329 ccccc BAKQ_8(1,L) = QOB*1E6 ! dak fix 2/27/13: can't be > bmiss
2330  ! else flting pt overflow in BUFRLIB
2331 ccc IF(QOB*1E6.LT.BMISS) BAKQ_8(1,L) = QOB*1E6
2332  ! dak add'l fix 4/12/13: don't allow
2333  ! calc. q to be < 0 which can occur
2334  ! on WCOSS for cases of horribly bad
2335  ! mesonet data
2336  IF(qob*1e6.LT.bmiss .AND. qob.GT.0.) bakq_8(1,l) = qob*1e6
2337  bakq_8(2,l) = qqm_8(l) ! Moist qm same as before for
2338  ! re-calc. q
2339  bakq_8(3,l) = vtcd
2340  bakq_8(4,l) = 0 ! Re-calc. q gets unique reason code 0
2341  evnq = .true.
2342 C If message type ADPUPA, test this level to see if at or above trop
2343 C (trop must be above 500 mb to pass test; if no trop level found
2344 C assume it's at 80 mb)
2345 C Don't calculate Tv on this level if at or above trop (doesn't affect
2346 C q calculation)
2347  trop = (subset.EQ.'ADPUPA ' .AND.
2348  $ ((cat.EQ.5 .AND. pob.LT.500.) .OR. pob.LT. 80. .OR. trop))
2349  IF(dovtmp .AND. .NOT.trop) THEN
2350  bakv_8(1,l) = (tob+273.16)*(1.+.61*qob)-273.16
2351  bakv_8(3,l) = vtcd
2352  IF(subset.EQ.'ADPUPA ') THEN
2353 C Message type ADPUPA comes here
2354  IF((qqm_8(l).LT.4.OR.qqm_8(l).EQ.9.OR.qqm_8(l).EQ.15)
2355  $ .OR. tqm_8(l).EQ.0 .OR. tqm_8(l).GT.3
2356  $ .OR. pob.LE.700.) THEN
2357  bakv_8(2,l) = tqm_8(l) ! Tv qm same as for T when
2358  ! q ok or q flagged by
2359  ! PREPRO (but not bad)
2360  bakv_8(4,l) = 0 ! Tv gets unique reason code 0
2361  ELSE
2362  bakv_8(2,l) = 3 !Tv qm susp for bad moist below
2363  ! 700 mb
2364  bakv_8(4,l) = 6 !Tv gets unique reason code 6
2365  ENDIF
2366  ELSE
2367 C All other message types come here
2368  IF(qqm_8(l).LT.4) THEN
2369  bakv_8(2,l) = tqm_8(l) ! Tv qm same as for T when
2370  ! q ok
2371  bakv_8(4,l) = 0 ! Tv gets unique reason code 0
2372  ELSE IF((qqm_8(l).EQ.9.OR.qqm_8(l).EQ.15).AND.
2373  $ (tqm_8(l).LE.3.OR.tqm_8(l).GE.15.OR.
2374  $ tqm_8(l).EQ.9)) THEN
2375 cdak print'(" %%% process tvirt on lvl ",I0," for missing moist obs ",
2376 cdak $ "error/high-up moist case when orig. T not ""bad"" (set TQM_8=",
2377 cdak $ "9)")', l
2378  bakv_8(2,l) = 9 ! Tv qm 9 for moist w/ missing obs
2379  ! error or moist flagged by
2380  ! PREPRO (but not bad) and T qm
2381  ! orig not "bad"
2382  bakv_8(4,l) = 4 ! Tv gets unique reason code 4
2383  ELSE
2384 cdak print'(" %%% process tvirt on lvl ",I0," for ""bad"" QQM_8 case ",
2385 cdak $ "or missing moist obs error/high-up moist w/ ""bad"" TQM_8 case",
2386 cdak $ " (set TQM_8=8)")', l
2387  bakv_8(2,l) = 8 ! Tv qm 8 (bad) for "bad" moist or
2388  ! moist w/ missing obs error or
2389  ! moist flagged by PREPRO (but
2390  ! not bad) and T qm orig "bad"
2391  bakv_8(4,l) = 2 ! Tv gets unique reason code 2
2392  ENDIF
2393  if(surf) then
2394 c For sfc rpts & dopmsl=T, derive Pmsl in cases when it is not reported
2395 c DAK: Note right now this can only happen for sfc rpts where a Tv is
2396 c set to be calculated and where it is successfully calculated.
2397 c Eventually this may be able to be relaxed such that PMSL can be
2398 c derived even if, e.g., no moisture were reported (in this case
2399 c Ts would have to be used, still a decent estimate of PMSL could
2400 c liekly be made). Might also be able to be derived if no T rpted.
2401  if(dopmsl) then
2402  if(ibfms(pmo_8).ne.0) then
2403  tv = bakv_8(1,1) + 273.16
2404  zob=zob_8
2405  pmsl_8=pmsl_fcn(pob,tv,zob)
2406  bakp_8(1) = pmsl_8
2407  bakp_8(2) = max(3,int(pqm_8))! qm>=3 b/c derived
2408  bakp_8(3) = 1. ! pmin=1 for derived value
2409 cccccccccc write(*,4000) stnid,bakp_8(3),bakp_8(2)
2410  4000 format('--> ID ',a8,' Pmsl missing - derived from Pstn; ',
2411  $ 'PMIN = ',f4.1,' PMQ = ',f4.1,'')
2412  evnp = .true.
2413 c Diagnostics for Pmsl values that are suspiciously high
2414  if(pmsl_8.gt.1060) then
2415  write(*,4001) stnid,pob,bakp_8(1)
2416  4001 format('--> ID ',a8,' Derived PMSL unrealistic; FLAG; ',
2417  $ 'POB = ',f7.2,' PMO = ',f7.2,'')
2418  end if
2419  end if
2420  end if
2421  end if
2422  ENDIF
2423  evnv = .true.
2424  ENDIF
2425  ENDIF
2426  ENDDO
2427  ENDIF
2428 
2429 C ENCODE EVENTS INTO REPORT
2430 C -------------------------
2431 
2432  IF(nlev.GT.0) THEN
2433  IF(evnq) CALL ufbint(iunitp,bakq_8,4,nlev,iret,evnstq)
2434  IF(evnv) CALL ufbint(iunitp,bakv_8,4,nlev,iret,evnstv)
2435  if(nlev.eq.1.and.evnp)
2436  $ call ufbint(iunitp,bakp_8,3,nlev,iret,evnstp)
2437  ENDIF
2438 
2439  RETURN
2440  END
2441 C#######################
2442 C#######################
2443 C#######################
2444 C#######################
2445 C#######################
2446 C#######################
2447 C***********************************************************************
2448 C***********************************************************************
2449 C> Do something
2450 C> @param IUNITF - 2-WORD ARRAY:
2451 C> For SIGIO input:
2452 C> - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
2453 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
2454 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
2455 C> THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
2456 C> IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
2457 C> AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
2458 C> IS A MULTIPLE OF 3
2459 C> - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
2460 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
2461 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
2462 C> THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
2463 C> OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
2464 C> THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
2465 C> For NEMSIO input:
2466 C> - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
2467 C> (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
2468 C> THE DATE IN IDATEP
2469 C> - WORD 2 - NOT USED, SHOULD BE A NULL FILE
2470 C>
2471 C> @param IDATEP - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
2472 C> @param im
2473 C> @param jm
2474 C> @param IDRT INTEGER GRID IDENTIFIER
2475 C> (IDRT=4 FOR GAUSSIAN GRID,
2476 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
2477 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
2478  SUBROUTINE gblevn10(IUNITF,IDATEP,IM,JM,IDRT) ! FORMERLY
2479  ! SUBROUTINE GESRES
2480  USE gblevn_module
2481  USE sigio_module
2482  USE sigio_r_module
2483 
2484  IMPLICIT NONE
2485  INTEGER IUNITF(2), IDATEP, IM, JM, IDRT
2486  REAL, PARAMETER :: PI180=.0174532
2487  INTEGER*4, PARAMETER :: ONE=1, ten=10
2488 
2489  TYPE(SIGIO_HEAD) :: HEAD(2)
2490  TYPE(SIGIO_DATS) :: DATS
2491  TYPE(SIGIO_DATM) :: DATM
2492 
2493  INTEGER*4 IRET,IRET1,IRETS,IMJM4,KM4,IDVM,NTRAC,IUNIT4(2)
2494 
2495  INTEGER KFILES,IFILE,JFILE,IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,
2496  $ mdima,mdimb,mdimc,iromb,maxwv,idir,ns,i,j,k,l,ii,jj,ib,ie
2497 
2498  INTEGER IDATE(8,2),JDATE(8,2),KDATE(8,2),KINDX(2)
2499 
2500  CHARACTER*6 COORD(3)
2501  CHARACTER*20 CFILE
2502 
2503  REAL FHOUR,RINC(5)
2504 
2505  DATA coord /'SIGMA ','HYBRID','GENHYB'/
2506 
2507  REAL, ALLOCATABLE :: cofs(:,:), cofv(:,:,:)
2508  REAL, ALLOCATABLE :: cofs_f(:,:,:), cofv_f(:,:,:,:)
2509  REAL (KIND=4),ALLOCATABLE :: grds(:,:,:), grdv(:,:,:,:),
2510  $ wrk1(:,:), wrk2(:,:)
2511 
2512  imax = im
2513  jmax = jm
2514  imjm4 = im*jm
2515  iunit4(:) = iunitf(:)
2516 
2517  IF(mod(mod(idatep,100),3).EQ.0) THEN
2518  kfiles = 1
2519  kindx = 0
2520  print 331, mod(idatep,100)
2521  331 FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',i2.2,
2522  $ ') IS A MULTIPLE OF 3 - ONLY ONE SIGIO (SIGMA OR HYBRID) ',
2523  $ 'INPUT GLOBAL FILE IS READ,'/16x,'NO TIME INTERPOLATION OF ',
2524  $ 'SPECTRAL COEFFICIENTS IS PERFORMED'/)
2525  ELSE
2526  kfiles = 2
2527  kindx(1) = mod(mod(idatep,100),3)
2528  kindx(2) = kindx(1) - 3
2529  print 332, mod(idatep,100)
2530  332 FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',i2.2,
2531  $ ') IS NOT A MULTIPLE OF 3 - TWO SPANNING SIGIO (SIGMA OR ',
2532  $ 'HYBRID) INPUT GLOBAL FILES'/16x,'ARE READ AND THE SPECTRAL ',
2533  $ 'COEFFICIENTS ARE INTERPOLATED TO THE PREPBUFR CENTER TIME'/)
2534  ENDIF
2535 
2536 C GET VALID-TIME DATE OF SIGMA OR HYBRID FILE(S), ALSO READ HEADERS
2537 C -----------------------------------------------------------------
2538 
2539  jfile = 0
2540  rinc = 0
2541  DO ifile=1,kfiles
2542  jfile = ifile
2543 
2544  WRITE(cfile,'("fort.",I2.2)') iunitf(ifile)
2545 
2546  CALL sigio_rropen(iunit4(ifile),cfile,iret)
2547  CALL sigio_rrhead(iunit4(ifile),head(ifile),iret1)
2548  print *,' cfile_sig=',cfile,'sigio_rropen: iret=',iret,
2549  $ 'sigio_rrhead: iret1=',iret1
2550 
2551  IF(iret.NE.0) GO TO 903
2552  IF(iret1.NE.0) GO TO 904
2553 
2554  idate(:,ifile) = 0
2555 
2556  idate(1,ifile) = head(ifile)%IDATE(4)
2557  idate(2:3,ifile) = head(ifile)%IDATE(2:3)
2558  idate(5,ifile) = head(ifile)%IDATE(1)
2559 
2560  fhour = head(ifile)%FHOUR
2561  print'(" idate=",I5,7I3.2," fhour=",F7.3)', idate(:,ifile),
2562  $ head(ifile)%fhour
2563 
2564  IF(idate(1,ifile).LT.100) THEN
2565 
2566 C IF 2-DIGIT YEAR FOUND IN GLOBAL SIMGA FILE INITIAL DATE
2567 C (IDATE(1,IFILE)), MUST USE "WINDOWING" TECHNIQUE TO CREATE A 4-DIGIT
2568 C YEAR (NOTE: THE T170 IMPLEMENTATION IN JUNE 1998 WAS TO INCLUDE THE
2569 C WRITING OF A 4-DIGIT YEAR HERE. PRIOR TO THIS, THE YEAR HERE WAS
2570 C 2-DIGIT.)
2571 
2572  print'(" ##GBLEVENTS/GBLEVN10 - 2-DIGIT YEAR FOUND IN ",
2573  $ "SIGIO (SIGMA OR HYBRID) INPUT GLOBAL FILE ",I0,
2574  $ "; INITIAL DATE (YEAR IS: ",I0,")")', ifile,idate(1,ifile)
2575  print'(" - USE WINDOWING TECHNIQUE TO OBTAIN 4-DIGIT",
2576  $ " YEAR")'
2577  IF(idate(1,ifile).GT.20) THEN
2578  idate(1,ifile) = 1900 + idate(1,ifile)
2579  ELSE
2580  idate(1,ifile) = 2000 + idate(1,ifile)
2581  ENDIF
2582  print'(" ##GBLEVENTS/GBLEVN10 - CORRECTED 4-DIGIT YEAR IS",
2583  $ " NOW: ",I0)', idate(1,ifile)
2584  ENDIF
2585 
2586  rinc(2) = fhour
2587  CALL w3movdat(rinc,idate(:,ifile),jdate(:,ifile))
2588 
2589  print 1, ifile,head(ifile)%FHOUR,
2590  $ (idate(ii,ifile),ii=1,3),idate(5,ifile),(jdate(ii,ifile),
2591  $ ii=1,3),jdate(5,ifile)
2592  1 FORMAT(' --> GBLEVENTS: SIGIO (SIGMA OR HYBRID) INPUT GLOBAL ',
2593  $ 'FILE',i2,' HERE IS A ',f5.1,' HOUR FORECAST FROM ',i5.4,
2594  $ 3i3.2,' VALID AT ',i5.4,3i3.2)
2595 
2596  kdate(:,ifile) = jdate(:,ifile)
2597 
2598  IF(kfiles.EQ.2) THEN
2599  rinc(2) = REAL(KINDX(IFILE))
2600  CALL w3movdat(rinc,jdate(:,ifile),kdate(:,ifile))
2601  ENDIF
2602 
2603  idatgs_cor = (kdate(1,ifile) * 1000000) + (kdate(2,ifile) *
2604  $ 10000) + (kdate(3,ifile) * 100) + kdate(5,ifile)
2605 
2606 C VALID DATES MUST MATCH
2607 C ----------------------
2608 
2609  IF(idatep.NE.idatgs_cor) GO TO 901
2610 
2611  ENDDO
2612 
2613 
2614 C EXTRACT HEADER INFO
2615 C -------------------
2616 
2617  jcap = head(1)%JCAP
2618  kmax = head(1)%LEVS
2619  km4 = kmax
2620  idvc = head(1)%IDVC
2621  idsl = head(1)%IDSL
2622  idvm = head(1)%IDVM
2623  ntrac = head(1)%NTRAC
2624  nvcoord = head(1)%NVCOORD
2625  ALLOCATE (vcoord(kmax+1,nvcoord))
2626  vcoord = 0.0
2627  vcoord = head(1)%VCOORD
2628 
2629  sfcpress_id = mod(head(1)%IDVM,ten)
2630  thermodyn_id = mod(head(1)%IDVM/ten,ten)
2631  IF(idvc == 3 .AND. thermodyn_id == 3) THEN
2632  kmaxs = (ntrac+1)*kmax + 2
2633  ELSE
2634  kmaxs = 2*kmax + 2
2635  ntrac = 1
2636  ENDIF
2637 
2638  ALLOCATE (iar12z(im,jm), iar13p(im,jm))
2639  ALLOCATE (iar14t(im,jm,kmax), iar15u(im,jm,kmax),
2640  $ iar16v(im,jm,kmax), iar17q(im,jm,kmax),
2641  $ iarpsl(im,jm,kmax), iarpsi(im,jm,kmax+1))
2642 
2643 
2644  if(idvc.eq.0) idvc = 1 ! Reset IDVC=0 to 1 (sigma coord.)
2645  IF(idvc < 0 .or. idvc > 3) THEN
2646  print *, '##GBLEVENTS/GBLEVN10: INVALID VERT COORD ID (=',idvc
2647  ENDIF
2648 
2649 
2650 C DEFINE THE OTHER RESOLUTION PARAMETERS
2651 C --------------------------------------
2652 
2653  jcap1 = jcap+1
2654  jcap2 = jcap+2
2655  jcap1x2 = jcap1*2
2656  mdima = jcap1*jcap2
2657  mdimb = mdima/2+jcap1
2658  mdimc = mdimb*2
2659  imax = 384
2660  jmax = imax/2+1
2661 
2662  dlat = 180./(jmax-1)
2663  dlon = 360./imax
2664 
2665  print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,coord(idvc)
2666  2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
2667  $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
2668  $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
2669  $ 'COORD: ',a)
2670 
2671  GO TO 902
2672 
2673  901 CONTINUE
2674  print 9901, jfile,(jdate(ii,jfile),ii=1,3),jdate(5,jfile),idatep
2675  9901 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,' DATE',
2676  $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH -OR SPAN- PREPBUFR FILE ',
2677  $ 'CENTER DATE (',i10,') -STOP 68'/)
2678  CALL errexit(68)
2679  903 CONTINUE
2680  print 9903, jfile,iret
2681  9903 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,
2682  $ ' RETURNED FROM SIGIO_RROPEN WITH R.C.',i3,' -STOP 70'/)
2683  CALL errexit(70)
2684  904 CONTINUE
2685  print 9904, jfile,iret1
2686  9904 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',i2,
2687  $ ' RETURNED FROM SIGIO_RRHEAD WITH R.C.',i3,' -STOP 71'/)
2688  CALL errexit(71)
2689  902 CONTINUE
2690  IF(kmax.GT.500) then
2691  print'(" ##GBLEVENTS/GBLEVN10 - KMAX TOO BIG = ",I0,
2692  $ " - UNABLE TO TRANSFORM GLOBAL SIGMA FILE(S) - STOP 69")',
2693  $ kmax
2694  CALL errexit(69)
2695  ENDIF
2696 
2697 C***********************************************************************
2698 C***********************************************************************
2699 
2700 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2701 
2702 C USAGE: CALL SPTEZM(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,WAVE,GRID,IDIR)
2703 C INPUT ARGUMENTS:
2704 C IROMB - INTEGER SPECTRAL DOMAIN SHAPE
2705 C (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
2706 C MAXWV - INTEGER SPECTRAL TRUNCATION
2707 C IDRT - INTEGER GRID IDENTIFIER
2708 C (IDRT=4 FOR GAUSSIAN GRID,
2709 C IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
2710 C IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
2711 C IMAX - INTEGER EVEN NUMBER OF LONGITUDES
2712 C JMAX - INTEGER NUMBER OF LATITUDES
2713 C KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM
2714 C WAVE - REAL (2*MX,KMAX) WAVE FIELD IF IDIR>0
2715 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2716 C GRID - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR<0
2717 C IDIR - INTEGER TRANSFORM FLAG
2718 C (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
2719 C OUTPUT ARGUMENTS:
2720 C WAVE - REAL (2*MX,KMAX) WAVE FIELD IF IDIR<0
2721 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2722 C GRID - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR>0
2723 
2724 
2725 C USAGE: CALL SPTEZMV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
2726 C & WAVED,WAVEZ,GRIDU,GRIDV,IDIR)
2727 C INPUT ARGUMENTS:
2728 C WAVED - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR>0
2729 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2730 C WAVEZ - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0
2731 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2732 C GRIDU - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR<0
2733 C GRIDV - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR<0
2734 C OUTPUT ARGUMENTS:
2735 C WAVED - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR<0
2736 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2737 C WAVEZ - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0
2738 C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
2739 C GRIDU - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR>0
2740 C GRIDV - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR>0
2741 
2742 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2743 
2744  allocate (cofs_f(mdima,kmaxs,2), cofv_f(mdima,kmax,2,2))
2745 
2746  iromb = 0
2747  maxwv = jcap
2748  if (idrt < 0 .or. idrt > 256) idrt = 0
2749  idir = 1
2750 
2751  IF(kindx(1).EQ.0) THEN
2752  kfiles = 1
2753  ELSE
2754  kfiles = 2
2755  ENDIF
2756 
2757 C Allocate for sigio read
2758 C -----------------------
2759 
2760  sfcpress_id = mod(head(1)%IDVM,ten)
2761  thermodyn_id = mod(head(1)%IDVM/ten,ten)
2762 
2763  print *,'in sig sfcpress_id=',sfcpress_id,' thermodyn_id=',
2764  $ thermodyn_id,' ntrac=',ntrac
2765  print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,
2766  $ ' nvcoord=',nvcoord
2767 
2768  DO ifile=1,kfiles
2769  CALL sigio_aldats(head(ifile),dats,irets)
2770  CALL sigio_aldatm(head(ifile),one,km4,datm,irets)
2771  ! Read surface fields
2772  CALL sigio_rrdats(iunit4(ifile),head(ifile),dats,irets)
2773 
2774  IF(irets.NE.0) THEN
2775  print *,' irets from sigio_rrdats = ', irets
2776  RETURN
2777  ENDIF
2778 
2779  DO i=1,mdima
2780  cofs_f(i,1,ifile) = dats%HS(i)
2781  cofs_f(i,2,ifile) = dats%PS(i)
2782  ENDDO
2783 
2784  ! Read fields on levels 1 through kmax
2785  CALL sigio_rrdatm(iunit4(ifile),head(ifile),datm,irets)
2786  IF(irets.NE.0) THEN
2787  print *,' irets from sigio_rrdatm = ', irets
2788  RETURN
2789  ENDIF
2790 ccccc print *,' aft sigio_rrdatm irets=',irets
2791 
2792  ie = kmax + 2
2793  cofs_f(:,3:ie,ifile) = datm%T
2794  DO i=1,ntrac
2795  ib = ie + 1
2796  ie = ib + kmax - 1
2797  cofs_f(:,ib:ie,ifile) = datm%Q(:,1:kmax,i)
2798  ENDDO
2799  cofv_f(:,:,1,ifile) = datm%D
2800  cofv_f(:,:,2,ifile) = datm%Z
2801 
2802  CALL sigio_axdats(dats,irets)
2803  CALL sigio_axdatm(datm,irets)
2804  ENDDO
2805 
2806 ccccc print *,' after sigio_axdatm'
2807 
2808  ALLOCATE (cofs(mdima,kmaxs), cofv(mdima,kmax,2))
2809  ALLOCATE (grds(imax,jmax,kmaxs), grdv(imax,jmax,kmax,2))
2810  ALLOCATE (wrk1(imax*jmax,kmax), wrk2(imax*jmax,kmax+1))
2811 
2812  IF(kfiles.EQ.1) THEN
2813  DO i = 1,mdima
2814  cofs(i,1:kmaxs) = cofs_f(i,1:kmaxs,1)
2815  cofv(i,1:kmax,1:2) = cofv_f(i,1:kmax,1:2,1)
2816  ENDDO
2817 
2818  ELSE
2819  cofs=
2820  $ ((abs(kindx(2))*cofs_f(:,:,1)) +(kindx(1)*cofs_f(:,:,2)))/3.
2821  cofv=
2822  $ ((abs(kindx(2))*cofv_f(:,:,:,1))+(kindx(1)*cofv_f(:,:,:,2)))/3.
2823  ENDIF
2824 
2825  DEALLOCATE (cofs_f, cofv_f)
2826 
2827  CALL sptezm(iromb,maxwv,idrt,imax,jmax,kmaxs,cofs,grds,idir)
2828  CALL sptezmv(iromb,maxwv,idrt,imax,jmax,kmax,
2829  & cofv(1,1,1),cofv(1,1,2),grdv(1,1,1,1),grdv(1,1,1,2),idir)
2830 
2831  IF( sfcpress_id == 2 ) THEN ! for enthalpy version
2832  grds(:,:,2) = 1000.0*grds(:,:,2) ! Now in Pa
2833  ELSE
2834  grds(:,:,2) = 1000.0*exp(grds(:,:,2)) ! Now in Pa
2835  ENDIF
2836 
2837  DO ns=1, kmaxs
2838  CALL gblevn11(imax,jmax,grds(1,1,ns))
2839  ENDDO
2840  DO j=1,jmax
2841  DO i=1,imax
2842  iar12z(i,j) = grds(i,j,1) ! Orography
2843  iar13p(i,j) = grds(i,j,2) * 0.01 ! Surface Pressure in hPa
2844  ENDDO
2845  ENDDO
2846 
2847  IF(thermodyn_id == 3 .AND. idvc == 3) THEN
2848  grds(:,:,3:2+kmax) = grds(:,:,3:2+kmax) / head(1)%CPI(1)
2849  print *,' cpi(0)=',head(1)%cpi(1)
2850  ENDIF
2851  CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
2852  $ grds(1,1,2),grds(1,1,3),pm=wrk1,pd=wrk2(1,2))
2853 
2854  DO j=1,jmax
2855  jj = (j-1)*imax
2856  DO i=1,imax
2857  wrk2(i+jj,1) = grds(i,j,2) ! in Pa
2858  ENDDO
2859  ENDDO
2860  DO l=1,kmax
2861  wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1) ! in Pa
2862  ENDDO
2863 
2864 ccccc print *,' wrk1=',wrk1(1001,:)
2865 ccccc print *,' wrk2=',wrk2(1001,:)
2866 
2867 ccccc CALL GBLEVN11(IMAX,JMAX,WRK2(1,KMAX+1))
2868 ccccc DO L=1,KMAX
2869 ccccc CALL GBLEVN11(IMAX,JMAX,WRK1(1,L))
2870 ccccc CALL GBLEVN11(IMAX,JMAX,WRK2(1,L))
2871 ccccc ENDDO
2872 
2873  IF(thermodyn_id == 3 .AND. idvc == 3) THEN
2874  ! Convert from enthalpy to temperature
2875  grds(:,:,3:2+kmax) = grds(:,:,3:2+kmax) * head(1)%CPI(1)
2876  CALL sigio_cnvtdv(imjm4,imjm4,km4,idvc,idvm,ntrac,iret,
2877  $ grds(1,1,3),grds(1,1,3+kmax),head(1)%CPI,1_4)
2878  ! Convert back to virtual temperature
2879  grds(:,:,3:kmax+2) = grds(:,:,3:kmax+2) *
2880  $ (1.+(461.50/287.05-1)*grds(:,:,3+kmax:2+kmax*2))
2881  ENDIF
2882 
2883  DO l=1,kmax
2884  DO k=1,2
2885  CALL gblevn11(imax,jmax,grdv(1,1,l,k))
2886  ENDDO
2887  DO j=1,jmax
2888  jj = (j-1)*imax
2889  DO i=1,imax
2890  iar14t(i,j,l) = grds(i,j,2+l) ! Temp (virtual)
2891  iar15u(i,j,l) = grdv(i,j,l,1) ! U component
2892  iar16v(i,j,l) = grdv(i,j,l,2) ! V component
2893  ! specific humidity
2894  iar17q(i,j,l) = max(0.0,grds(i,j,2+kmax+l)*1.0e6)
2895  iarpsl(i,j,l) = wrk1(i+jj,l)*0.01 ! 3D layer pres(hPa)
2896  ENDDO
2897  ENDDO
2898  ENDDO
2899  DO l=1,kmax+1
2900  DO j=1,jmax
2901  jj = (j-1)*imax
2902  DO i=1,imax
2903  iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
2904  ! (hPa)
2905  ENDDO
2906  ENDDO
2907  ENDDO
2908 
2909 
2910 ccccc print *,' iar14t=',iar14t(1,80,:)
2911 ccccc print *,' iar15u=',iar15u(1,80,:)
2912 ccccc print *,' iar16v=',iar16v(1,80,:)
2913 ccccc print *,' iarpsi=',iarpsi(1,80,:)
2914 ccccc print *,' iarpsl=',iarpsl(1,80,:)
2915 
2916  DEALLOCATE (cofs, cofv)
2917  DEALLOCATE (grds, grdv, wrk1, wrk2)
2918 
2919  print *,' RETURNING from GBLENV10'
2920 
2921  RETURN
2922 
2923  END
2924 C***********************************************************************
2925 C***********************************************************************
2926 C> North-south swap
2927  subroutine gblevn11(imax,jmax,grid) ! formerly subroutine n_s_swap
2928  implicit none
2929  integer imax, jmax
2930  real grid(imax,jmax)
2931  real temp (imax)
2932  integer i, j, jj
2933 
2934  do j=1,jmax/2
2935  jj = jmax-j+1
2936  do i=1,imax
2937  temp(i) = grid(i,j)
2938  grid(i,j) = grid(i,jj)
2939  grid(i,jj) = temp(i)
2940  enddo
2941  enddo
2942  return
2943  end
2944 C***********************************************************************
2945 C***********************************************************************
2946 C> Does something
2947  subroutine gblevn11d(imax,jmax,grid)
2948  implicit none
2949  integer imax, jmax
2950  real(kind=8) grid(imax,jmax)
2951  real(kind=8) temp (imax)
2952  integer i, j, jj
2953 
2954  do j=1,jmax/2
2955  jj = jmax-j+1
2956  do i=1,imax
2957  temp(i) = grid(i,j)
2958  grid(i,j) = grid(i,jj)
2959  grid(i,jj) = temp(i)
2960  enddo
2961  enddo
2962  return
2963  end
2964 !---------------------------------------------------------
2965 C> Does something
2966 C> @param iunitf
2967 C> @param idatep
2968 C> @param im
2969 C> @param jm
2970 C> @param IUNITF - 2-WORD ARRAY:
2971 C> For SIGIO input:
2972 C> - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
2973 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
2974 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
2975 C> THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
2976 C> IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
2977 C> AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
2978 C> IS A MULTIPLE OF 3
2979 C> - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
2980 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
2981 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
2982 C> THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
2983 C> OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
2984 C> THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
2985 C> For NEMSIO input:
2986 C> - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
2987 C> (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
2988 C> THE DATE IN IDATEP
2989 C> - WORD 2 - NOT USED, SHOULD BE A NULL FILE
2990 C> @param IDATEP - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
2991 C> @param im
2992 C> @param jm
2993 C> @param IDRT INTEGER GRID IDENTIFIER
2994 C> (IDRT=4 FOR GAUSSIAN GRID,
2995 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
2996 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
2997  SUBROUTINE gblevn12(IUNITF,IDATEP,IM,JM,IDRT)
2998 !---------------------------------------------------------
2999 !--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NEMSIO
3000 !--AND SUBSEQUENTLY MODIFIED BY FANGLIN YANG AND RUSS TREADON
3001 
3002  USE gblevn_module
3003  USE nemsio_module
3004  USE nemsio_openclose
3005  USE nemsio_read
3006  USE nemsio_write
3007  use sigio_module
3008 
3009  IMPLICIT NONE
3010  INTEGER IUNITF(2), IDATEP, IM, JM, IDRT
3011  REAL(KIND=8),PARAMETER:: CON_RV=4.6150e+2,con_rd=2.8705e+2,
3012  $ fv=con_rv/con_rd-1.,oner=1.0,qmin=1.0e-10
3013  INTEGER*4, PARAMETER :: TEN=10
3014 
3015  INTEGER*4 IDRTNEMS
3016  INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC
3017 
3018  INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC,
3019  $ iromb,maxwv,idir,i,j,k,l,ii,jj
3020 
3021  INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT
3022  INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4
3023  REAL(8) tfac
3024  REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:)
3025  REAL,ALLOCATABLE :: tmp(:)
3026  TYPE(NEMSIO_GFILE) :: GFILE
3027 
3028  INTEGER IDATE(8),JDATE(8)
3029  CHARACTER*20 CFILE2
3030  REAL FHOUR,RINC(5)
3031  REAL (KIND=4),ALLOCATABLE :: psfc(:,:), tv(:,:,:),
3032  $ wrk1(:,:), wrk2(:,:)
3033 
3034  imax = im
3035  jmax = jm
3036  imjm4 = im*jm
3037 
3038 ! no time interpolation needed for nemsio input since files are
3039 ! available every hour - original logic to perform time interpolation
3040 ! is flawed and has been removed in this revision
3041 ! --------------------------------------------------------------------
3042 
3043  print 331, mod(idatep,100)
3044  331 format(/' --> GBLEVENTS: ONLY ONE NEMSIO INPUT GLOBAL FILE IS ',
3045  $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ',
3046  $ 'INTERPOLATE'/16x,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR',
3047  $ ' (',i2.2,') IS NOT A MULTIPLE OF 3)'/)
3048 
3049 C GET VALID-TIME DATE OF NEMSIO INPUT FILE, ALSO READ HEADERS
3050 C -----------------------------------------------------------
3051  CALL nemsio_init()
3052 
3053  rinc = 0
3054  idate = 0
3055  WRITE(cfile2,'("fort.",I2.2)') iunitf(1)
3056 
3057  CALL nemsio_open(gfile,trim(cfile2),'read',iret=iret)
3058  print *,' cfile_nems2=',cfile2,'nemsio open,iret=',iret
3059 
3060  idrtnems=idrt
3061  CALL nemsio_getfilehead(gfile,idate=idate7,
3062  & jcap=jcap4,dimx=im4,dimy=jm4,
3063  & dimz=dimz4,idvc=idvc4,ntrac=ntrac,idrt=idrtnems,
3064  & nfhour=nfhour,nfminute=nfminute,nfsecondn=nfsecondn,
3065  & nfsecondd=nfsecondd,idsl=idsl4,idvm=idvm4,iret=iret)
3066  jcap=jcap4
3067  idvc=idvc4
3068  idsl=idsl4
3069  idvm=idvm4
3070  imax=im4
3071  jmax=jm4
3072  kmax=dimz4
3073  if(idrt==-9999) idrt=4 ! default for gaussian grid
3074  idate(1:3)=idate7(1:3)
3075  idate(5)=idate7(4)
3076 
3077  ALLOCATE(vcoord(kmax+1,3))
3078  vcoord=0.0
3079 
3080  IF(nfsecondd/=0.AND. nfsecondd/=-9999) THEN
3081  fhour=REAL(NFHOUR,8)+REAL(NFMINUTE/60.,8)+
3082  $ REAL(NFSECONDN*1./(NFSECONDD*360.),8)
3083  ELSE
3084  fhour=REAL(NFHOUR,8)+REAL(NFMINUTE/60.,8)
3085  ENDIF
3086  print'(" idate=",I5,7I3.2," fhour=",F7.3)', idate,fhour
3087 
3088  ALLOCATE(vcoord4(kmax+1,3,2))
3089  CALL nemsio_getfilehead(gfile,vcoord=vcoord4,iret=iret)
3090  nvcoord=3
3091  IF(maxval(vcoord4(:,3,1))==0..AND.minval(vcoord4(:,3,1))==0.) THEN
3092  nvcoord=2
3093  IF(maxval(vcoord4(:,2,1))==0..AND.minval(vcoord4(:,2,1))==0.)
3094  $ nvcoord=1
3095  ENDIF
3096 
3097  vcoord(1:kmax+1,1:nvcoord)=vcoord4(1:kmax+1,1:nvcoord,1)
3098  DEALLOCATE(vcoord4)
3099 
3100  ALLOCATE(cpi(ntrac+1))
3101  CALL nemsio_getheadvar(gfile,'CPI',cpi,iret=iret)
3102 
3103  rinc(2) = fhour
3104  CALL w3movdat(rinc,idate,jdate)
3105 
3106  print 1, fhour,(idate(ii),ii=1,3),idate(5),(jdate(ii),ii=1,3),
3107  $ jdate(5)
3108  1 FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',f5.1,
3109  $ ' HOUR FORECAST FROM ',i5.4,3i3.2,' VALID AT ',i5.4,3i3.2)
3110 
3111  idatgs_cor = (jdate(1) * 1000000) + (jdate(2) * 10000) +
3112  $ (jdate(3) * 100) + jdate(5)
3113 
3114 C VALID DATES MUST MATCH
3115 C ----------------------
3116 
3117  IF(idatep.NE.idatgs_cor) GO TO 901
3118 
3119 C------------------------------------------
3120 
3121  sfcpress_id = mod(idvm,ten)
3122  thermodyn_id = mod(idvm/ten,ten)
3123 
3124  IF(idvc == 3 .AND. thermodyn_id == 3) THEN
3125  kmaxs = (ntrac+1)*kmax + 2
3126  ELSE
3127  kmaxs = 2*kmax + 2
3128  ntrac = 1
3129  ENDIF
3130 
3131  ALLOCATE (iar12z(imax,jmax), iar13p(imax,jmax))
3132  ALLOCATE (iar14t(imax,jmax,kmax), iar15u(imax,jmax,kmax),
3133  $ iar16v(imax,jmax,kmax), iar17q(imax,jmax,kmax),
3134  $ iarpsl(imax,jmax,kmax), iarpsi(imax,jmax,kmax+1),
3135  $ iarpsd(imax,jmax,kmax))
3136 
3137 !!! DAK: is this right???????
3138  if(idvc.eq.0) idvc = 1 ! Reset IDVC=0 to 1 (sigma coord.)
3139  IF(idvc < 0 .or. idvc > 3) THEN
3140  print *, '##GBLEVENTS/GBLEVN12: INVALID VERT COORD ID (=',idvc
3141  ENDIF
3142 
3143 
3144 C DEFINE THE OTHER RESOLUTION PARAMETERS
3145 C --------------------------------------
3146 
3147  jcap1 = jcap+1
3148  jcap2 = jcap+2
3149  jcap1x2 = jcap1*2
3150  mdima = jcap1*jcap2
3151  mdimb = mdima/2+jcap1
3152  mdimc = mdimb*2
3153 
3154  dlat = 180./(jmax-1)
3155  dlon = 360./imax
3156 
3157  print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,idvc
3158 
3159 ! DAK: is below correct??????????
3160  2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
3161  $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
3162  $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
3163  $ 'COORD ID IS: ',i0,' (not sure what this means!')
3164 
3165 C***********************************************************************
3166 
3167  iromb = 0
3168  maxwv = jcap
3169  if (idrt < 0 .or. idrt > 256) idrt = 0
3170  idir = 1
3171 
3172  print *,'in nem sfcpress_id=',sfcpress_id,' thermodyn_id=',
3173  $ thermodyn_id,' ntrac=',ntrac
3174  print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=',
3175  $ nvcoord
3176 
3177  call w3kind(kindreal,kindint)
3178 
3179 ! print *,'in gblevent, get fields,imax,jmax
3180  allocate(tmp(imax*jmax))
3181 !
3182 ! hgt
3183  k4=1
3184  if(kindreal==4) then
3185  CALL nemsio_readrecvw34(gfile,'hgt','sfc',k4,tmp,iret=iret)
3186  else
3187  CALL nemsio_readrecv(gfile,'hgt','sfc',k4,tmp,iret=iret)
3188  endif
3189 ! print *,'in getnemsio,hgtsfc,',maxval(tmp),minval(tmp),iret
3190  if(iret==0) THEN
3191  iar12z(:,:)=reshape(tmp,(/imax,jmax/))
3192  call gblevn11d(imax,jmax,iar12z)
3193  ENDIF
3194 
3195 !sfc pres, input in Pa, out in hPa
3196  if(kindreal==4) then
3197  CALL nemsio_readrecvw34(gfile,'pres','sfc',k4,tmp,iret=iret)
3198  else
3199  CALL nemsio_readrecv(gfile,'pres','sfc',k4,tmp,iret=iret)
3200  endif
3201 ! print *,'in getnemsio,pressfc,',maxval(tmp),minval(tmp),iret
3202  if(iret==0) THEN
3203  iar13p(:,:)=reshape(tmp*0.01,(/imax,jmax/))
3204  call gblevn11d(imax,jmax,iar13p)
3205  ENDIF
3206 !dpres
3207 ! DO K4=1,KMAX
3208 ! if(kindreal==4) then
3209 ! CALL NEMSIO_READRECVw34(GFILE,'dpres','mid layer',K4,tmp,
3210 ! $ iret=iret)
3211 ! else
3212 ! CALL NEMSIO_READRECV(GFILE,'dpres','mid layer',K4,tmp,
3213 ! $ iret=iret)
3214 ! endif
3215 ! print *,'in getnemsio,dpres,k4=',k4,maxval(tmp),minval(tmp),
3216 ! $ iret
3217 ! if(iret==0) THEN
3218 ! IARPSI(:,:,K4+1)=reshape(tmp*0.01,(/imax,jmax/))
3219 ! ENDIF
3220 ! ENDDO
3221 
3222 !pres
3223 ! DO K4=1,KMAX
3224 ! if(kindreal==4) then
3225 ! CALL NEMSIO_READRECVw34(GFILE,'pres','mid layer',K4,tmp,
3226 ! $ iret=iret)
3227 ! else
3228 ! CALL NEMSIO_READRECV(GFILE,'pres','mid layer',K4,tmp,
3229 ! $ iret=iret)
3230 ! endif
3231 ! print *,'in getnemsio,pres,k4=',k4,maxval(tmp),minval(tmp),iret
3232 ! if(iret==0) THEN
3233 ! IARPSL(:,:,k4)=reshape(tmp*0.01,(/imax,jmax/))
3234 ! ENDIF
3235 ! ENDDO
3236 ! tmp
3237  DO k4=1,kmax
3238  if(kindreal==4) then
3239  CALL nemsio_readrecvw34(gfile,'tmp','mid layer',k4,tmp,
3240  $ iret=iret)
3241  else
3242  CALL nemsio_readrecv(gfile,'tmp','mid layer',k4,tmp,
3243  $ iret=iret)
3244  endif
3245 ! print *,'in getnemsio,tmp,k4=',k4,maxval(tmp),minval(tmp),iret
3246  if(iret==0) THEN
3247  iar14t(:,:,k4)=reshape(tmp,(/imax,jmax/))
3248  call gblevn11d(imax,jmax,iar14t(1,1,k4))
3249  ENDIF
3250  ENDDO
3251 !u
3252  DO k4=1,kmax
3253  if(kindreal==4) then
3254  CALL nemsio_readrecvw34(gfile,'ugrd','mid layer',k4,tmp,
3255  $ iret=iret)
3256  else
3257  CALL nemsio_readrecv(gfile,'ugrd','mid layer',k4,tmp,
3258  $ iret=iret)
3259  endif
3260 ! print *,'in getnemsio,utmp,k4=',k4,maxval(tmp),minval(tmp),iret
3261  if(iret==0) THEN
3262  iar15u(:,:,k4)=reshape(tmp,(/imax,jmax/))
3263  call gblevn11d(imax,jmax,iar15u(1,1,k4))
3264  ENDIF
3265  ENDDO
3266 ! do k=1,kmax
3267 ! print *,'in getnemsio,i15u,k=',k,maxval(iar15u(:,:,k)),
3268 ! $ minval(iar15u(:,:,k)),iret
3269 ! end do
3270 ! print *,'in getnemsio,IAR15U.1,',maxval(IAR15U),minval(IAR15U)
3271 
3272 !v
3273  DO k4=1,kmax
3274  if(kindreal==4) then
3275  CALL nemsio_readrecvw34(gfile,'vgrd','mid layer',k4,tmp,
3276  $ iret=iret)
3277  else
3278  CALL nemsio_readrecv(gfile,'vgrd','mid layer',k4,tmp,
3279  $ iret=iret)
3280  endif
3281  if(iret==0) THEN
3282  iar16v(:,:,k4)=reshape(tmp,(/imax,jmax/))
3283  call gblevn11d(imax,jmax,iar16v(1,1,k4))
3284  ENDIF
3285  ENDDO
3286 ! print *,'in getnemsio,IAR15U.2,',maxval(IAR15U),minval(IAR15U)
3287 !q
3288  DO k4=1,kmax
3289  if(kindreal==4) then
3290  CALL nemsio_readrecvw34(gfile,'spfh','mid layer',k4,tmp,
3291  $ iret=iret)
3292  else
3293  CALL nemsio_readrecv(gfile,'spfh','mid layer',k4,tmp,
3294  $ iret=iret)
3295  endif
3296  if(iret==0) THEN
3297  iar17q(:,:,k4)=max(0.0,reshape(tmp*1.0e6,(/imax,jmax/)) )
3298  call gblevn11d(imax,jmax,iar17q(1,1,k4))
3299  endif
3300  ENDDO
3301 ! print *,'in getnemsio,IAR15U.3,',maxval(IAR15U),minval(IAR15U)
3302 
3303 !
3304  deallocate(tmp)
3305  CALL nemsio_close(gfile,iret=iret)
3306 ! print *,'end of nemsio'
3307 !
3308 
3309 ! print *,'in getnemsio,IAR15U.4,',maxval(IAR15U),minval(IAR15U)
3310 
3311 !--compute Tv from temp
3312  DO k=1,kmax
3313  DO j=1,jmax
3314  DO i=1,imax
3315  tfac=oner+fv*max(iar17q(i,j,k)*1.0e-6,qmin)
3316  iar14t(i,j,k)=iar14t(i,j,k)*tfac
3317  ENDDO
3318  ENDDO
3319  ENDDO
3320 
3321 ! print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U)
3322 
3323 !--COMPUTE PM AND PD
3324  ALLOCATE (psfc(imax,jmax), tv(imax,jmax,kmax))
3325  ALLOCATE (wrk1(imax*jmax,kmax), wrk2(imax*jmax,kmax+1))
3326 
3327  imjm4=imax*jmax
3328  km4=kmax
3329  psfc(:,:) = iar13p(:,:)*100.
3330  tv(:,:,:) = iar14t(:,:,:)
3331 
3332  IF(thermodyn_id == 3 .AND. idvc == 3) THEN
3333  tv(:,:,:) = tv(:,:,:) / cpi(1)
3334  print *,' cpi(1)=',cpi(1)
3335  ENDIF
3336 
3337  CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
3338  $ psfc,tv,pm=wrk1,pd=wrk2(1,2))
3339  DO j=1,jmax
3340  jj = (j-1)*imax
3341  DO i=1,imax
3342  wrk2(i+jj,1) = psfc(i,j) ! in Pa
3343  ENDDO
3344  ENDDO
3345  DO l=1,kmax
3346  wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1) ! in Pa
3347  ENDDO
3348 
3349  DO l=1,kmax
3350  DO j=1,jmax
3351  jj = (j-1)*imax
3352  DO i=1,imax
3353  iarpsl(i,j,l) = wrk1(i+jj,l)*0.01 ! 3D layer pres(hPa)
3354  ENDDO
3355  ENDDO
3356  ENDDO
3357  DO l=1,kmax+1
3358  DO j=1,jmax
3359  jj = (j-1)*imax
3360  DO i=1,imax
3361  iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
3362  ! (hPa)
3363  ENDDO
3364  ENDDO
3365  ENDDO
3366 
3367 ! print*,'GBLEVN12 IARPSI,',maxval(IARPSI),minval(IARPSI)
3368 ! print*,'GBLEVN12 IARPSL,',maxval(IARPSL),minval(IARPSL)
3369 
3370  CALL nemsio_finalize()
3371 !
3372  CALL getlats(idrt)
3373 
3374  RETURN
3375 
3376  901 CONTINUE
3377  print 9901, (jdate(ii),ii=1,3),jdate(5),idatep
3378  9901 FORMAT(/' ##GBLEVENTS/GBLEVN12 - NEMSIO INPUT GLOBAL FILE DATE',
3379  $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ',
3380  $ 'DATE (',i10,') -STOP 68'/)
3381  CALL errexit(68)
3382 
3383  END
3384 ! --------------------------------------------------------------------
3385 C> NetCDF Input
3386 C> @param IUNITF - 2-WORD ARRAY:
3387 C> For SIGIO input:
3388 C> - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL
3389 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
3390 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
3391 C> THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN
3392 C> IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID
3393 C> AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT
3394 C> IS A MULTIPLE OF 3
3395 C> - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL
3396 C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR
3397 C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN
3398 C> THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE
3399 C> OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER
3400 C> THE DATE IN IDATEP THAT IS A MULTIPLE OF 3
3401 C> For NEMSIO input:
3402 C> - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE
3403 C> (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT
3404 C> THE DATE IN IDATEP
3405 C> - WORD 2 - NOT USED, SHOULD BE A NULL FILE
3406 C> @param IDATEP - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH
3407 C> @param IMX
3408 C> @param JMX
3409 C> @param IDRT INTEGER GRID IDENTIFIER
3410 C> (IDRT=4 FOR GAUSSIAN GRID,
3411 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
3412 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
3413  SUBROUTINE gblevn13(IUNITF,IDATEP,IMX,JMX,IDRT)
3414 !---------------------------------------------------------
3415 !--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NETCDF INPUT
3416 
3417  USE gblevn_module
3418  use sigio_module
3419  use netcdf
3420 
3421  IMPLICIT NONE
3422  INTEGER IUNITF(2), IDATEP, IDRT,IMX,JMX
3423  INTEGER yyyy,mm,dd,hh
3424  INTEGER*4 ncid
3425  integer*4 error, id_var, dimid, len
3426  integer*4 im,jm,km, lm,n, k,nargs
3427  REAL(KIND=8),PARAMETER:: CON_RV=4.6150e+2,con_rd=2.8705e+2,
3428  $ fv=con_rv/con_rd-1.,oner=1.0,qmin=1.0e-10
3429  INTEGER*4, PARAMETER :: TEN=10
3430 
3431  INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC
3432 
3433  INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC,
3434  $ iromb,maxwv,idir,i,j,l,ii,jj
3435 
3436  INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT
3437  INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4,kr
3438  REAL(8) tfac, time
3439  REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:),ak(:),bk(:)
3440  REAL,ALLOCATABLE :: temp(:,:), temp3d(:,:,:)
3441 ! CHARACTER,ALLOCATABLE :: attone(:)
3442  character (len = 80) :: attone
3443  INTEGER nt1, nt2
3444  character(len=10) :: dim_nam, grid
3445  CHARACTER*20 gfname
3446 
3447  INTEGER(4) IDATE(8),JDATE(8)
3448  REAL(4) FHOUR,RINC(5)
3449 
3450  REAL (KIND=4),ALLOCATABLE :: psfc(:,:), tv(:,:,:),
3451  $ wrk1(:,:), wrk2(:,:)
3452 
3453  imax = im
3454  jmax = jm
3455 
3456 ! no time interpolation needed for netcdf input since files are
3457 ! available every hour - original logic to perform time interpolation
3458 ! is flawed and has been removed in this revision
3459 ! --------------------------------------------------------------------
3460 
3461  print 331, mod(idatep,100)
3462  331 format(/' --> GBLEVENTS: ONLY ONE NETCDF INPUT GLOBAL FILE IS ',
3463  $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ',
3464  $ 'INTERPOLATE'/16x,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR',
3465  $ ' (',i2.2,') IS NOT A MULTIPLE OF 3)'/)
3466 
3467 C GET VALID-TIME DATE OF NETCDF INPUT FILE, ALSO READ HEADERS
3468 C -----------------------------------------------------------
3469 
3470  WRITE(gfname,'("fort.",I2.2)') iunitf(1)
3471 
3472  error=nf90_open(trim(gfname),nf90_nowrite,ncid)
3473  error=nf90_inq_dimid(ncid,"grid_xt",dimid)
3474  error=nf90_inquire_dimension(ncid,dimid,dim_nam,im)
3475  error=nf90_inq_dimid(ncid,"grid_yt",dimid)
3476  error=nf90_inquire_dimension(ncid,dimid,dim_nam,jm)
3477  error=nf90_inq_dimid(ncid,"pfull",dimid)
3478  error=nf90_inquire_dimension(ncid,dimid,dim_nam,km)
3479  error=nf90_inq_dimid(ncid,"phalf",dimid)
3480  error=nf90_inquire_dimension(ncid,dimid,dim_nam,lm)
3481  print*, "im,jm,kmi,lm:",im,jm,km,lm
3482  print*, "IMX, JMX:", imx, jmx
3483  IF (imx .NE. im .OR. jmx .NE. jm) print*, "Different Resolutions"
3484  imax = im
3485  jmax = jm
3486  kmax = km
3487 
3488 ! error=nf90_inq_varid(ncid, 'pfull', id_var)
3489 ! error=nf90_get_var(ncid, id_var, )
3490 
3491 
3492 ! CALL W3MOVDAT(RINC,IDATE,JDATE)
3493 
3494 
3495 
3496 C VALID DATES MUST MATCH
3497 C ----------------------
3498  error=nf90_inq_varid(ncid, "time", id_var)
3499  error=nf90_inquire_attribute(ncid, id_var, "units", len=len)
3500 ! ALLOCATE (attone(len))
3501  error=nf90_get_att(ncid, id_var, "units", attone)
3502  read(attone(len-18:len-15),'(i4)') yyyy
3503  read(attone(len-13:len-12),'(i2)') mm
3504  read(attone(len-10:len-9),'(i2)') dd
3505  read(attone(len-7:len-6),'(i2)') hh
3506  IF(attone(1:5) .NE. "hours") THEN
3507  print*, "checkout the time unit, not hourly",attone
3508  ELSE
3509  print*, "base time", yyyy,mm,dd,hh
3510  ENDIF
3511 ! DEALLOCATE (attone)
3512  error=nf90_get_var(ncid, id_var, time)
3513  fhour=time
3514 
3515  idate = 0
3516  idate(1)=yyyy
3517  idate(2)=mm
3518  idate(3)=dd
3519  idate(5)=hh
3520 
3521  rinc=0.0
3522  rinc(2)=fhour
3523 
3524  CALL w3movdat(rinc,idate,jdate)
3525 
3526  print 1, fhour,(idate(ii),ii=1,3),idate(5),(jdate(ii),ii=1,3),
3527  $ jdate(5)
3528  1 FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',f5.1,
3529  $ ' HOUR FORECAST FROM ',i5.4,3i3.2,' VALID AT ',i5.4,3i3.2)
3530 
3531  idatgs_cor = (jdate(1) * 1000000) + (jdate(2) * 10000) +
3532  $ (jdate(3) * 100) + jdate(5)
3533 
3534 C VALID DATES MUST MATCH
3535 C ----------------------
3536 
3537  IF(idatep.NE.idatgs_cor) GO TO 901
3538 
3539 C------------------------------------------
3540 
3541  ALLOCATE (iar12z(im,jm), iar13p(im,jm))
3542  ALLOCATE (iar14t(im,jm,km), iar15u(im,jm,km),
3543  $ iar16v(im,jm,km), iar17q(im,jm,km),
3544  $ iarpsl(im,jm,km), iarpsi(im,jm,km+1),
3545  $ iarpsd(im,jm,km))
3546 
3547  error=nf90_get_att(ncid, nf90_global, "grid", grid)
3548  IF (grid == "gaussian")THEN
3549  idrt=4
3550  ENDIF
3551 
3552  nvcoord=2 !for idvc=2, nvcoord=2: hybrid interface a and b
3553  idvc=2 !(1 for sigma and 2 for hybrid)
3554  idsl=1 !(1 for phillips or 2 for mean)
3555  idvm=1
3556  jcap = -9999
3557 
3558  sfcpress_id = mod(idvm,ten)
3559  thermodyn_id = mod(idvm/ten,ten)
3560 
3561  IF(idvc == 3 .AND. thermodyn_id == 3) THEN
3562  kmaxs = (ntrac+1)*kmax + 2
3563  ELSE
3564  kmaxs = 2*kmax + 2
3565  ntrac = 1
3566  ENDIF
3567 
3568 C DEFINE THE OTHER RESOLUTION PARAMETERS
3569 C --------------------------------------
3570 
3571  dlat = 180./(jmax-1)
3572  dlon = 360./imax
3573 
3574  print 2, jcap,imax,jmax,kmax,kmaxs,dlat,dlon,idvc
3575 
3576  2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',i5,' ',
3577  $ i5,' x ',i5,' GRID WITH ',i3,' LEVELS ',i3,
3578  $' SCALARS -------> ',f5.2,' X ',f5.2,' VERT. ',
3579  $ 'COORD ID IS: ',i0)
3580 
3581  print *,'in netcdf sfcpress_id=',sfcpress_id,' thermodyn_id=',
3582  $ thermodyn_id,' ntrac=',ntrac
3583  print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=',
3584  $ nvcoord
3585 
3586  ALLOCATE(vcoord(km+1,3))
3587  vcoord=0.0
3588  allocate(ak(km+1))
3589  allocate(bk(km+1))
3590 
3591  error=nf90_get_att(ncid, nf90_global, "ak", ak)
3592  error=nf90_get_att(ncid, nf90_global, "bk", bk)
3593 
3594  do k=1,km+1
3595  kr=km+2-k
3596  vcoord(k,1) = ak(kr)
3597  vcoord(k,2) = bk(kr)
3598  end do
3599 
3600  deallocate(ak,bk)
3601 
3602  allocate(temp(im,jm))
3603  error=nf90_inq_varid(ncid, 'hgtsfc', id_var)
3604  error=nf90_get_var(ncid, id_var, temp)
3605 
3606  iar12z(:,:)=temp(:,:)
3607  call gblevn11d(imax,jmax,iar12z)
3608 
3609 
3610  error=nf90_inq_varid(ncid, 'pressfc', id_var)
3611  error=nf90_get_var(ncid, id_var, temp)
3612 
3613  iar13p(:,:)=temp(:,:)*0.01
3614  call gblevn11d(imax,jmax,iar13p)
3615 
3616  deallocate(temp)
3617 
3618 !! write(6,111)'hgtsfc', 1,minval(iar12z),maxval(iar12z)
3619 !! write(6,111)'pressfc',1,minval(iar13p),maxval(iar13p)
3620 !! 111 format(a8,1x,i3,1x,2(g13.6,1x))
3621 
3622 
3623  allocate(temp3d(im,jm,km))
3624  error=nf90_inq_varid(ncid, 'tmp', id_var)
3625  error=nf90_get_var(ncid, id_var, temp3d)
3626  DO k4=1,km
3627  kr=km+1-k4
3628  iar14t(:,:,k4)=temp3d(:,:,kr)
3629  call gblevn11d(imax,jmax,iar14t(1,1,k4))
3630 !! write(6,111)'tmp',k4,minval(iar14t(:,:,k4)),
3631 !! & maxval(iar14t(:,:,k4))
3632  ENDDO
3633 
3634  error=nf90_inq_varid(ncid, 'ugrd', id_var)
3635  error=nf90_get_var(ncid, id_var, temp3d)
3636  DO k4=1,km
3637  kr=km+1-k4
3638  iar15u(:,:,k4)=temp3d(:,:,kr)
3639  call gblevn11d(imax,jmax,iar15u(1,1,k4))
3640 !! write(6,111)'ugrd',k4,minval(iar15u(:,:,k4)),
3641 !! & maxval(iar15u(:,:,k4))
3642  ENDDO
3643 
3644  error=nf90_inq_varid(ncid, 'vgrd', id_var)
3645  error=nf90_get_var(ncid, id_var, temp3d)
3646  DO k4=1,km
3647  kr=km+1-k4
3648  iar16v(:,:,k4)=temp3d(:,:,kr)
3649  call gblevn11d(imax,jmax,iar16v(1,1,k4))
3650 !! write(6,111)'vgrd',k4,minval(iar16v(:,:,k4)),
3651 !! & maxval(iar16v(:,:,k4))
3652  ENDDO
3653 
3654  error=nf90_inq_varid(ncid, 'spfh', id_var)
3655  error=nf90_get_var(ncid, id_var, temp3d)
3656  DO k4=1,km
3657  kr=km+1-k4
3658  iar17q(:,:,k4)=max(0.0,temp3d(:,:,kr)*1.e6)
3659  call gblevn11d(imax,jmax,iar17q(1,1,k4))
3660 !! write(6,111)'spfh',k4,minval(iar17q(:,:,k4)),
3661 !! & maxval(iar17q(:,:,k4))
3662  ENDDO
3663 
3664  deallocate(temp3d)
3665 !
3666 
3667 
3668 !--compute Tv from temp
3669  DO k=1,km
3670  DO j=1,jm
3671  DO i=1,im
3672  tfac=oner+fv*max(iar17q(i,j,k)*1.0e-6,qmin)
3673  iar14t(i,j,k)=iar14t(i,j,k)*tfac
3674  ENDDO
3675  ENDDO
3676  ENDDO
3677 
3678 ! print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U)
3679 
3680 
3681 
3682 !--COMPUTE PM AND PD
3683  ALLOCATE (psfc(im,jm), tv(im,jm,km))
3684  ALLOCATE (wrk1(im*jm,km), wrk2(im*jm,km+1))
3685 
3686  imjm4=im*jm
3687  km4=km
3688  psfc(:,:) = iar13p(:,:)*100.
3689  tv(:,:,:) = iar14t(:,:,:)
3690 
3691  IF(thermodyn_id == 3 .AND. idvc == 3) THEN
3692  tv(:,:,:) = tv(:,:,:) / cpi(1)
3693  print *,' cpi(1)=',cpi(1)
3694  ENDIF
3695 
3696  CALL sigio_modpr(imjm4,imjm4,km4,nvcoord,idvc,idsl,vcoord,iret,
3697  $ psfc,tv,pm=wrk1,pd=wrk2(1,2))
3698  DO j=1,jm
3699  jj = (j-1)*im
3700  DO i=1,im
3701  wrk2(i+jj,1) = psfc(i,j) ! in Pa
3702  ENDDO
3703  ENDDO
3704  DO l=1,km
3705  wrk2(:,l+1) = wrk2(:,l) - wrk2(:,l+1) ! in Pa
3706  ENDDO
3707 
3708  DO l=1,km
3709  DO j=1,jm
3710  jj = (j-1)*im
3711  DO i=1,im
3712  iarpsl(i,j,l) = wrk1(i+jj,l)*0.01 ! 3D layer pres(hPa)
3713  ENDDO
3714  ENDDO
3715  ENDDO
3716  DO l=1,km+1
3717  DO j=1,jm
3718  jj = (j-1)*im
3719  DO i=1,im
3720  iarpsi(i,j,l) = wrk2(i+jj,l)*0.01 ! 3D interface pressure
3721  ! (hPa)
3722  ENDDO
3723  ENDDO
3724  ENDDO
3725 
3726  error=nf90_close(ncid)
3727 ! print*,'GBLEVN13 IARPSI,',maxval(IARPSI),minval(IARPSI)
3728 ! print*,'GBLEVN13 IARPSL,',maxval(IARPSL),minval(IARPSL)
3729 
3730 !
3731  CALL getlats(idrt)
3732 
3733  RETURN
3734 
3735  901 CONTINUE
3736  print 9901, idatep,idatgs_cor
3737  9901 FORMAT(/' ##GBLEVENTS/GBLEVN13 - NETCDF INPUT GLOBAL FILE DATE',
3738  $ ' (',i4.4,3(i2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ',
3739  $ 'DATE (',i10,') -STOP 68'/)
3740  CALL errexit(68)
3741 
3742  END
3743 
3744 !-----------------------------------------------------------------
3745 !-----------------------------------------------------------------
3746 C> get latitudes
3747 C> @param IDRT INTEGER GRID IDENTIFIER
3748 C> (IDRT=4 FOR GAUSSIAN GRID,
3749 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
3750 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
3751  subroutine getlats(idrt)
3752 
3753  USE gblevn_module
3754 ! integer :: jmax
3755 ! real,allocatable :: slat(:),wlat(:)
3756 ! real :: rad2deg
3757  real,allocatable :: slats(:)
3758  real(4) slat(jmax),wlat(jmax),rad2deg
3759 !get gaussian or regular latitude array based on idrt
3760 
3761 ! print *,'in getlats,idrt=',idrt
3762  call splat(idrt,jmax,slat,wlat)
3763  rad2deg=180./acos(-1.)
3764  allocate(slats(jmax));
3765  rad2deg=180./acos(-1.)
3766  slats(:)=-asin(slat(:))*rad2deg
3767  dlat=180./float(jmax-1)
3768 
3769  return
3770  end
getlats
subroutine getlats(idrt)
get latitudes
Definition: gblevents.f:3752
gblevn11d
subroutine gblevn11d(imax, jmax, grid)
Does something.
Definition: gblevents.f:2948
gblevn11
subroutine gblevn11(imax, jmax, grid)
North-south swap.
Definition: gblevents.f:2928