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