Line data Source code
1 : !BOC
2 : #if defined ( MODEL_GCCLASSIC ) || defined( MODEL_WRF ) || defined( MODEL_CESM ) || defined( HEMCO_STANDALONE )
3 : ! The 'standard' HEMCO I/O module is used for:
4 : ! - HEMCO Standalone (HEMCO_STANDALONE)
5 : ! - GEOS-Chem 'Classic' (MODEL_GCCLASSIC)
6 : ! - WRF-GC (MODEL_WRF)
7 : ! - CESM-GC and CAM-Chem / HEMCO-CESM (MODEL_CESM)
8 : !EOC
9 : !------------------------------------------------------------------------------
10 : ! Harmonized Emissions Component (HEMCO) !
11 : !------------------------------------------------------------------------------
12 : !BOP
13 : !
14 : ! !MODULE: hcoio_read_std_mod.F90
15 : !
16 : ! !DESCRIPTION: Module HCOIO\_read\_mod controls data processing
17 : ! (file reading, unit conversion, regridding) for HEMCO in the
18 : ! 'standard' environment (i.e. non-ESMF).
19 : !
20 : ! This module implements the 'standard' environment (i.e. non-ESMF).
21 : !\\
22 : !\\
23 : ! !INTERFACE:
24 : !
25 : MODULE HCOIO_Read_Mod
26 : !
27 : ! !USES:
28 : !
29 : USE HCO_Types_Mod
30 : USE HCO_Error_Mod
31 : USE HCO_CharTools_Mod
32 : USE HCO_State_Mod, ONLY : Hco_State
33 : USE HCOIO_Util_Mod
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 : !
38 : ! !PUBLIC MEMBER FUNCTIONS:
39 : !
40 : PUBLIC :: HCOIO_Read
41 : PUBLIC :: HCOIO_CloseAll
42 : !
43 : ! !REMARKS:
44 : ! Beginning with HEMCO 3.0.0, all I/O modules use the same module names,
45 : ! and their compilation depends on pre-processor flags defined at the top
46 : ! of the file.
47 : !
48 : ! This is to streamline the implementation of one unified Data Input Layer,
49 : ! that can be switched in and out at compile time, and reduce branching of
50 : ! code paths elsewhere.
51 : !
52 : ! !REVISION HISTORY:
53 : ! 22 Aug 2013 - C. Keller - Initial version
54 : ! See https://github.com/geoschem/hemco for complete history
55 : !EOP
56 : !------------------------------------------------------------------------------
57 : !BOC
58 : !
59 : ! !DEFINED PARAMETERS
60 : !
61 : ! Parameter used for difference testing of floating points
62 : REAL(dp), PRIVATE, PARAMETER :: EPSILON = 1.0e-5_dp
63 :
64 : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
65 : REAL(hp), PRIVATE :: GC_72_EDGE_SIGMA(73) = (/ &
66 : 1.000000E+00, 9.849998E-01, 9.699136E-01, 9.548285E-01, 9.397434E-01, 9.246593E-01, &
67 : 9.095741E-01, 8.944900E-01, 8.794069E-01, 8.643237E-01, 8.492406E-01, 8.341584E-01, &
68 : 8.190762E-01, 7.989697E-01, 7.738347E-01, 7.487007E-01, 7.235727E-01, 6.984446E-01, &
69 : 6.733175E-01, 6.356319E-01, 5.979571E-01, 5.602823E-01, 5.226252E-01, 4.849751E-01, &
70 : 4.473417E-01, 4.097261E-01, 3.721392E-01, 3.345719E-01, 2.851488E-01, 2.420390E-01, &
71 : 2.055208E-01, 1.746163E-01, 1.484264E-01, 1.261653E-01, 1.072420E-01, 9.115815E-02, &
72 : 7.748532E-02, 6.573205E-02, 5.565063E-02, 4.702097E-02, 3.964964E-02, 3.336788E-02, &
73 : 2.799704E-02, 2.341969E-02, 1.953319E-02, 1.624180E-02, 1.346459E-02, 1.112953E-02, &
74 : 9.171478E-03, 7.520355E-03, 6.135702E-03, 4.981002E-03, 4.023686E-03, 3.233161E-03, &
75 : 2.585739E-03, 2.057735E-03, 1.629410E-03, 1.283987E-03, 1.005675E-03, 7.846040E-04, &
76 : 6.089317E-04, 4.697755E-04, 3.602270E-04, 2.753516E-04, 2.082408E-04, 1.569208E-04, &
77 : 1.184308E-04, 8.783617E-05, 6.513694E-05, 4.737232E-05, 3.256847E-05, 1.973847E-05, &
78 : 9.869233E-06/)
79 : #endif
80 :
81 : CONTAINS
82 : !EOC
83 : !------------------------------------------------------------------------------
84 : ! Harmonized Emissions Component (HEMCO) !
85 : !------------------------------------------------------------------------------
86 : !BOP
87 : !
88 : ! !IROUTINE: HCOIO_Read
89 : !
90 : ! !DESCRIPTION: Reads a netCDF file and returns the regridded array in proper
91 : ! units. This routine uses the HEMCO generic data reading and regridding
92 : ! routines.
93 : !\\
94 : !\\
95 : ! Two different regridding algorithm are used: NCREGRID for 3D data with
96 : ! vertical regridding, and map\_a2a for all other data. map\_a2a also
97 : ! supports index-based remapping, while this feature is currently not
98 : ! possible in combination with NCREGRID.
99 : !\\
100 : !\\
101 : ! 3D data is vertically regridded onto the simulation grid on the sigma
102 : ! interface levels. In order to calculate these levels correctly, the netCDF
103 : ! vertical coordinate description must adhere to the CF - conventions. See
104 : ! routine NC\_Get\_Sigma\_Levels in Ncdf\_Mod for more details.
105 : !\\
106 : !\\
107 : ! A simpler vertical interpolation scheme is used if (a) the number of
108 : ! vertical levels of the input data corresponds to the number of levels
109 : ! on the simulation grid (direct mapping, no remapping), (b) the vertical
110 : ! level variable name (long\_name) contains the word "GEOS-Chem level". In
111 : ! the latter case, the vertical levels of the input data is interpreted as
112 : ! GEOS vertical levels and mapped onto the simulation grid using routine
113 : ! ModelLev\_Interpolate.
114 : !\\
115 : !\\
116 : ! !INTERFACE:
117 : !
118 0 : SUBROUTINE HCOIO_Read( HcoState, Lct, RC )
119 : !
120 : ! !USES:
121 : !
122 : USE HCO_Ncdf_Mod, ONLY : NC_Open
123 : USE HCO_Ncdf_Mod, ONLY : NC_Close
124 : USE HCO_Ncdf_Mod, ONLY : NC_Read_Var
125 : USE HCO_Ncdf_Mod, ONLY : NC_Read_Arr
126 : USE HCO_Ncdf_Mod, ONLY : NC_Get_Grid_Edges
127 : USE HCO_Ncdf_Mod, ONLY : NC_Get_Sigma_Levels
128 : USE HCO_CHARPAK_MOD, ONLY : TRANLC
129 : USE HCO_Unit_Mod, ONLY : HCO_Unit_Change
130 : USE HCO_Unit_Mod, ONLY : HCO_Unit_ScalCheck
131 : USE HCO_Unit_Mod, ONLY : HCO_IsUnitless
132 : USE HCO_Unit_Mod, ONLY : HCO_IsIndexData
133 : USE HCO_Unit_Mod, ONLY : HCO_UnitTolerance
134 : USE HCO_GeoTools_Mod, ONLY : HCO_ValidateLon
135 : USE HCO_FileData_Mod, ONLY : FileData_ArrCheck
136 : USE HCO_FileData_Mod, ONLY : FileData_ArrInit
137 : USE HCO_FileData_Mod, ONLY : FileData_Cleanup
138 : USE HCOIO_MESSY_MOD, ONLY : HCO_MESSY_REGRID
139 : USE HCO_INTERP_MOD, ONLY : REGRID_MAPA2A
140 : USE HCO_INTERP_MOD, ONLY : ModelLev_Check
141 : USE HCO_CLOCK_MOD, ONLY : HcoClock_Get
142 : USE HCO_DIAGN_MOD, ONLY : Diagn_Update
143 : USE HCO_EXTLIST_MOD, ONLY : HCO_GetOpt
144 : USE HCO_TIDX_MOD, ONLY : tIDx_IsInRange
145 :
146 : include "netcdf.inc"
147 : !
148 : ! !INPUT PARAMETERS:
149 : !
150 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
151 : TYPE(ListCont), POINTER :: Lct ! HEMCO list container
152 : !
153 : ! !INPUT/OUTPUT PARAMETERS:
154 : !
155 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
156 : !
157 : ! !REVISION HISTORY:
158 : ! 13 Mar 2013 - C. Keller - Initial version
159 : ! See https://github.com/geoschem/hemco for complete history
160 : !EOP
161 : !------------------------------------------------------------------------------
162 : !BOC
163 : !
164 : ! !LOCAL VARIABLES:
165 : !
166 : CHARACTER(LEN=255) :: thisUnit, LevUnit, LevName
167 : CHARACTER(LEN=1023) :: MSG, LOC
168 : CHARACTER(LEN=1023) :: srcFile, srcFile2
169 : INTEGER :: NX, NY
170 : INTEGER :: NCRC, Flag, AS
171 : INTEGER :: ncLun, ncLun2
172 : INTEGER :: ierr, v_id
173 : INTEGER :: nlon, nlat, nlev, nTime
174 : INTEGER :: lev1, lev2, dir
175 : INTEGER :: tidx1, tidx2, ncYr, ncMt
176 : INTEGER :: tidx1b, tidx2b, ncYr2, ncMt2
177 : INTEGER :: HcoID
178 : INTEGER :: ArbIdx
179 : INTEGER :: nlatEdge, nlonEdge
180 : INTEGER :: Direction
181 : REAL(hp) :: MW_g
182 : REAL(sp) :: wgt1, wgt2
183 0 : REAL(sp), POINTER :: ncArr(:,:,:,:)
184 0 : REAL(sp), POINTER :: ncArr2(:,:,:,:)
185 0 : REAL(hp), POINTER :: SigEdge(:,:,:)
186 0 : REAL(hp), POINTER :: SigLev (:,:,:)
187 0 : REAL(hp), POINTER :: LonMid (:)
188 0 : REAL(hp), POINTER :: LatMid (:)
189 0 : REAL(hp), POINTER :: LevMid (:)
190 0 : REAL(hp), POINTER :: LonEdge (:)
191 0 : REAL(hp), POINTER :: LatEdge (:)
192 : REAL(hp) :: UnitFactor
193 : LOGICAL :: KeepSpec
194 : LOGICAL :: FOUND
195 : LOGICAL :: IsModelLevel
196 : LOGICAL :: DoReturn
197 : INTEGER :: UnitTolerance
198 : INTEGER :: AreaFlag, TimeFlag
199 : REAL(dp) :: YMDhma, YMDhmb, YMDhm1
200 : REAL(dp) :: oYMDhm1, oYMDhm2
201 : INTEGER :: cYr, cMt, cDy, cHr, Yr1, Yr2
202 : INTEGER :: nYears, iYear
203 : INTEGER :: I, J
204 :
205 : ! Use MESSy regridding routines?
206 : LOGICAL :: UseMESSy
207 :
208 : ! SAVEd scalars
209 : LOGICAL, SAVE :: doPrintWarning = .TRUE.
210 :
211 : !=================================================================
212 : ! HCOIO_READ begins here
213 : !=================================================================
214 0 : LOC = 'HCOIO_READ (HCOIO_READ_STD_MOD.F90)'
215 :
216 : ! Do not try to read a mask file where the mask bounding box limits
217 : ! are given in the srcFile location, as there is no file to read.
218 : ! This should fix https://github.com/geoschem/HEMCO/issues/153.
219 : ! -- Bob Yantosca (12 Jul 2022)
220 0 : IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
221 0 : IF ( .not. Lct%Dct%Dta%ncRead ) RETURN
222 : ENDIF
223 :
224 : ! Enter
225 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
226 0 : IF ( RC /= HCO_SUCCESS ) THEN
227 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
228 0 : RETURN
229 : ENDIF
230 :
231 : ! Initialize pointers
232 0 : ncArr => NULL()
233 0 : ncArr2 => NULL()
234 0 : SigEdge => NULL()
235 0 : SigLev => NULL()
236 0 : LonMid => NULL()
237 0 : LatMid => NULL()
238 0 : LevMid => NULL()
239 0 : LonEdge => NULL()
240 0 : LatEdge => NULL()
241 :
242 : ! Zero local variables for safety's sake
243 0 : dir = 0
244 0 : lev1 = 0
245 0 : lev2 = 0
246 0 : ncYr = 0
247 0 : ncMt = 0
248 0 : ncYr2 = 0
249 0 : ncMt2 = 0
250 0 : nLon = 0
251 0 : nLat = 0
252 0 : nLev = 0
253 0 : nTime = 0
254 0 : tIdx1 = 0
255 0 : tIdx2 = 0
256 0 : tidx1b = 0
257 0 : tidx2b = 0
258 0 : wgt1 = 0.0_sp
259 0 : wgt2 = 0.0_sp
260 :
261 : ! Get unit tolerance set in configuration file
262 0 : UnitTolerance = HCO_UnitTolerance( HcoState%Config )
263 :
264 : ! For convenience, copy horizontal grid dimensions from HEMCO
265 : ! state object
266 0 : NX = HcoState%NX
267 0 : NY = HcoState%NY
268 :
269 : ! Verbose
270 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
271 0 : WRITE(MSG,*) 'Processing container: ', TRIM(Lct%Dct%cName)
272 0 : CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-' )
273 : ENDIF
274 :
275 : ! If the file has cycle flag "E" (e.g. it's a restart file), then we will
276 : ! read it only once and then never again. If the file has already been
277 : ! read on a previous call, then don't call HCOIO_READ. (bmy, 10/4/18)
278 : !
279 : ! Moved this handling from hcoio_dataread_mod, as it is non-MAPL specific
280 : ! (hplin, 4/5/21)
281 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT .and. &
282 0 : Lct%Dct%Dta%UpdtFlag == HCO_UFLAG_ONCE .and. &
283 : Lct%Dct%Dta%isTouched ) THEN
284 :
285 : ! Print a warning message only once
286 0 : IF ( doPrintWarning ) THEN
287 0 : doPrintWarning = .FALSE.
288 : MSG = 'No further attempts will be made to read file: ' // &
289 0 : TRIM( Lct%Dct%Dta%NcFile )
290 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
291 : ENDIF
292 :
293 : ! Return without reading
294 0 : CALL HCO_LEAVE( HcoState%Config%Err, RC )
295 0 : RETURN
296 : ENDIF
297 :
298 : ! ----------------------------------------------------------------
299 : ! Parse source file name. This will replace all tokens ($ROOT,
300 : ! ($YYYY), etc., with valid values.
301 : ! ----------------------------------------------------------------
302 0 : CALL SrcFile_Parse( HcoState, Lct, srcFile, FOUND, RC )
303 0 : IF ( RC /= HCO_SUCCESS ) THEN
304 : MSG = 'Error encountered in routine "SrcFile_Parse", located ' // &
305 0 : 'module src/Core/hcoio_read_std_mod.F90!'
306 0 : CALL HCO_ERROR( MSG, RC )
307 0 : RETURN
308 : ENDIF
309 :
310 : ! Handle found or not in the standard way if HEMCO is in regular run mode.
311 0 : IF ( .NOT. HcoState%Options%isDryRun ) THEN
312 :
313 : !====================================================================
314 : ! HEMCO is in regular simulation mode (not dry-run)!
315 : !====================================================================
316 :
317 : ! If file not found, return w/ error. No error if cycling attribute is
318 : ! select to range. In that case, just make sure that array is empty.
319 0 : IF ( .NOT. FOUND ) THEN
320 0 : IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) .OR. &
321 : ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) ) THEN
322 :
323 : ! If MustFind flag is enabled, return with error if field is not
324 : ! found
325 0 : IF ( Lct%Dct%Dta%MustFind ) THEN
326 : MSG = 'Cannot find file for current simulation time: ' // &
327 : TRIM(srcFile) // ' - Cannot get field ' // &
328 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
329 0 : 'and time (incl. time range flag) in the config. file'
330 0 : CALL HCO_ERROR( MSG, RC )
331 0 : RETURN
332 :
333 : ! If MustFind flag is not enabled, ignore this field and return
334 : ! with a warning.
335 : ELSE
336 0 : CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
337 : MSG = 'No valid file found for current simulation time - data '// &
338 0 : 'will be ignored for time being - ' // TRIM(Lct%Dct%cName)
339 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
340 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
341 0 : RETURN
342 : ENDIF
343 :
344 : ELSE
345 : MSG = 'Cannot find file for current simulation time: ' // &
346 : TRIM(srcFile) // ' - Cannot get field ' // &
347 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
348 0 : 'and time (incl. time range flag) in the config. file'
349 0 : CALL HCO_ERROR( MSG, RC )
350 0 : RETURN
351 : ENDIF
352 : ENDIF
353 :
354 : ELSE
355 :
356 : !====================================================================
357 : ! HEMCO is in a "dry-run" mode!
358 : !====================================================================
359 :
360 : ! Simulate file read buffer
361 0 : IF ( TRIM(HcoState%ReadLists%FileInArchive) == TRIM(srcFile) ) THEN
362 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
363 0 : RETURN
364 : ENDIF
365 :
366 : ! If file exists, print the result. If NOT, then handle accordingly
367 : ! But NEVER error out (HCO_ERROR), as we want to get a list of all
368 : ! files. (hplin, 11/2/19)
369 0 : IF ( .NOT. FOUND ) THEN
370 0 : IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) .OR. &
371 : ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) ) THEN
372 :
373 : ! If MustFind flag is enabled, return with error if field is not
374 : ! found
375 0 : IF ( Lct%Dct%Dta%MustFind ) THEN
376 : MSG = 'Cannot find file for current simulation time: ' // &
377 : TRIM(srcFile) // ' - Cannot get field ' // &
378 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
379 0 : 'and time (incl. time range flag) in the config. file'
380 0 : CALL HCO_Warning( HcoState%Config%Err, MSG, RC )
381 :
382 : ! Write a msg to stdout (NOT FOUND)
383 0 : WRITE( 6, 300 ) TRIM( srcFile )
384 : 300 FORMAT( 'HEMCO: REQUIRED FILE NOT FOUND ', a )
385 :
386 : ! If MustFind flag is not enabled, ignore this field and return
387 : ! with a warning.
388 : ELSE
389 0 : CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
390 : MSG = 'No valid file found for current simulation time - '// &
391 : 'data will be ignored for time being - ' // &
392 0 : TRIM(Lct%Dct%cName)
393 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
394 :
395 : ! Write a msg to stdout (OPTIONAL)
396 0 : WRITE( 6, 310 ) TRIM( srcFile )
397 : 310 FORMAT( 'HEMCO: OPTIONAL FILE NOT FOUND ', a )
398 :
399 : ENDIF
400 :
401 : ! Not range or exact
402 : ELSE
403 : MSG = 'Cannot find file for current simulation time: ' // &
404 : TRIM(srcFile) // ' - Cannot get field ' // &
405 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
406 0 : 'and time (incl. time range flag) in the config. file'
407 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
408 :
409 : ! Write a msg to stdout (NOT FOUND)
410 0 : WRITE( 6, 300 ) TRIM(srcFile)
411 :
412 : ENDIF
413 : ELSE
414 :
415 : ! Write a mesage to stdout (HEMCO: Opening...)
416 0 : WRITE( 6, 100 ) TRIM( srcFile )
417 :
418 : ENDIF
419 :
420 : ! It is safe to leave now, we do not need to handle opening the file.
421 : ! This may be changed in the future if the "dry-run" mode requires
422 : ! a check of the file contents... this may be beyond the scope for now.
423 :
424 : ! Simulate the "reading" in netCDF to prevent duplicate entries
425 : ! in the log
426 0 : HcoState%ReadLists%FileInArchive = TRIM(srcFile)
427 :
428 : ! Skip further processing
429 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
430 0 : RETURN
431 : ENDIF ! End of dry-run mode else clause
432 :
433 : ! ----------------------------------------------------------------
434 : ! Open netCDF
435 : ! ----------------------------------------------------------------
436 :
437 : ! Check if file is already in buffer. In that case use existing
438 : ! open stream. Otherwise open new file. At any given time there
439 : ! can only be one file in buffer.
440 0 : ncLun = -1
441 0 : IF ( HcoState%ReadLists%FileLun > 0 ) THEN
442 0 : IF ( TRIM(HcoState%ReadLists%FileInArchive) == TRIM(srcFile) ) THEN
443 0 : ncLun = HcoState%ReadLists%FileLun
444 : ELSE
445 0 : CALL NC_CLOSE ( HcoState%ReadLists%FileLun )
446 0 : HcoState%ReadLists%FileLun = -1
447 : ENDIF
448 : ENDIF
449 :
450 : ! To read from existing stream:
451 0 : IF ( ncLun > 0 ) THEN
452 :
453 : ! Verbose mode
454 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
455 0 : WRITE(MSG,*) 'Reading from existing stream: ', TRIM(srcFile)
456 0 : CALL HCO_MSG( HcoState%Config%Err, MSG )
457 : ENDIF
458 :
459 : ! To open a new file:
460 : ELSE
461 0 : CALL NC_OPEN ( TRIM(srcFile), ncLun )
462 :
463 : ! Verbose mode
464 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
465 0 : WRITE(MSG,*) 'Opening file: ', TRIM(srcFile)
466 0 : CALL HCO_MSG( HcoState%Config%Err, MSG )
467 : ENDIF
468 :
469 : ! Also write to standard output
470 0 : WRITE( 6, 100 ) TRIM( srcFile )
471 : 100 FORMAT( 'HEMCO: Opening ', a )
472 :
473 : ! This is now the file in archive
474 0 : HcoState%ReadLists%FileInArchive = TRIM(srcFile)
475 0 : HcoState%ReadLists%FileLun = ncLun
476 : ENDIF
477 :
478 : ! ----------------------------------------------------------------
479 : ! Extract time slice information
480 : ! This determines the lower and upper time slice index (tidx1
481 : ! and tidx2) to be read based upon the time slice information
482 : ! extracted from the file and the time stamp settings set in the
483 : ! HEMCO configuration file. Multiple time slices are only selected
484 : ! for weekdaily data or for 'autodetected' hourly data (using the
485 : ! wildcard character in the configuration file time attribute) or
486 : ! if data shall be interpolated between two (consecutive) time
487 : ! slices. The weights to be assigned to those two time slices is
488 : ! also calculated in GET_TIMEIDX and returned as variables wgt1
489 : ! and wgt2, respectively.
490 : ! ----------------------------------------------------------------
491 : CALL GET_TIMEIDX ( HcoState, Lct, &
492 : ncLun, tidx1, tidx2, &
493 : wgt1, wgt2, oYMDhm1, &
494 0 : YMDhma, YMDhm1, RC )
495 0 : IF ( RC /= HCO_SUCCESS ) THEN
496 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
497 0 : RETURN
498 : ENDIF
499 :
500 : !-----------------------------------------------------------------
501 : ! Check for negative tidx1. tidx1 can still be negative if:
502 : ! (a) CycleFlag is set to range and the current simulation
503 : ! time is outside of the data time range. In this case, we
504 : ! prompt a warning and make sure that there is no data
505 : ! associated with this FileData container.
506 : ! (b) CycleFlag is set to exact and none of the data time
507 : ! stamps matches the current simulation time exactly. Return
508 : ! with error!
509 : !-----------------------------------------------------------------
510 0 : IF ( tidx1 < 0 ) THEN
511 0 : DoReturn = .FALSE.
512 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_CYCLE ) THEN
513 0 : MSG = 'Invalid time index in ' // TRIM(srcFile)
514 0 : CALL HCO_ERROR( MSG, RC )
515 : DoReturn = .TRUE.
516 0 : ELSEIF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) .OR. &
517 : ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) ) THEN
518 0 : IF ( Lct%Dct%Dta%MustFind ) THEN
519 : MSG = 'Cannot find field with valid time stamp in ' // &
520 : TRIM(srcFile) // ' - Cannot get field ' // &
521 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
522 0 : 'and time (incl. time range flag) in the config. file'
523 0 : CALL HCO_ERROR( MSG, RC )
524 : DoReturn = .TRUE.
525 : ELSE
526 0 : CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE.)
527 : MSG = 'Simulation time is outside of time range provided for '//&
528 0 : TRIM(Lct%Dct%cName) // ' - field is ignored for the time being!'
529 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
530 0 : DoReturn = .TRUE.
531 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
532 : ENDIF
533 : ENDIF
534 :
535 : ! Eventually return here
536 : IF ( DoReturn ) THEN
537 0 : RETURN
538 : ENDIF
539 : ENDIF
540 :
541 : ! ----------------------------------------------------------------
542 : ! Check if variable is in file
543 : ! ----------------------------------------------------------------
544 0 : ierr = Nf_Inq_Varid( ncLun, Lct%Dct%Dta%ncPara, v_id )
545 0 : IF ( ierr /= NF_NOERR ) THEN
546 :
547 : ! If MustFind flag is enabled, return with error if field is not
548 : ! found
549 0 : IF ( Lct%Dct%Dta%MustFind ) THEN
550 : MSG = 'Cannot find field ' // TRIM(Lct%Dct%cName) // &
551 0 : '. Please check variable name in the config. file'
552 0 : CALL HCO_ERROR( MSG, RC )
553 0 : RETURN
554 :
555 : ! If MustFind flag is not enabled, ignore this field and return
556 : ! with a warning.
557 : ELSE
558 0 : CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
559 : MSG = 'Cannot find field ' // TRIM(Lct%Dct%cName) // &
560 0 : '. Will be ignored for time being.'
561 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
562 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
563 0 : RETURN
564 : ENDIF
565 : ENDIF
566 :
567 : ! ----------------------------------------------------------------
568 : ! Read grid
569 : ! ----------------------------------------------------------------
570 :
571 : ! Extract longitude midpoints
572 0 : CALL NC_READ_VAR ( ncLun, 'lon', nlon, thisUnit, LonMid, NCRC )
573 0 : IF ( NCRC /= 0 ) THEN
574 0 : CALL HCO_ERROR( 'NC_READ_VAR: lon', RC )
575 0 : RETURN
576 : ENDIF
577 :
578 0 : IF ( nlon == 0 ) THEN
579 0 : CALL NC_READ_VAR ( ncLun, 'longitude', nlon, thisUnit, LonMid, NCRC )
580 : ENDIF
581 0 : IF ( NCRC /= 0 ) THEN
582 0 : CALL HCO_ERROR( 'NC_READ_VAR: longitude', RC )
583 0 : RETURN
584 : ENDIF
585 :
586 0 : IF ( nlon == 0 ) THEN
587 0 : CALL NC_READ_VAR ( ncLun, 'Longitude', nlon, thisUnit, LonMid, NCRC )
588 : ENDIF
589 0 : IF ( NCRC /= 0 ) THEN
590 0 : CALL HCO_ERROR( 'NC_READ_LON: Longitude', RC )
591 0 : RETURN
592 : ENDIF
593 :
594 0 : IF ( nlon == 0 ) THEN
595 : MSG = 'Cannot find longitude variable in ' // TRIM(srcFile) // &
596 0 : ' - Must be one of `lon`, `longitude`, `Longitude`'
597 0 : CALL HCO_ERROR( MSG, RC )
598 0 : RETURN
599 : ENDIF
600 :
601 : ! Unit must be degrees_east
602 0 : CALL TRANLC( thisUnit)
603 0 : IF ( INDEX( thisUnit, 'degrees_east' ) == 0 ) THEN
604 : MSG = 'illegal longitude unit in ' // TRIM(srcFile) // &
605 0 : ' - Must be `degrees_east`.'
606 0 : CALL HCO_ERROR( MSG, RC )
607 0 : RETURN
608 : ENDIF
609 :
610 : ! Make sure longitude is steadily increasing.
611 0 : CALL HCO_ValidateLon( HcoState, nlon, LonMid, RC )
612 0 : IF ( RC /= HCO_SUCCESS ) THEN
613 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
614 0 : RETURN
615 : ENDIF
616 :
617 : ! Extract latitude midpoints
618 0 : CALL NC_READ_VAR ( ncLun, 'lat', nlat, thisUnit, LatMid, NCRC )
619 0 : IF ( NCRC /= 0 ) THEN
620 0 : CALL HCO_ERROR( 'NC_READ_LON: lat', RC )
621 0 : RETURN
622 : ENDIF
623 :
624 0 : IF ( nlat == 0 ) THEN
625 0 : CALL NC_READ_VAR ( ncLun, 'latitude', nlat, thisUnit, LatMid, NCRC )
626 : ENDIF
627 0 : IF ( NCRC /= 0 ) THEN
628 0 : CALL HCO_ERROR( 'NC_READ_LON: latitude', RC )
629 0 : RETURN
630 : ENDIF
631 :
632 0 : IF ( nlat == 0 ) THEN
633 0 : CALL NC_READ_VAR ( ncLun, 'Latitude', nlat, thisUnit, LatMid, NCRC )
634 : ENDIF
635 0 : IF ( NCRC /= 0 ) THEN
636 0 : CALL HCO_ERROR( 'NC_READ_LON: Latitude', RC )
637 0 : RETURN
638 : ENDIF
639 :
640 0 : IF ( nlat == 0 ) THEN
641 : MSG = 'Cannot find latitude variable in ' // TRIM(srcFile) // &
642 0 : ' - Must be one of `lat`, `latitude`, `Latitude`'
643 0 : CALL HCO_ERROR( MSG, RC )
644 0 : RETURN
645 : ENDIF
646 :
647 : ! Unit must be degrees_north
648 0 : CALL TRANLC( thisUnit)
649 0 : IF ( INDEX( thisUnit, 'degrees_north' ) == 0 ) THEN
650 : MSG = 'illegal latitude unit in ' // TRIM(srcFile) // &
651 0 : ' - Must be `degrees_north`.'
652 0 : CALL HCO_ERROR( MSG, RC )
653 0 : RETURN
654 : ENDIF
655 :
656 : ! Get level index if we are dealing with 3D data
657 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
658 :
659 : ! Try to extract level midpoints
660 0 : LevName = 'lev'
661 0 : CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
662 0 : IF ( NCRC /= 0 ) THEN
663 0 : CALL HCO_ERROR( 'NC_READ_VAR: lev', RC )
664 0 : RETURN
665 : ENDIF
666 0 : IF ( nlev == 0 ) THEN
667 0 : LevName = 'height'
668 0 : CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
669 0 : IF ( NCRC /= 0 ) THEN
670 0 : CALL HCO_ERROR( 'NC_READ_VAR: height', RC )
671 0 : RETURN
672 : ENDIF
673 : ENDIF
674 0 : IF ( nlev == 0 ) THEN
675 0 : LevName = 'level'
676 0 : CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
677 0 : IF ( NCRC /= 0 ) THEN
678 0 : CALL HCO_ERROR( 'NC_READ_VAR: level', RC )
679 0 : RETURN
680 : ENDIF
681 : ENDIF
682 :
683 : ! Error check
684 0 : IF ( nlev == 0 ) THEN
685 : MSG = 'Cannot find vertical coordinate variable in ' // &
686 0 : TRIM(SrcFile) // ' - Must be one of `lev`, `level`, `height`.'
687 0 : CALL HCO_ERROR( MSG, RC )
688 0 : RETURN
689 : ENDIF
690 :
691 : ! Are these model levels? This will only return true if
692 : ! (1) the variable is on 72/73 levels and you are going to 47
693 : ! levels, (2) if you are on 102/103 levels and you are going
694 : ! to 74 levels, (3) if you are on 47/48 levels and you are
695 : ! going to 72 levels. Otherwise, use MESSy (nbalasus, 8/24/2023).
696 0 : IF ( Lct%Dct%Dta%Levels == 0 ) THEN
697 :
698 : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
699 :
700 : ! In WRF/CESM, IsModelLevel has a different meaning of "GEOS-Chem levels"
701 : ! because the models in WRF and CESM are user-defined and thus fixed input
702 : ! files would never be on the model level. In this case, a check is added
703 : ! in order to match the file with known GEOS-Chem levels, and if so, the
704 : ! data will be handled later in this file accordingly to be vertically
705 : ! regridded to the runtime model levels using MESSy.
706 : ! This fixes a regression from the vertical regridding fixes in 3.7.1.
707 : !
708 : ! The meaning of "is model levels" in WRF and CESM are different.
709 : ! Model levels can be changed and thus data is never on the model level.
710 : ! In this case, IsModelLevel means that the data is on standard
711 : ! GEOS-Chem levels, and if so, the data will be handled accordingly
712 : ! using a hard-coded set of GEOS-Chem levels to be interpolated using MESSy.
713 : ! (hplin, 10/15/23)
714 0 : IF ( TRIM(LevUnit) == "level" .or. TRIM(LevUnit) == "GEOS-Chem level" ) THEN
715 : ! the below check will be obsolete and is unmaintainable, but would be consistent with ModelLev_Check.
716 : ! it is more robust to check for the explicit intention of LevUnit
717 : ! nlev == 47 .or. nlev == 48 .or. nlev == 36 .or. nlev == 72 .or. nlev == 73 ) THEN
718 0 : IsModelLevel = .true.
719 : ENDIF
720 :
721 : #else
722 :
723 : CALL ModelLev_Check( HcoState, nlev, IsModelLevel, RC )
724 : IF ( RC /= HCO_SUCCESS ) THEN
725 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
726 : RETURN
727 : ENDIF
728 :
729 : #endif
730 :
731 : ! Set level indeces to be read
732 0 : lev1 = 1
733 0 : lev2 = nlev
734 :
735 : ! If levels are explicitly given:
736 : ELSE
737 :
738 : ! Number of levels to be read must be smaller or equal to total
739 : ! number of available levels
740 0 : IF ( ABS(Lct%Dct%Dta%Levels) > nlev ) THEN
741 0 : WRITE(MSG,*) Lct%Dct%Dta%Levels, ' levels requested but file ', &
742 0 : 'has only ', nlev, ' levels: ', TRIM(Lct%Dct%cName)
743 0 : CALL HCO_ERROR( MSG, RC )
744 0 : RETURN
745 : ENDIF
746 :
747 : ! Set levels to be read
748 0 : IF ( Lct%Dct%Dta%Levels > 0 ) THEN
749 0 : lev1 = 1
750 0 : lev2 = Lct%Dct%Dta%Levels
751 :
752 : ! Reverse axis!
753 : ELSE
754 0 : lev1 = nlev
755 0 : lev2 = nlev + Lct%Dct%Dta%Levels + 1
756 : ENDIF
757 :
758 : ! Use MESSy regridding
759 0 : IsModelLevel = .FALSE.
760 :
761 : ENDIF
762 :
763 : ! Verbose
764 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
765 0 : WRITE(MSG,*) 'Will read vertical levels ', lev1, ' to ', lev2
766 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
767 : ENDIF
768 :
769 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) .AND. IsModelLevel ) THEN
770 0 : WRITE(MSG,*) 'Data is assumed to already be on the model level grid'
771 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
772 : ENDIF
773 :
774 : ! For 2D data, set lev1 and lev2 to zero. This will ignore
775 : ! the level dimension in the netCDF reading call that follows.
776 : ELSE
777 0 : nlev = 0
778 0 : lev1 = 0
779 0 : lev2 = 0
780 0 : IsModelLevel = .FALSE.
781 : ENDIF
782 :
783 : ! ----------------------------------------------------------------
784 : ! Check for arbitrary additional dimension. Will return -1 if not
785 : ! set.
786 : ! ----------------------------------------------------------------
787 0 : CALL GetArbDimIndex( HcoState, ncLun, Lct, ArbIdx, RC )
788 0 : IF ( RC /= HCO_SUCCESS ) THEN
789 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
790 0 : RETURN
791 : ENDIF
792 :
793 : ! ----------------------------------------------------------------
794 : ! Read data
795 : ! ----------------------------------------------------------------
796 :
797 : ! Verbose mode
798 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
799 0 : WRITE(MSG,*) 'Reading variable ', TRIM(Lct%Dct%Dta%ncPara)
800 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
801 : ENDIF
802 :
803 : CALL NC_READ_ARR( fID = ncLun, &
804 : ncVar = Lct%Dct%Dta%ncPara, &
805 : lon1 = 1, &
806 : lon2 = nlon, &
807 : lat1 = 1, &
808 : lat2 = nlat, &
809 : lev1 = lev1, &
810 : lev2 = lev2, &
811 : time1 = tidx1, &
812 : time2 = tidx2, &
813 : ncArr = ncArr, &
814 : varUnit = thisUnit, &
815 : wgt1 = wgt1, &
816 : wgt2 = wgt2, &
817 : MissVal = HCO_MISSVAL, &
818 : ArbIdx = ArbIdx, &
819 0 : RC = NCRC )
820 :
821 0 : IF ( NCRC /= 0 ) THEN
822 0 : CALL HCO_ERROR( 'NC_READ_ARRAY', RC )
823 0 : RETURN
824 : ENDIF
825 :
826 : ! Check for missing values: set base emissions and masks to 0, and
827 : ! scale factors to 1. This will make sure that these entries will
828 : ! be ignored.
829 : !!! CALL CheckMissVal ( Lct, ncArr )
830 :
831 : !-----------------------------------------------------------------
832 : ! Eventually do interpolation between files. This is a pretty
833 : ! crude implementation for data interpolation between different
834 : ! files. It is only applied to data that is marked as interpolated
835 : ! data and if no appropriate interpolation date could be found in
836 : ! the first file. This will only be the case if the preferred date-
837 : ! time is outside the file range.
838 : !-----------------------------------------------------------------
839 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .AND. wgt1 < 0.0_sp ) THEN
840 :
841 : ! No need to read another file if the previous file had exactly the
842 : ! time stamp we were looking for.
843 0 : IF ( oYMDhm1 == YMDhma ) THEN
844 0 : FOUND = .FALSE.
845 : ELSE
846 : ! Check if there exists another file for a future/previous date
847 : ! oYMDhm1 is the originally preferred date, YMDhma is the date read
848 : ! from the file. If YMDhma is in the future, we need to look for a
849 : ! file in the past. Else, need to look for a file in the future.
850 0 : IF ( oYMDhm1 < YMDhma ) THEN
851 0 : Direction = -1
852 : ELSE
853 0 : Direction = +1
854 : ENDIF
855 : CALL SrcFile_Parse ( HcoState, Lct, srcFile2, &
856 0 : FOUND, RC, Direction = Direction )
857 0 : IF ( RC /= HCO_SUCCESS ) THEN
858 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
859 0 : RETURN
860 : ENDIF
861 : ENDIF
862 :
863 : ! If found, read data. Assume that all meta-data is the same.
864 0 : IF ( FOUND ) THEN
865 :
866 : ! Open file
867 0 : CALL NC_OPEN ( TRIM(srcFile2), ncLun2 )
868 :
869 : ! Define time stamp to be read. Use this call only
870 : ! to get the datetime of the first time slice (YMDhm1).
871 : ! All other values will be ignored and reset below.
872 : CALL GET_TIMEIDX ( HcoState, Lct, &
873 : ncLun2, tidx1, tidx2, &
874 : wgt1, wgt2, oYMDhm2, &
875 0 : YMDhmb, YMDhm1, RC )
876 0 : IF ( RC /= HCO_SUCCESS ) THEN
877 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
878 0 : RETURN
879 : ENDIF
880 :
881 : ! Always read first time slice
882 0 : tidx1 = 1
883 0 : tidx2 = 1
884 0 : wgt1 = -1.0_sp
885 0 : wgt2 = -1.0_sp
886 :
887 : ! Read data and write into array ncArr2
888 : CALL NC_READ_ARR( fID = ncLun2, &
889 : ncVar = Lct%Dct%Dta%ncPara, &
890 : lon1 = 1, &
891 : lon2 = nlon, &
892 : lat1 = 1, &
893 : lat2 = nlat, &
894 : lev1 = lev1, &
895 : lev2 = lev2, &
896 : time1 = tidx1, &
897 : time2 = tidx2, &
898 : ncArr = ncArr2, &
899 : varUnit = thisUnit, &
900 : wgt1 = wgt1, &
901 : wgt2 = wgt2, &
902 : MissVal = HCO_MISSVAL, &
903 : ArbIdx = ArbIdx, &
904 0 : RC = NCRC )
905 0 : IF ( NCRC /= 0 ) THEN
906 0 : CALL HCO_ERROR( 'NC_READ_ARRAY (2)', RC )
907 0 : RETURN
908 : ENDIF
909 :
910 : ! Eventually fissing values
911 : !!! CALL CheckMissVal ( Lct, ncArr2 )
912 :
913 : ! Calculate weights to be applied to ncArr2 and ncArr1. These
914 : ! weights are calculated based on the originally preferred
915 : ! datetime oYMDh1 and the selected datetime of file 1 (YMDhma)
916 : ! and file 2 (YMDhm1)
917 : ! If date on file 1 < date on file 2:
918 0 : IF ( YMDhma < YMDhm1 ) THEN
919 0 : CALL GetWeights ( YMDhma, YMDhm1, oYMDhm1, wgt1, wgt2 )
920 : ! If date on file 1 > date on file 2:
921 0 : ELSEIF ( YMDhma > YMDhm1 ) THEN
922 0 : CALL GetWeights ( YMDhm1, YMDhma, oYMDhm1, wgt2, wgt1 )
923 : ! If both datetimes are for some reason the same (this should
924 : ! not happen!)
925 : ELSE
926 0 : wgt1 = 0.5_sp
927 0 : wgt2 = 0.5_sp
928 : ENDIF
929 :
930 : ! Apply weights
931 0 : ncArr = (wgt1 * ncArr) + (wgt2 * ncArr2)
932 :
933 : ! Verbose
934 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
935 0 : MSG = 'Interpolated data between two files:'
936 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
937 0 : MSG = '- File 1: ' // TRIM(srcFile)
938 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
939 0 : WRITE(MSG,*) ' Time stamp used: ', YMDhma
940 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
941 0 : WRITE(MSG,*) ' Applied weight: ', wgt1
942 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
943 0 : MSG = '- File 2: ' // TRIM(srcFile2)
944 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
945 0 : WRITE(MSG,*) ' Time stamp used: ', YMDhm1
946 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
947 0 : WRITE(MSG,*) ' Applied weight: ', wgt2
948 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
949 : ENDIF
950 :
951 : ! Cleanup
952 0 : IF ( ASSOCIATED(ncArr2) ) DEALLOCATE(ncArr2)
953 :
954 : ! Close file
955 0 : CALL NC_CLOSE ( ncLun2 )
956 : ENDIF !FOUND
957 :
958 : !-----------------------------------------------------------------
959 : ! Eventually calculate averages. Currently, averages are only
960 : ! calculated on the year dimension, e.g. over years.
961 : !-----------------------------------------------------------------
962 0 : ELSEIF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_AVERG .OR. &
963 : Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGEAVG ) THEN
964 :
965 : ! cYr is the current simulation year
966 : CALL HcoClock_Get( HcoState%Clock, cYYYY=cYr, cMM=cMt, cDD=cDy, &
967 0 : cH=cHr, RC=RC )
968 0 : IF ( RC /= HCO_SUCCESS ) THEN
969 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
970 0 : RETURN
971 : ENDIF
972 :
973 : ! Determine year range to be read:
974 : ! By default, we would like to average between the year range given
975 : ! in the time attribute
976 0 : Yr1 = Lct%Dct%Dta%ncYrs(1)
977 0 : Yr2 = Lct%Dct%Dta%ncYrs(2)
978 :
979 : ! If averaging shall only be performed if outside the given
980 : ! range, check if current simulation date is within the range
981 : ! provied in the configuration file. If so, set year range to
982 : ! be read to current year only.
983 0 : IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGEAVG ) ) THEN
984 0 : IF ( tIDx_IsInRange(Lct,cYr,cMt,cDy,cHr) ) THEN
985 0 : Yr1 = cYr
986 0 : Yr2 = cYr
987 : ENDIF
988 : ENDIF
989 :
990 : ! Total number of years to be read
991 0 : nYears = Yr2 - Yr1 + 1
992 :
993 : ! Read and add annual data if there is more than one year to be
994 : ! used.
995 0 : IF ( nYears > 1 ) THEN
996 :
997 : ! Cleanup ncArr. This is refilled again
998 0 : ncArr = 0.0_sp
999 :
1000 0 : DO iYear = Yr1, Yr2
1001 :
1002 : ! Get file name for this year
1003 : CALL SrcFile_Parse ( HcoState, Lct, srcFile2, &
1004 0 : FOUND, RC, Year=iYear )
1005 0 : IF ( RC /= HCO_SUCCESS ) THEN
1006 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
1007 0 : RETURN
1008 : ENDIF
1009 :
1010 : ! If found, read data. Assume that all meta-data is the same.
1011 0 : IF ( .NOT. FOUND ) THEN
1012 0 : WRITE(MSG,*) 'Cannot find file for year ', iYear, ' - needed ', &
1013 0 : 'to perform time-averaging on file ', TRIM(Lct%Dct%Dta%ncFile)
1014 0 : CALL HCO_ERROR( MSG, RC )
1015 0 : RETURN
1016 : ENDIF
1017 :
1018 : ! Open file
1019 0 : CALL NC_OPEN ( TRIM(srcFile2), ncLun2 )
1020 :
1021 : ! Define time stamp to be read.
1022 : CALL GET_TIMEIDX ( HcoState, Lct, &
1023 : ncLun2, tidx1, tidx2, &
1024 : wgt1, wgt2, oYMDhm2, &
1025 : YMDhmb, YMDhm1, RC, &
1026 0 : Year=iYear )
1027 0 : IF ( RC /= HCO_SUCCESS ) THEN
1028 0 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
1029 0 : RETURN
1030 : ENDIF
1031 :
1032 : ! Do not perform weights
1033 0 : wgt1 = -1.0_sp
1034 0 : wgt2 = -1.0_sp
1035 :
1036 : ! Read data and write into array ncArr2
1037 : CALL NC_READ_ARR( fID = ncLun2, &
1038 : ncVar = Lct%Dct%Dta%ncPara, &
1039 : lon1 = 1, &
1040 : lon2 = nlon, &
1041 : lat1 = 1, &
1042 : lat2 = nlat, &
1043 : lev1 = lev1, &
1044 : lev2 = lev2, &
1045 : time1 = tidx1, &
1046 : time2 = tidx2, &
1047 : ncArr = ncArr2, &
1048 : varUnit = thisUnit, &
1049 : wgt1 = wgt1, &
1050 : wgt2 = wgt2, &
1051 : MissVal = HCO_MISSVAL, &
1052 : ArbIdx = ArbIdx, &
1053 0 : RC = NCRC )
1054 0 : IF ( NCRC /= 0 ) THEN
1055 0 : CALL HCO_ERROR( 'NC_READ_ARRAY (3)', RC )
1056 0 : RETURN
1057 : ENDIF
1058 :
1059 : ! Eventually fissing values
1060 : !!! CALL CheckMissVal ( Lct, ncArr2 )
1061 :
1062 : ! Add all values to ncArr
1063 0 : ncArr = ncArr + ncArr2
1064 :
1065 : ! Cleanup
1066 0 : IF ( ASSOCIATED(ncArr2) ) DEALLOCATE(ncArr2)
1067 :
1068 : ! Close file
1069 0 : CALL NC_CLOSE ( ncLun2 )
1070 :
1071 : ENDDO !iYear
1072 :
1073 : ! Now calculate average
1074 0 : ncArr = ncArr / REAL(nYears,sp)
1075 :
1076 : ! Verbose
1077 0 : IF ( HcoState%amIRoot .AND. HCO_IsVerb( HcoState%Config%Err ) ) THEN
1078 0 : WRITE(MSG,110) TRIM(Lct%Dct%cName), Yr1, Yr2
1079 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1080 : ENDIF
1081 : 110 FORMAT( 'Field ', a, ': Average data over years ', I4.4, ' to ', I4.4 )
1082 :
1083 : ENDIF !nYears>1
1084 : ENDIF !Averaging
1085 :
1086 : !-----------------------------------------------------------------
1087 : ! Convert to HEMCO units
1088 : ! HEMCO data are all in kg/m2/s for fluxes and kg/m3 for
1089 : ! concentrations. Unit conversion is performed based on the
1090 : ! unit on the input file and the srcUnit attribute given in the
1091 : ! configuration file. By default, HEMCO will attempt to convert
1092 : ! the units found in the input file to the standard quantities
1093 : ! for mass (kg), area (m2 or m3), and time (s). For instance,
1094 : ! g/cm2/hr will be converted to kg/m2/s. The exceptions to this
1095 : ! rule are:
1096 : ! 1. If srcUnit is set to '1', the input data are expected to
1097 : ! be unitless. If the units string on the input file is none
1098 : ! of the units recognized by HEMCO as unitless, an error is
1099 : ! returned if the unit tolerance setting is set to zero, or
1100 : ! a warning is prompted if unit tolerance is greater than zero.
1101 : ! 2. If srcUnit is set to 'count', no unit conversion is performed
1102 : ! and data will be treated as 'index' data, e.g. regridding will
1103 : ! preserve the absolute values.
1104 : !-----------------------------------------------------------------
1105 :
1106 : ! If OrigUnit is set to wildcard character: use unit from source file
1107 0 : IF ( TRIM(Lct%Dct%Dta%OrigUnit) == &
1108 : TRIM(HCO_GetOpt(HcoState%Config%ExtList,'Wildcard')) ) THEN
1109 0 : Lct%Dct%Dta%OrigUnit = TRIM(thisUnit)
1110 : ENDIF
1111 :
1112 : ! If OrigUnit is set to '1' or to 'count', perform no unit
1113 : ! conversion.
1114 0 : IF ( HCO_IsUnitLess(Lct%Dct%Dta%OrigUnit) .OR. &
1115 : HCO_IsIndexData(Lct%Dct%Dta%OrigUnit) ) THEN
1116 :
1117 : ! Check if file unit is also unitless. This will return 0 for
1118 : ! unitless, 1 for HEMCO emission unit, 2 for HEMCO conc. unit,
1119 : ! -1 otherwise.
1120 0 : Flag = HCO_UNIT_SCALCHECK( thisUnit )
1121 :
1122 : ! Return with error if: (1) thisUnit is recognized as HEMCO unit and
1123 : ! unit tolerance is set to zero; (2) thisUnit is neither unitless nor
1124 : ! a HEMCO unit and unit tolerance is set to zero or one.
1125 : ! The unit tolerance is defined in the configuration file.
1126 0 : IF ( Flag /= 0 .AND. UnitTolerance == 0 ) THEN
1127 : MSG = 'Illegal unit: ' // TRIM(thisUnit) // '. File: ' // &
1128 0 : TRIM(srcFile)
1129 0 : CALL HCO_ERROR( MSG, RC )
1130 0 : RETURN
1131 : ENDIF
1132 :
1133 : ! Prompt a warning if thisUnit is not recognized as unitless.
1134 0 : IF ( Flag /= 0 ) THEN
1135 : MSG = 'Data is treated as unitless, but file attribute suggests ' // &
1136 0 : 'it is not: ' // TRIM(thisUnit) // '. File: ' // TRIM(srcFile)
1137 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
1138 : ENDIF
1139 :
1140 : ! Verbose mode
1141 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1142 0 : WRITE(MSG,*) 'Based on srcUnit attribute (', TRIM(Lct%Dct%Dta%OrigUnit), &
1143 0 : '), no unit conversion is performed.'
1144 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1145 : ENDIF
1146 :
1147 : ! Convert to HEMCO units in all other cases.
1148 : ELSE
1149 :
1150 : ! For zero unit tolerance, make sure that thisUnit matches
1151 : ! with unit set in configuration file. For higher unit
1152 : ! tolerances, prompt a level 3 warning.
1153 0 : IF ( TRIM(Lct%Dct%Dta%OrigUnit) /= TRIM(thisUnit) ) THEN
1154 : MSG = 'File units do not match: ' // TRIM(thisUnit) // &
1155 : ' vs. ' // TRIM(Lct%Dct%Dta%OrigUnit) // &
1156 0 : '. File: ' // TRIM(srcFile)
1157 :
1158 0 : IF ( UnitTolerance == 0 ) THEN
1159 0 : CALL HCO_ERROR( MSG, RC )
1160 0 : RETURN
1161 : ELSE
1162 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
1163 : ENDIF
1164 : ENDIF
1165 :
1166 : ! Shadow species properties needed for unit conversion
1167 0 : HcoID = Lct%Dct%HcoID
1168 0 : IF ( HcoID > 0 ) THEN
1169 0 : MW_g = HcoState%Spc(HcoID)%MW_g
1170 : ELSE
1171 0 : MW_g = -999.0_hp
1172 : ENDIF
1173 :
1174 : ! Now convert to HEMCO units. This attempts to convert mass,
1175 : ! area/volume and time to HEMCO standards (kg, m2/m3, s).
1176 0 : ncYr = FLOOR( MOD( oYMDhm1, 1.0e12_dp ) / 1.0e8_dp )
1177 0 : ncMt = FLOOR( MOD( oYMDhm1, 1.0e8_dp ) / 1.0e6_dp )
1178 :
1179 0 : IF ( ncYr == 0 ) THEN
1180 0 : CALL HcoClock_Get( HcoState%Clock, cYYYY = ncYr, RC=RC )
1181 0 : IF ( RC /= HCO_SUCCESS ) THEN
1182 0 : CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
1183 0 : RETURN
1184 : ENDIF
1185 : ENDIF
1186 0 : IF ( ncMt == 0 ) THEN
1187 0 : CALL HcoClock_Get( HcoState%Clock, cMM = ncMt, RC=RC )
1188 0 : IF ( RC /= HCO_SUCCESS ) THEN
1189 0 : CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
1190 0 : RETURN
1191 : ENDIF
1192 : ENDIF
1193 :
1194 : ! Verbose mode
1195 0 : IF ( HcoState%amIRoot .and. HCO_IsVerb( HcoState%Config%Err ) ) THEN
1196 0 : WRITE(MSG,*) 'Unit conversion settings: '
1197 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1198 0 : WRITE(MSG,*) '- Year, month : ', ncYr, ncMt
1199 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1200 : ENDIF
1201 :
1202 : CALL HCO_UNIT_CHANGE( &
1203 : HcoConfig = HcoState%Config, &
1204 : Array = ncArr, &
1205 : Units = thisUnit, &
1206 : MW = MW_g, &
1207 : YYYY = ncYr, &
1208 : MM = ncMt, &
1209 : AreaFlag = AreaFlag, &
1210 : TimeFlag = TimeFlag, &
1211 : FACT = UnitFactor, &
1212 0 : RC = RC )
1213 0 : IF ( RC /= HCO_SUCCESS ) THEN
1214 0 : MSG = 'Cannot convert units for ' // TRIM(Lct%Dct%cName)
1215 0 : CALL HCO_ERROR( MSG , RC )
1216 0 : RETURN
1217 : ENDIF
1218 :
1219 : ! Verbose mode
1220 0 : IF ( UnitFactor /= 1.0_hp ) THEN
1221 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1222 0 : WRITE(MSG,*) 'Data was in units of ', TRIM(thisUnit), &
1223 0 : ' - converted to HEMCO units by applying ', &
1224 0 : 'scale factor ', UnitFactor
1225 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1226 : ENDIF
1227 : ELSE
1228 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1229 0 : WRITE(MSG,*) 'Data was in units of ', TRIM(thisUnit), &
1230 0 : ' - unit conversion factor is ', UnitFactor
1231 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1232 : ENDIF
1233 : ENDIF
1234 :
1235 : ! Check for valid unit combinations, i.e. emissions must be kg/m2/s,
1236 : ! concentrations kg/m3. Eventually multiply by emission time step
1237 : ! or divide by area to obtain those values.
1238 :
1239 : ! Concentration data
1240 0 : IF ( AreaFlag == 3 .AND. TimeFlag == 0 ) THEN
1241 0 : Lct%Dct%Dta%IsConc = .TRUE.
1242 :
1243 : ! If concentration data is per second (kg/m3/s), multiply by emission
1244 : ! time step to get concentration (kg/m3).
1245 0 : ELSEIF ( AreaFlag == 3 .AND. TimeFlag == 1 ) THEN
1246 0 : Lct%Dct%Dta%IsConc = .TRUE.
1247 :
1248 0 : ncArr = ncArr * HcoState%TS_EMIS
1249 : MSG = 'Data converted from kg/m3/s to kg/m3: ' // &
1250 0 : TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
1251 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
1252 :
1253 : ! Unitless data
1254 0 : ELSEIF ( AreaFlag == -1 .AND. TimeFlag == -1 ) THEN
1255 : ! nothing do to
1256 :
1257 : ! Emission data
1258 0 : ELSEIF ( AreaFlag == 2 .AND. TimeFlag == 1 ) THEN
1259 : ! nothing do to
1260 :
1261 : ! Emission data that is not per time (kg/m2): convert to kg/m2/s
1262 0 : ELSEIF ( AreaFlag == 2 .AND. TimeFlag == 0 ) THEN
1263 0 : ncArr = ncArr / HcoState%TS_EMIS
1264 : MSG = 'Data converted from kg/m2 to kg/m2/s: ' // &
1265 0 : TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
1266 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
1267 :
1268 : ! Emission data that is not per area (i.e. kg/s) needs to be converted
1269 : ! to per area manually.
1270 0 : ELSEIF ( AreaFlag == 0 .AND. TimeFlag == 1 ) THEN
1271 :
1272 : ! Get lat edges: those are read from file if possible, otherwise
1273 : ! calculated from the lat midpoints.
1274 : ! ==> Sine of lat is needed. Do conversion right here.
1275 : CALL NC_GET_GRID_EDGES ( ncLun, 2, LatMid, nlat, &
1276 0 : LatEdge, nlatEdge, NCRC )
1277 0 : IF ( NCRC /= 0 ) THEN
1278 0 : MSG = 'Cannot read lat edge of ' // TRIM(srcFile)
1279 0 : CALL HCO_ERROR( MSG, RC )
1280 0 : RETURN
1281 : ENDIF
1282 :
1283 : ! Now normalize data by area calculated from lat edges.
1284 : CALL NORMALIZE_AREA( HcoState, ncArr, nlon, &
1285 0 : LatEdge, srcFile, RC )
1286 0 : IF ( RC /= HCO_SUCCESS ) THEN
1287 0 : CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
1288 0 : RETURN
1289 : ENDIF
1290 :
1291 : ! All other combinations are invalid
1292 : ELSE
1293 : MSG = 'Unit must be unitless, emission or concentration: ' // &
1294 0 : TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
1295 0 : CALL HCO_ERROR( MSG, RC )
1296 0 : RETURN
1297 : ENDIF
1298 : ENDIF ! Unit conversion
1299 :
1300 : !-----------------------------------------------------------------
1301 : ! Get horizontal grid edges
1302 : !-----------------------------------------------------------------
1303 :
1304 : ! Get longitude edges and make sure they are steadily increasing.
1305 : CALL NC_GET_GRID_EDGES ( ncLun, 1, LonMid, nlon, &
1306 0 : LonEdge, nlonEdge, NCRC )
1307 0 : IF ( NCRC /= 0 ) THEN
1308 0 : MSG = 'Cannot read lon edge of ' // TRIM(srcFile)
1309 0 : CALL HCO_ERROR( MSG, RC )
1310 0 : RETURN
1311 : ENDIF
1312 0 : CALL HCO_ValidateLon( HcoState, nlonEdge, LonEdge, RC )
1313 0 : IF ( RC /= HCO_SUCCESS ) THEN
1314 0 : CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
1315 0 : RETURN
1316 : ENDIF
1317 :
1318 : ! Get latitude edges (only if they have not been read yet
1319 : ! for unit conversion)
1320 0 : IF ( .NOT. ASSOCIATED( LatEdge ) ) THEN
1321 : CALL NC_GET_GRID_EDGES ( ncLun, 2, LatMid, nlat, &
1322 0 : LatEdge, nlatEdge, NCRC )
1323 0 : IF ( NCRC /= 0 ) THEN
1324 0 : MSG = 'Cannot read lat edge of ' // TRIM(srcFile)
1325 0 : CALL HCO_ERROR( MSG, RC )
1326 0 : RETURN
1327 : ENDIF
1328 : ENDIF
1329 :
1330 : !-----------------------------------------------------------------
1331 : ! Determine regridding algorithm to be applied: use NCREGRID from
1332 : ! MESSy only if we need to regrid vertical levels. For all other
1333 : ! fields, use the much faster map_a2a.
1334 : ! Perform no vertical regridding if the vertical levels are model
1335 : ! levels. Model levels are assumed to start at the surface, i.e.
1336 : ! the first input level must correspond to the surface level. The
1337 : ! total number of vertical levels must not match the number of
1338 : ! vertical levels on the simulation grid. Data is not extrapolated
1339 : ! beyond the existing levels.
1340 : ! Vertical regridding based on NCREGRID will always map the input
1341 : ! data onto the entire simulation grid (no extrapolation beyond
1342 : ! the vertical input coordinates).
1343 : ! Index-based remapping can currently not be done with the MESSy
1344 : ! routines, i.e. it is not possible to vertically regrid index-
1345 : ! based data.
1346 : !-----------------------------------------------------------------
1347 :
1348 0 : UseMESSy = .FALSE.
1349 0 : IF ( nlev > 1 .AND. .NOT. IsModelLevel ) THEN
1350 0 : UseMESSy = .TRUE.
1351 : ENDIF
1352 :
1353 : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
1354 : ! If in WRF or the CESM environment, the vertical grid is arbitrary.
1355 : ! MESSy regridding ALWAYS has to be used.
1356 0 : IF ( nlev > 1 ) THEN
1357 0 : UseMESSy = .TRUE.
1358 :
1359 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1360 0 : WRITE(MSG,*) ' ==> WRF/CESM: Always forcing MESSy regridding for number of verticals', nlev, IsModelLevel
1361 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1362 : ENDIF
1363 : ENDIF
1364 : #endif
1365 :
1366 0 : IF ( HCO_IsIndexData(Lct%Dct%Dta%OrigUnit) .AND. UseMESSy ) THEN
1367 : MSG = 'Cannot do MESSy regridding for index data: ' // &
1368 0 : TRIM(srcFile)
1369 0 : CALL HCO_ERROR( MSG, RC )
1370 0 : RETURN
1371 : ENDIF
1372 :
1373 : !-----------------------------------------------------------------
1374 : ! Use MESSy regridding
1375 : !-----------------------------------------------------------------
1376 0 : IF ( UseMESSy ) THEN
1377 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1378 0 : WRITE(MSG,*) ' ==> Use MESSy regridding (NCREGRID)'
1379 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1380 : ENDIF
1381 :
1382 : #if !defined( MODEL_CESM ) && !defined( MODEL_WRF )
1383 : ! If we do MESSy regridding, we can only do one time step
1384 : ! at a time at the moment!
1385 : IF ( tidx1 /= tidx2 ) THEN
1386 : MSG = 'Cannot do MESSy regridding for more than one time step; ' &
1387 : // TRIM(srcFile)
1388 : CALL HCO_ERROR( MSG, RC )
1389 : RETURN
1390 : ENDIF
1391 :
1392 : ! Note: This seems to be a soft restriction - removing this
1393 : ! does not conflict with MESSy regridding. Need to check (hplin, 5/30/20)
1394 : ! This has to be used for WRF-GC and CESM so ifdefd out
1395 : #endif
1396 :
1397 : #if defined( MODEL_WRF ) || defined( MODEL_CESM )
1398 : !--------------------------------------------------------------
1399 : ! Eventually get sigma levels
1400 : ! For files that have hardcoded GEOS-Chem "index"-based levels,
1401 : ! translate these levels back into a sigma representation
1402 : ! of the GEOS-Chem levels (sigma = p/ps on INTERFACE)
1403 : !
1404 : ! There are caveats with this. This is essentially a copy of the
1405 : ! hardcoded hPa lists from
1406 : ! http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_vertical_grids
1407 : ! hard-coded by hand, and we only assume that the data is either
1408 : ! 47-levels or 72-levels.
1409 : !
1410 : ! Parse the 72 list using regex like so: ^ ?\d{1,2} then remove the lines
1411 : ! Then you have the 73 edges.
1412 : !
1413 : ! psfc = PEDGE(0) = 1013.250 hPa
1414 : !
1415 : ! Ported from the original WRF-GC implementation (hplin, 5/27/20)
1416 : !--------------------------------------------------------------
1417 0 : IF ( nlev > 1 .AND. IsModelLevel ) THEN
1418 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1419 0 : WRITE(MSG,*) ' ==> WRF/CESM: Writing in fixed sigma coordinates for GEOS-Chem levels', nlon, nlat
1420 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1421 : ENDIF
1422 :
1423 0 : ALLOCATE(SigEdge(nlon, nlat, nlev))
1424 :
1425 0 : DO I = 1, nlon
1426 0 : DO J = 1, nlat
1427 : ! Fill with pre-defined, hard coded sigma levels computed.
1428 0 : SigEdge(I, J, :) = GC_72_EDGE_SIGMA(1:nlev)
1429 : ENDDO
1430 : ENDDO
1431 : ENDIF
1432 : #endif
1433 :
1434 : !--------------------------------------------------------------
1435 : ! Eventually get sigma levels
1436 : ! Vertical regridding is performed on sigma interface levels:
1437 : ! sigma(i,j,l) = p(i,j,l) / ps(i,j)
1438 : ! NC_Get_Sigma_Levels attempts to create the sigma levels from
1439 : ! the content of the netCDF file.
1440 : ! For now, it is assumed that all input data is on vertical
1441 : ! mid-point levels, and the interface values are calculated
1442 : ! by linear interpolation of the mid-point values in a second
1443 : ! step.
1444 : ! For model levels, the sigma levels don't need to be known
1445 : ! as vertical interpolation will be done based on subroutine
1446 : ! ModelLev_Interpolate (within HCO_MESSY_REGRID).
1447 : !--------------------------------------------------------------
1448 0 : IF ( nlev > 1 .AND. .NOT. IsModelLevel ) THEN
1449 :
1450 : ! Get sigma levels
1451 : CALL NC_Get_Sigma_Levels ( fID = ncLun, &
1452 : ncFile = srcFile, &
1453 : levName = LevName, &
1454 : lon1 = 1, &
1455 : lon2 = nlon, &
1456 : lat1 = 1, &
1457 : lat2 = nlat, &
1458 : lev1 = 1, &
1459 : lev2 = nlev, &
1460 : time = tidx1, &
1461 : SigLev = SigLev, &
1462 : Dir = dir, &
1463 0 : RC = NCRC )
1464 0 : IF ( NCRC /= 0 ) THEN
1465 0 : CALL HCO_ERROR( 'Cannot read sigma levels of '//TRIM(srcFile), RC )
1466 0 : RETURN
1467 : ENDIF
1468 :
1469 : ! Interpolate onto edges
1470 0 : CALL SigmaMidToEdges ( HcoState, SigLev, SigEdge, RC )
1471 0 : IF ( RC /= HCO_SUCCESS ) THEN
1472 0 : CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
1473 0 : RETURN
1474 : ENDIF
1475 :
1476 : ! Sigma levels are not needed anymore
1477 0 : IF ( ASSOCIATED(SigLev) ) DEALLOCATE(SigLev)
1478 :
1479 : !-----------------------------------------------------------
1480 : ! Flip vertical axis if positive axis is 'down', i.e. level
1481 : ! index 1 is the top of the atmosphere
1482 : !-----------------------------------------------------------
1483 0 : IF ( dir == -1 ) THEN
1484 0 : SigEdge(:,:,: ) = SigEdge(:,:,nlev+1:1:-1 )
1485 0 : NcArr (:,:,:,:) = NcArr (:,:,nlev :1:-1,:)
1486 : ENDIF
1487 :
1488 : ENDIF ! nlev>1
1489 :
1490 : #if defined( MODEL_WRF ) || defined( MODEL_CESM )
1491 : ! Input data is "never" on model levels because model levels can change! (hplin, 5/29/20)
1492 0 : IsModelLevel = .false.
1493 : #endif
1494 :
1495 : ! Now do the regridding
1496 : CALL HCO_MESSY_REGRID ( HcoState, NcArr, &
1497 : LonEdge, LatEdge, SigEdge, &
1498 0 : Lct, IsModelLevel, RC )
1499 0 : IF ( RC /= HCO_SUCCESS ) THEN
1500 0 : CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
1501 0 : RETURN
1502 : ENDIF
1503 :
1504 : ! Cleanup
1505 0 : IF ( ASSOCIATED(SigEdge) ) DEALLOCATE(SigEdge)
1506 :
1507 : !-----------------------------------------------------------------
1508 : ! Use map_a2a regridding
1509 : !-----------------------------------------------------------------
1510 : ELSE
1511 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1512 0 : WRITE(MSG,*) ' ==> Use map_a2a regridding'
1513 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1514 : ENDIF
1515 :
1516 0 : CALL REGRID_MAPA2A ( HcoState, NcArr, LonEdge, LatEdge, Lct, RC )
1517 0 : IF ( RC /= HCO_SUCCESS ) THEN
1518 0 : CALL HCO_ERROR( 'ERROR 16', RC, THISLOC=LOC )
1519 0 : RETURN
1520 : ENDIF
1521 :
1522 : ENDIF
1523 :
1524 : !-----------------------------------------------------------------
1525 : ! Add to diagnostics (if it exists)
1526 : !-----------------------------------------------------------------
1527 0 : IF ( HcoState%Options%Field2Diagn ) THEN
1528 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 .AND. ASSOCIATED(Lct%Dct%Dta%V3) ) THEN
1529 0 : IF ( ASSOCIATED(Lct%Dct%Dta%V3(1)%Val) ) THEN
1530 : CALL Diagn_Update ( HcoState, cName=TRIM(Lct%Dct%cName), &
1531 0 : Array3D=Lct%Dct%Dta%V3(1)%Val, COL=-1, RC=RC )
1532 0 : IF ( RC /= HCO_SUCCESS ) THEN
1533 0 : CALL HCO_ERROR( 'ERROR 17', RC, THISLOC=LOC )
1534 0 : RETURN
1535 : ENDIF
1536 : ENDIF
1537 0 : ELSEIF ( Lct%Dct%Dta%SpaceDim == 2 .AND. ASSOCIATED(Lct%Dct%Dta%V2) ) THEN
1538 0 : IF ( ASSOCIATED(Lct%Dct%Dta%V2(1)%Val) ) THEN
1539 : CALL Diagn_Update ( HcoState, cName=TRIM(Lct%Dct%cName), &
1540 0 : Array2D=Lct%Dct%Dta%V2(1)%Val, COL=-1, RC=RC )
1541 0 : IF ( RC /= HCO_SUCCESS ) THEN
1542 0 : CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
1543 0 : RETURN
1544 : ENDIF
1545 : ENDIF
1546 : ENDIF
1547 : ENDIF
1548 :
1549 : !-----------------------------------------------------------------
1550 : ! Cleanup and leave
1551 : !-----------------------------------------------------------------
1552 0 : IF ( ASSOCIATED ( ncArr ) ) DEALLOCATE ( ncArr )
1553 0 : IF ( ASSOCIATED ( LonMid ) ) DEALLOCATE ( LonMid )
1554 0 : IF ( ASSOCIATED ( LatMid ) ) DEALLOCATE ( LatMid )
1555 0 : IF ( ASSOCIATED ( LevMid ) ) DEALLOCATE ( LevMid )
1556 0 : IF ( ASSOCIATED ( LonEdge ) ) DEALLOCATE ( LonEdge )
1557 0 : IF ( ASSOCIATED ( LatEdge ) ) DEALLOCATE ( LatEdge )
1558 :
1559 : ! Return w/ success
1560 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
1561 :
1562 0 : END SUBROUTINE HCOIO_Read
1563 : !EOC
1564 : !------------------------------------------------------------------------------
1565 : ! Harmonized Emissions Component (HEMCO) !
1566 : !------------------------------------------------------------------------------
1567 : !BOP
1568 : !
1569 : ! !IROUTINE: HCOIO_CloseAll
1570 : !
1571 : ! !DESCRIPTION: Subroutine HCOIO\_CloseAll makes sure that there is no open
1572 : ! netCDF file left in the stream.
1573 : !\\
1574 : !\\
1575 : ! !INTERFACE:
1576 : !
1577 0 : SUBROUTINE HCOIO_CloseAll( HcoState, RC )
1578 : !
1579 : ! !USES:
1580 : !
1581 : USE HCO_Ncdf_Mod, ONLY : NC_CLOSE
1582 : !
1583 : ! !INPUT PARAMTERS:
1584 : !
1585 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
1586 : !
1587 : ! !INPUT/OUTPUT PARAMETERS:
1588 : !
1589 : INTEGER, INTENT(INOUT) :: RC
1590 : !
1591 : ! !REVISION HISTORY:
1592 : ! 24 Mar 2016 - C. Keller: Initial version
1593 : ! See https://github.com/geoschem/hemco for complete history
1594 : !EOP
1595 : !------------------------------------------------------------------------------
1596 : !BOC
1597 : !
1598 : ! !LOCAL VARIABLES:
1599 : !
1600 : !======================================================================
1601 : ! HCOIO_CloseAll begins here
1602 : !======================================================================
1603 0 : IF ( HcoState%ReadLists%FileLun > 0 ) THEN
1604 0 : CALL NC_CLOSE( HcoState%ReadLists%FileLun )
1605 0 : HcoState%ReadLists%FileLun = -1
1606 : ENDIF
1607 :
1608 : ! Return w/ success
1609 0 : RC = HCO_SUCCESS
1610 :
1611 0 : END SUBROUTINE HCOIO_CloseAll
1612 : !EOC
1613 : END MODULE HCOIO_Read_Mod
1614 : #endif
|