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