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