Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hcoio_util_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCOIO\_Util\_Mod contains utility functions
9 : ! for use in data processing including file reading, unit conversions,
10 : ! and regridding.
11 : !\\
12 : !\\
13 : ! !INTERFACE:
14 : !
15 : MODULE HCOIO_Util_Mod
16 : !
17 : ! !USES:
18 : !
19 : USE HCO_Types_Mod
20 : USE HCO_Error_Mod
21 : USE HCO_CharTools_Mod
22 : USE HCO_State_Mod, ONLY : Hco_State
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 : !
27 : ! !PUBLIC MEMBER FUNCTIONS:
28 : !
29 : #if !defined(ESMF_)
30 : PUBLIC :: GET_TIMEIDX
31 : PUBLIC :: Check_AvailYMDhm
32 : PUBLIC :: prefYMDhm_Adjust
33 : PUBLIC :: Set_tIdx2
34 : PUBLIC :: IsClosest
35 : PUBLIC :: GetIndex2Interp
36 : PUBLIC :: GetWeights
37 : PUBLIC :: YMDhm2hrs
38 : PUBLIC :: Normalize_Area
39 : PUBLIC :: SrcFile_Parse
40 : PUBLIC :: SigmaMidToEdges
41 : PUBLIC :: CheckMissVal
42 : PUBLIC :: GetArbDimIndex
43 : #endif
44 : PUBLIC :: HCOIO_ReadOther
45 : PUBLIC :: HCOIO_ReadCountryValues
46 : PUBLIC :: HCOIO_ReadFromConfig
47 : PUBLIC :: GetDataVals
48 : PUBLIC :: GetSliceIdx
49 : PUBLIC :: FillMaskBox
50 : PUBLIC :: ReadMath
51 : !
52 : ! !REVISION HISTORY:
53 : ! 12 Jun 2020 - E. Lundgren - Initial version, created from subset of
54 : ! hcoio_util_mod.F90
55 : ! See https://github.com/geoschem/hemco for complete history
56 : !EOP
57 : !------------------------------------------------------------------------------
58 : !BOC
59 : !
60 : ! !DEFINED PARAMETERS
61 : !
62 : ! Parameter used for difference testing of floating points
63 : REAL(dp), PRIVATE, PARAMETER :: EPSILON = 1.0e-5_dp
64 :
65 : CONTAINS
66 : !EOC
67 : #if !defined( ESMF_ )
68 : !------------------------------------------------------------------------------
69 : ! Harmonized Emissions Component (HEMCO) !
70 : !------------------------------------------------------------------------------
71 : !BOP
72 : !
73 : ! !IROUTINE: Get_TimeIdx
74 : !
75 : ! !DESCRIPTION: Returns the lower and upper time slice index (tidx1
76 : ! and tidx2, respectively) to be read. These values are determined
77 : ! based upon the time slice information extracted from the netCDF file,
78 : ! the time stamp settings set in the config. file, and the current
79 : ! simulation date.
80 : !\\
81 : !\\
82 : ! Return arguments wgt1 and wgt2 denote the weights to be given to
83 : ! the two time slices. This is only of relevance for data that shall
84 : ! be interpolated between two (not necessarily consecutive) time slices.
85 : ! In all other cases, the returned weights are negative and will be
86 : ! ignored.
87 : !\\
88 : !\\
89 : ! Also returns the time slice year and month, as these values may be
90 : ! used for unit conversion.
91 : !\\
92 : !\\
93 : ! !INTERFACE:
94 : !
95 0 : SUBROUTINE GET_TIMEIDX( HcoState, Lct, &
96 : ncLun, tidx1, tidx2, &
97 : wgt1, wgt2, oYMDhm, &
98 : YMDhm, YMDhm1, RC, &
99 : Year )
100 : !
101 : ! !USES:
102 : !
103 : USE HCO_Ncdf_Mod, ONLY : NC_Read_Time_YYYYMMDDhhmm
104 : USE HCO_tIdx_Mod, ONLY : HCO_GetPrefTimeAttr
105 : !
106 : ! !INPUT PARAMETERS:
107 : !
108 : TYPE(HCO_State), POINTER :: HcoState ! HcoState object
109 : TYPE(ListCont), POINTER :: Lct ! List container
110 : INTEGER, INTENT(IN ) :: ncLun ! open ncLun
111 : INTEGER, INTENT(IN ), OPTIONAL :: Year ! year to be used
112 : !
113 : ! !OUTPUT PARAMETERS:
114 : !
115 : INTEGER, INTENT( OUT) :: tidx1 ! lower time idx
116 : INTEGER, INTENT( OUT) :: tidx2 ! upper time idx
117 : REAL(sp), INTENT( OUT) :: wgt1 ! weight to tidx1
118 : REAL(sp), INTENT( OUT) :: wgt2 ! weight to tidx2
119 : REAL(dp), INTENT( OUT) :: oYMDhm ! preferred time slice
120 : REAL(dp), INTENT( OUT) :: YMDhm ! selected time slice
121 : REAL(dp), INTENT( OUT) :: YMDhm1 ! 1st time slice in file
122 : !
123 : ! !INPUT/OUTPUT PARAMETERS:
124 : !
125 : INTEGER, INTENT(INOUT) :: RC
126 : !
127 : ! !REVISION HISTORY:
128 : ! 13 Mar 2013 - C. Keller - Initial version
129 : ! See https://github.com/geoschem/hemco for complete history
130 : !EOP
131 : !------------------------------------------------------------------------------
132 : !BOC
133 : !
134 : ! !LOcAL VARIABLES:
135 : !
136 : CHARACTER(LEN=255) :: MSG, LOC
137 : CHARACTER(LEN=1023) :: MSG_LONG
138 : INTEGER :: tidx1a
139 : INTEGER :: nTime, T, CNT, NCRC
140 : INTEGER :: prefYr, prefMt, prefDy, prefHr, prefMn
141 : INTEGER :: refYear
142 : REAL(dp) :: origYMDhm, prefYMDhm
143 0 : REAL(dp), POINTER :: availYMDhm(:)
144 : LOGICAL :: ExitSearch
145 : LOGICAL :: verb
146 :
147 : !=================================================================
148 : ! GET_TIMEIDX begins here
149 : !=================================================================
150 :
151 : ! Initialize
152 0 : LOC = 'GET_TIMEIDX (HCOIO_UTIL_MOD.F90)'
153 :
154 : ! Officially enter Get_TimeIdx
155 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
156 0 : IF ( RC /= HCO_SUCCESS ) THEN
157 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
158 0 : RETURN
159 : ENDIF
160 0 : verb = HCO_IsVerb( HcoState%Config%Err )
161 :
162 : ! Initialize local variables for safety's sake
163 0 : nTime = 0
164 0 : cnt = 0
165 0 : prefYr = 0
166 0 : prefMt = 0
167 0 : prefDy = 0
168 0 : prefHr = 0
169 0 : prefMn = 0
170 0 : refYear = 0
171 0 : origYMDhm = 0
172 0 : prefYMDhm = 0
173 0 : tidx1 = 0
174 0 : tidx2 = 0
175 0 : tidx1a = 0
176 0 : wgt1 = -1.0_sp
177 0 : wgt2 = -1.0_sp
178 0 : oYMDhm = 0.0_dp
179 0 : YMDhm = 0.0_dp
180 0 : YMDhm1 = 0.0_dp
181 0 : ExitSearch = .FALSE.
182 0 : availYMDhm => NULL()
183 :
184 : ! ----------------------------------------------------------------
185 : ! Extract netCDF time slices (YYYYMMDDhhmm)
186 : ! ----------------------------------------------------------------
187 : CALL NC_READ_TIME_YYYYMMDDhhmm( ncLun, nTime, availYMDhm, &
188 0 : refYear=refYear, RC=NCRC )
189 0 : IF ( NCRC /= 0 ) THEN
190 0 : CALL HCO_ERROR( 'NC_READ_TIME_YYYYMMDDhhmm', RC )
191 0 : RETURN
192 : ENDIF
193 :
194 : ! Return warning if netCDF reference year prior to 1901: it seems
195 : ! like there are some problems with that and the time slices can be
196 : ! off by one day!
197 0 : IF ( (refYear <= 1900) .AND. (nTime > 0) ) THEN
198 : MSG = 'ncdf reference year is prior to 1901 - ' // &
199 0 : 'time stamps may be wrong!'
200 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
201 : ENDIF
202 :
203 : ! verbose mode
204 0 : IF ( verb ) THEN
205 0 : write(MSG,*) 'Number of time slices found: ', nTime
206 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
207 0 : IF ( nTime > 0 ) THEN
208 0 : write(MSG,*) 'Time slice range : ', &
209 0 : availYMDhm(1), availYMDhm(nTime)
210 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
211 : ENDIF
212 : ENDIF
213 :
214 : ! ----------------------------------------------------------------
215 : ! Select time slices to read
216 : ! ----------------------------------------------------------------
217 :
218 : ! ----------------------------------------------------------------
219 : ! Get preferred time stamp to read based upon the specs set in the
220 : ! config. file.
221 : ! This can return value -1 for prefHr, indicating that all
222 : ! corresponding time slices shall be read.
223 : ! This call will return -1 for all date attributes if the
224 : ! simulation date is outside of the data range given in the
225 : ! configuration file.
226 : ! ----------------------------------------------------------------
227 : CALL HCO_GetPrefTimeAttr ( HcoState, Lct, &
228 0 : prefYr, prefMt, prefDy, prefHr, prefMn, RC )
229 0 : IF ( RC /= HCO_SUCCESS ) THEN
230 : MSG = &
231 0 : 'Error encountered in HCO_GetPrefTimeAttr for ' // TRIM(Lct%Dct%cName)
232 0 : CALL HCO_ERROR( MSG, RC )
233 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
234 0 : DEALLOCATE(availYMDhm)
235 : availYMDhm => NULL()
236 : ENDIF
237 0 : RETURN
238 : ENDIF
239 :
240 : ! Eventually force preferred year to passed value
241 0 : IF ( PRESENT(Year) ) prefYr = Year
242 :
243 : ! Check if we are outside of provided range
244 0 : IF ( prefYr < 0 .OR. prefMt < 0 .OR. prefDy < 0 ) THEN
245 :
246 : ! This should only happen for 'range' data
247 0 : IF ( Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_RANGE ) THEN
248 0 : MSG = 'Cannot get preferred datetime for ' // TRIM(Lct%Dct%cName)
249 0 : CALL HCO_ERROR( MSG, RC )
250 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
251 0 : DEALLOCATE(availYMDhm)
252 : availYMDhm => NULL()
253 : ENDIF
254 0 : RETURN
255 : ENDIF
256 :
257 : ! If this part of the code gets executed, the data associated
258 : ! with this container shall not be used at the current date.
259 : ! To do so, set the time indeces to -1 and leave right here.
260 0 : tidx1 = -1
261 0 : tidx2 = -1
262 :
263 : ! Leave w/ success
264 0 : CALL HCO_LEAVE( HcoState%Config%Err, RC )
265 0 : RETURN
266 : ENDIF
267 :
268 : ! origYMDhm is the preferred datetime. Store into shadow variable
269 : ! prefYMDhm. prefYMDhm may be adjusted if origYMDhm is outside of the
270 : ! netCDF datetime range.
271 : ! Now put origYMDhm, prefYMDhm in YYYYMMDDhhmm format (bmy, 4/10/17)
272 : origYMDhm = ( DBLE( prefYr ) * 1.0e8_dp ) + &
273 : ( DBLE( prefMt ) * 1.0e6_dp ) + &
274 : ( DBLE( prefDy ) * 1.0e4_dp ) + &
275 : ( DBLE( MAX( prefHr, 0 ) ) * 1.0e2_dp ) + &
276 0 : ( DBLE( MAX( prefMn, 0 ) ) )
277 0 : prefYMDhm = origYMDhm
278 :
279 : ! verbose mode
280 0 : IF ( verb ) THEN
281 0 : write(MSG,'(A30,f14.0)') 'preferred datetime: ', prefYMDhm
282 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
283 : ENDIF
284 :
285 : ! ================================================================
286 : ! Case 1: Only one time slice available.
287 : ! ================================================================
288 0 : IF ( nTime == 1 ) THEN
289 0 : tidx1 = 1
290 0 : tidx2 = 1
291 :
292 : ! ================================================================
293 : ! Case 2: More than one time slice available. Determine lower
294 : ! and upper time slice index from file & HEMCO settings.
295 : ! ================================================================
296 0 : ELSEIF ( nTime > 1 ) THEN
297 :
298 : ! Init
299 0 : tidx1 = -1
300 0 : tidx2 = -1
301 :
302 : ! -------------------------------------------------------------
303 : ! Check if preferred datetime prefYMDhm is within the range
304 : ! available time slices, e.g. it falls within the interval
305 : ! of availYMDhm. In this case, set tidx1 to the index of the
306 : ! closest time slice that is not in the future.
307 : ! -------------------------------------------------------------
308 0 : CALL Check_AvailYMDhm ( Lct, nTime, availYMDhm, prefYMDhm, tidx1a )
309 :
310 : ! -------------------------------------------------------------
311 : ! Check if we need to continue search. Even if the call above
312 : ! returned a time slice, it may be possible to continue looking
313 : ! for a better suited time stamp. This is only the case if
314 : ! there are discontinuities in the time stamps, e.g. if a file
315 : ! contains monthly data for 2005 and 2020. In that case, the
316 : ! call above would return the index for Dec 2005 for any
317 : ! simulation date between 2005 and 2010 (e.g. July 2010),
318 : ! whereas it makes more sense to use July 2005 (and eventually
319 : ! interpolate between the July 2005 and July 2020 data).
320 : ! The IsClosest command checks if there are any netCDF time
321 : ! stamps (prior to the selected one) that are closer to each
322 : ! other than the difference between the preferred time stamp
323 : ! prefYMDhm and the currently selected time stamp
324 : ! availYMDhm(tidx1a). In that case, it continues the search by
325 : ! updating prefYMDhm so that it falls within the range of the
326 : ! 'high-frequency' interval.
327 : ! -------------------------------------------------------------
328 0 : ExitSearch = .FALSE.
329 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) THEN
330 : ExitSearch = .TRUE.
331 0 : ELSE IF ( tidx1a > 0 ) THEN
332 0 : ExitSearch = IsClosest( prefYMDhm, availYMDhm, nTime, tidx1a )
333 : ENDIF
334 :
335 : ! When using the interpolation flag, use the first or last timestep
336 : ! when outside of the available date range
337 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .and. tidx1a < 0 ) THEN
338 0 : IF ( prefYMDhm < availYMDhm(1) ) THEN
339 0 : tidx1a = 1
340 0 : ELSE IF ( prefYMDhm > availYMDhm(nTime) ) THEN
341 0 : tidx1a = nTime
342 : ENDIF
343 : ENDIF
344 :
345 : ! Do not continue search if data is to be interpolated and is
346 : ! not discontinuous (mps, 10/23/19)
347 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .and. &
348 : .not. Lct%Dct%Dta%Discontinuous ) THEN
349 : ExitSearch = .TRUE.
350 : ENDIF
351 :
352 : ! Write to tidx1 if this is the best match.
353 0 : IF ( ExitSearch ) THEN
354 0 : tidx1 = tidx1a
355 :
356 : ! -------------------------------------------------------------
357 : ! If search shall be continued, adjust preferred year, then
358 : ! month, then day to the closest available year (month, day)
359 : ! in the time slices, and check if this is a better match.
360 : ! -------------------------------------------------------------
361 : ELSE
362 :
363 : ! Adjust year, month, and day (in this order).
364 0 : CNT = 0
365 : DO
366 0 : CNT = CNT + 1
367 0 : IF ( ExitSearch .OR. CNT > 3 ) EXIT
368 :
369 : ! Adjust prefYMDhm at the given level (1=Y, 2=M, 3=D)
370 0 : CALL prefYMDhm_Adjust ( nTime, availYMDhm, prefYMDhm, CNT, tidx1a )
371 :
372 : ! verbose mode
373 0 : IF ( verb ) THEN
374 0 : write(MSG,'(A30,f14.0)') 'adjusted preferred datetime: ', &
375 0 : prefYMDhm
376 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
377 : ENDIF
378 :
379 : ! check for time stamp with updated date/time
380 0 : CALL Check_AvailYMDhm ( Lct, nTime, availYMDhm, prefYMDhm, tidx1a )
381 :
382 : ! Can we leave now?
383 0 : ExitSearch = IsClosest( prefYMDhm, availYMDhm, nTime, tidx1a )
384 0 : IF ( ExitSearch ) tidx1 = tidx1a
385 :
386 : ENDDO
387 : ENDIF
388 :
389 : ! -------------------------------------------------------------
390 : ! If tidx1 still isn't defined, i.e. prefYMDhm is still
391 : ! outside the range of availYMDhm, set tidx1 to the closest
392 : ! available date. This must be 1 or nTime!
393 : ! -------------------------------------------------------------
394 0 : IF ( .NOT. ExitSearch ) THEN
395 0 : IF ( prefYMDhm < availYMDhm(1) ) THEN
396 0 : tidx1 = 1
397 : ELSE
398 0 : tidx1 = nTime
399 : ENDIF
400 : ENDIF
401 :
402 : ! -------------------------------------------------------------
403 : ! If we are dealing with 3-hourly or hourly data, select all timesteps
404 : ! -------------------------------------------------------------
405 :
406 : ! Hour flag is -1: wildcard
407 0 : IF ( Lct%Dct%Dta%ncHrs(1) == -1 .AND. nTime == 8 ) THEN
408 0 : tidx1 = 1
409 0 : tidx2 = nTime
410 :
411 : ! verbose mode
412 0 : IF ( verb ) THEN
413 0 : WRITE(MSG,*) 'Data is 3-hourly. Entire day will be read.'
414 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
415 : ENDIF
416 : ENDIF
417 0 : IF ( Lct%Dct%Dta%ncHrs(1) == -1 .AND. nTime == 24 ) THEN
418 0 : tidx1 = 1
419 0 : tidx2 = nTime
420 :
421 : ! verbose mode
422 0 : IF ( verb ) THEN
423 0 : WRITE(MSG,*) 'Data is hourly. Entire day will be read.'
424 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
425 : ENDIF
426 : ENDIF
427 :
428 : ! -------------------------------------------------------------
429 : ! If we are dealing with weekday data, pick the slice to be
430 : ! used based on the current day of week.
431 : ! The ncDys flag has been set in subroutine HCO_ExtractTime
432 : ! (hco_tidx_mod.F90) based upon the time attributes set in the
433 : ! configuration file. It can have the following values:
434 : ! >0 : specific days are given.
435 : ! -1 : wildcard (autodetect)
436 : ! -10 : WD (weekday).
437 : ! -999: determine from current simulation day.
438 : ! For specific days or if determined from the current datetime
439 : ! (flags >0 or -999), the weekday is not taken into account.
440 : ! If auto-detection is enabled, days are treated as weekday if
441 : ! (and only if) there are exactly 7 time slices. Otherwise, they
442 : ! are interpreted as 'regular' day data.
443 : ! If flag is set to -10, e.g. time attribute is 'WD', the current
444 : ! time index is assumed to hold Sunday data, with the following
445 : ! six slices being Mon, Tue, ..., Sat. For weekdaily data, all
446 : ! seven time slices will be read into memory so that at any given
447 : ! time, the local weekday can be taken (weekdaily data is always
448 : ! assumed to be in local time).
449 : ! -------------------------------------------------------------
450 :
451 : ! Day flag is -1: wildcard
452 0 : IF ( Lct%Dct%Dta%ncDys(1) == -1 .AND. nTime == 7 ) THEN
453 0 : tidx1 = 1
454 0 : tidx2 = nTime
455 :
456 : ! Make sure data is treated in local time
457 0 : Lct%Dct%Dta%IsLocTime = .TRUE.
458 :
459 : ! Day flag is -10: WD
460 0 : ELSEIF ( Lct%Dct%Dta%ncDys(1) == -10 ) THEN
461 :
462 : ! There must be at least 7 time slices
463 0 : IF ( nTime < 7 ) THEN
464 : MSG = 'Data must have exactly 7 time slices '// &
465 0 : 'if you set day attribute to WD: '//TRIM(Lct%Dct%cName)
466 0 : CALL HCO_ERROR( MSG, RC )
467 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
468 0 : DEALLOCATE(availYMDhm)
469 : availYMDhm => NULL()
470 : ENDIF
471 0 : RETURN
472 : ENDIF
473 :
474 : ! If there are exactly seven time slices, interpret them as
475 : ! the seven weekdays.
476 0 : IF ( nTime == 7 ) THEN
477 0 : tidx1 = 1
478 0 : tidx2 = 7
479 :
480 : ! If there are more than 7 time slices, interpret the current
481 : ! selected index as sunday of the current time frame (e.g. sunday
482 : ! data of current month), and select the time slice index
483 : ! accordingly. This requires that there are at least 6 more time
484 : ! slices following the current one.
485 : ELSE
486 0 : IF ( tidx1 < 0 ) THEN
487 0 : WRITE(MSG,*) 'Cannot get weekday slices for: ', &
488 0 : TRIM(Lct%Dct%cName), '. Cannot find first time slice.'
489 0 : CALL HCO_ERROR( MSG, RC )
490 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
491 0 : DEALLOCATE(availYMDhm)
492 : availYMDhm => NULL()
493 : ENDIF
494 0 : RETURN
495 : ENDIF
496 :
497 0 : IF ( (tidx1+6) > nTime ) THEN
498 0 : WRITE(MSG,*) 'Cannot get weekday for: ',TRIM(Lct%Dct%cName), &
499 0 : '. There are less than 6 additional time slices after ', &
500 0 : 'selected start date ', availYMDhm(tidx1)
501 0 : CALL HCO_ERROR( MSG, RC )
502 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
503 0 : DEALLOCATE(availYMDhm)
504 : availYMDhm => NULL()
505 : ENDIF
506 0 : RETURN
507 : ENDIF
508 0 : tidx2 = tidx1 + 6
509 : ENDIF
510 :
511 : ! Make sure data is treated in local time
512 0 : Lct%Dct%Dta%IsLocTime = .TRUE.
513 :
514 : ENDIF
515 :
516 : ! -------------------------------------------------------------
517 : ! Now need to set upper time slice index tidx2. This index
518 : ! is only different from tidx1 if:
519 : ! (1) We interpolate between two time slices, i.e. TimeCycle
520 : ! attribute is set to 'I'. In this case, we simply pick
521 : ! the next higher time slice index and calculate the
522 : ! weights for time1 and time2 based on the current time.
523 : ! (2) Multiple hourly slices are read (--> prefHr = -1 or -10,
524 : ! e.g. hour attribute in config. file was set to wildcard
525 : ! character or data is in local hours). In this case,
526 : ! check if there are multiple time slices for the selected
527 : ! date (y/m/d).
528 : ! tidx2 has already been set to proper value above if it's
529 : ! weekday data.
530 : ! -------------------------------------------------------------
531 0 : IF ( tidx2 < 0 ) THEN
532 :
533 : ! Interpolate between dates
534 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER ) THEN
535 :
536 : CALL GetIndex2Interp( HcoState, Lct, nTime, &
537 : availYMDhm, prefYMDhm, origYMDhm, &
538 : tidx1, tidx2, wgt1, &
539 0 : wgt2, RC )
540 0 : IF ( RC /= HCO_SUCCESS ) THEN
541 : MSG = 'Error encountered in GetIndex2Interp for: ' // &
542 0 : TRIM(Lct%Dct%Cname)
543 0 : CALL HCO_ERROR( MSG, RC )
544 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
545 0 : DEALLOCATE(availYMDhm)
546 : availYMDhm => NULL()
547 : ENDIF
548 0 : RETURN
549 : ENDIF
550 :
551 : ! Check for multiple hourly data
552 0 : ELSEIF ( tidx1 > 0 .AND. prefHr < 0 ) THEN
553 0 : CALL SET_TIDX2 ( nTime, availYMDhm, tidx1, tidx2 )
554 :
555 : ! Denote as local time if necessary
556 0 : IF ( Lct%Dct%Dta%ncHrs(1) == -10 ) THEN
557 0 : Lct%Dct%Dta%IsLocTime = .TRUE.
558 : ENDIF
559 : ELSE
560 0 : tidx2 = tidx1
561 : ENDIF
562 : ENDIF
563 :
564 : ! ================================================================
565 : ! Case 3: No time slice available. Set both indeces to zero. Data
566 : ! with no time stamp must have CycleFlag 'Cycling'.
567 : ! ================================================================
568 : ELSE
569 0 : IF ( Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_CYCLE ) THEN
570 : MSG = 'Field has no time/date variable - cycle flag must' // &
571 : 'be set to `C` in the HEMCO configuration file:' // &
572 0 : TRIM(Lct%Dct%cName)
573 0 : CALL HCO_ERROR( MSG, RC )
574 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
575 0 : DEALLOCATE(availYMDhm)
576 : availYMDhm => NULL()
577 : ENDIF
578 0 : RETURN
579 : ENDIF
580 :
581 0 : tidx1 = 0
582 0 : tidx2 = 0
583 : ENDIF
584 :
585 : !-----------------------------------------------------------------
586 : ! Sanity check: if CycleFlag is set to 'Exact', the file time stamp
587 : ! must exactly match the current time.
588 : !-----------------------------------------------------------------
589 0 : IF ( (Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT) .AND. (tidx1 > 0) ) THEN
590 0 : IF ( availYMDhm(tidx1) /= prefYMDhm ) THEN
591 0 : tidx1 = -1
592 0 : tidx2 = -1
593 : ENDIF
594 : ENDIF
595 :
596 : !-----------------------------------------------------------------
597 : ! If multiple time slices are read, extract time interval between
598 : ! time slices in memory (in hours). This is to make sure that the
599 : ! cycling between the slices will be done at the correct rate
600 : ! (e.g. every hour, every 3 hours, ...).
601 : !-----------------------------------------------------------------
602 0 : IF ( (tidx2>tidx1) .AND. (Lct%Dct%Dta%CycleFlag/=HCO_CFLAG_INTER) ) THEN
603 0 : Lct%Dct%Dta%DeltaT = YMDhm2hrs( availYMDhm(tidx1+1) - availYMDhm(tidx1) )
604 : ELSE
605 0 : Lct%Dct%Dta%DeltaT = 0
606 : ENDIF
607 :
608 : ! verbose mode
609 0 : IF ( verb ) THEN
610 0 : WRITE(MSG,'(A30,I14)') 'selected tidx1: ', tidx1
611 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
612 0 : IF ( tidx1 > 0 ) THEN
613 0 : WRITE(MSG,'(A30,f14.0)') 'corresponding datetime 1: ', &
614 0 : availYMDhm(tidx1)
615 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
616 0 : IF ( wgt1 >= 0.0_sp ) THEN
617 0 : WRITE(MSG,*) 'weight1: ', wgt1
618 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
619 : ENDIF
620 : ENDIF
621 :
622 0 : IF ( (tidx2 /= tidx1) ) THEN
623 0 : WRITE(MSG,'(A30,I14)') 'selected tidx2: ', tidx2
624 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
625 0 : WRITE(MSG,'(A30,f14.0)') 'corresponding datetime 2: ', &
626 0 : availYMDhm(tidx2)
627 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
628 0 : IF ( wgt1 >= 0.0_sp ) THEN
629 0 : WRITE(MSG,*) 'weight2: ', wgt2
630 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
631 : ENDIF
632 : ENDIF
633 :
634 0 : WRITE(MSG,'(A30,I14)') 'assigned delta t [h]: ', Lct%Dct%Dta%DeltaT
635 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
636 0 : WRITE(MSG,*) 'local time? ', Lct%Dct%Dta%IsLocTime
637 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
638 : ENDIF
639 :
640 : ! ----------------------------------------------------------------
641 : ! TODO: set time brackets
642 : ! --> In future, we may want to set time brackets denoting the
643 : ! previous and next time slice available in the netCDF file. This
644 : ! may become useful for temporal interpolations and more efficient
645 : ! data update calls (only update if new time slice is available).
646 : ! ----------------------------------------------------------------
647 :
648 : !-----------------------------------------------------------------
649 : ! Prepare output, cleanup and leave
650 : !-----------------------------------------------------------------
651 :
652 : ! ncYr and ncMt are the year and month fo the time slice to be
653 : ! used. These values may be required to convert units to 'per
654 : ! seconds'.
655 0 : IF ( tidx1 > 0 ) THEN
656 0 : YMDhm = availYMDhm(tidx1)
657 0 : YMDhm1 = availYMDhm(1)
658 0 : oYMDhm = origYMDhm
659 : ENDIF
660 :
661 : ! Deallocate and nullify the pointer
662 0 : IF ( ASSOCIATED(availYMDhm) ) THEN
663 0 : DEALLOCATE(availYMDhm)
664 : availYMDhm => NULL()
665 : ENDIF
666 :
667 : ! Return w/ success
668 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
669 :
670 0 : END SUBROUTINE GET_TIMEIDX
671 : !EOC
672 : !------------------------------------------------------------------------------
673 : ! Harmonized Emissions Component (HEMCO) !
674 : !------------------------------------------------------------------------------
675 : !BOP
676 : !
677 : ! !IROUTINE: Check_AvailYMDhm
678 : !
679 : ! !DESCRIPTION: Checks if prefYMDhm is within the range of availYMDhm
680 : ! and returns the location of the closest vector element that is in
681 : ! the past (--> tidx1). tidx1 is set to -1 otherwise.
682 : !\\
683 : !\\
684 : ! !INTERFACE:
685 : !
686 0 : SUBROUTINE Check_AvailYMDhm( Lct, N, availYMDhm, prefYMDhm, tidx1 )
687 : !
688 : ! !INPUT PARAMETERS:
689 : !
690 : TYPE(ListCont), POINTER :: Lct
691 : INTEGER, INTENT(IN) :: N
692 : REAL(dp), INTENT(IN) :: availYMDhm(N)
693 : REAL(dp), INTENT(IN) :: prefYMDhm
694 : !
695 : ! !OUTPUT PARAMETERS:
696 : !
697 : INTEGER, INTENT(OUT) :: tidx1
698 : !
699 : ! !REVISION HISTORY:
700 : ! 13 Mar 2013 - C. Keller - Initial version
701 : ! See https://github.com/geoschem/hemco for complete history
702 : !EOP
703 : !------------------------------------------------------------------------------
704 : !BOC
705 : !
706 : ! !LOCAL VARIABLES:
707 : !
708 : INTEGER :: I, nTime
709 :
710 : !=================================================================
711 : ! Check_availYMDhm begins here
712 : !=================================================================
713 :
714 : ! Init
715 0 : tidx1 = -1
716 :
717 : ! Return if preferred datetime not within the vector range
718 0 : IF ( prefYMDhm < availYMDhm(1) .OR. prefYMDhm > availYMDhm(N) ) RETURN
719 :
720 : ! To avoid out-of-bounds error in the loop below:
721 : ! (1) For interpolated data, the upper loop limit should be N;
722 : ! (2) Otherwise, the upper loop limit should be N-1.
723 : ! (bmy, 4/28/21)
724 0 : nTime = N - 1
725 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER ) nTime = N
726 :
727 : ! Get closest index that is not in the future
728 0 : DO I = 1, nTime
729 :
730 : ! NOTE: Epsilon test is more robust than an equality test
731 : ! for double-precision variables (bmy, 4/11/17)
732 0 : IF ( ABS( availYMDhm(I) - prefYMDhm ) < EPSILON ) THEN
733 0 : tidx1 = I
734 0 : EXIT
735 : ENDIF
736 :
737 : ! Check if next time slice is in the future, in which case the
738 : ! current slice is selected. Don't do this for a CycleFlag of
739 : ! 3 (==> exact match).
740 0 : IF ( (availYMDhm(I+1) > prefYMDhm ) .AND. &
741 0 : (Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_EXACT) ) THEN
742 0 : tidx1 = I
743 0 : EXIT
744 : ENDIF
745 : ENDDO
746 :
747 : END SUBROUTINE Check_AvailYMDhm
748 : !EOC
749 : !------------------------------------------------------------------------------
750 : ! Harmonized Emissions Component (HEMCO) !
751 : !------------------------------------------------------------------------------
752 : !BOP
753 : !
754 : ! !IROUTINE: prefYMDhm_Adjust
755 : !
756 : ! !DESCRIPTION: Adjusts prefYMDhm to the closest available time attribute. Can
757 : ! be adjusted for year (level=1), month (level=2), or day (level=3).
758 : !\\
759 : !\\
760 : ! !INTERFACE:
761 : !
762 0 : SUBROUTINE prefYMDhm_Adjust( N, availYMDhm, prefYMDhm, level, tidx1 )
763 : !
764 : ! !INPUT PARAMETERS:
765 : !
766 : INTEGER , INTENT(IN) :: N
767 : REAL(dp) , INTENT(IN) :: availYMDhm(N)
768 : INTEGER , INTENT(IN) :: level
769 : INTEGER , INTENT(IN) :: tidx1
770 : !
771 : ! !INPUT/OUTPUT PARAMETERS:
772 : !
773 : REAL(dp) , INTENT(INOUT) :: prefYMDhm
774 : !
775 : ! !REVISION HISTORY:
776 : ! 13 Mar 2013 - C. Keller - Initial version
777 : ! See https://github.com/geoschem/hemco for complete history
778 : !EOP
779 : !------------------------------------------------------------------------------
780 : !BOC
781 : !
782 : ! !LOCAL VARIABLES:
783 : !
784 : ! Scalars
785 : INTEGER :: I, IMIN, IMAX
786 : REAL(dp) :: origYr, origMt, origDy, origHr, origMi
787 : REAL(dp) :: refAttr, tmpAttr, newAttr
788 : REAL(dp) :: iDiff, minDiff
789 : REAL(dp) :: modVal
790 : REAL(dp) :: div
791 :
792 : !=================================================================
793 : ! prefYMDhm_Adjust begins here!
794 : !=================================================================
795 :
796 : ! Get original Yr, Mt, Day, Hr, Mi
797 : ! Time values are now in YYYYMMDDhhmm format (bmy, 4/11/17)
798 0 : origYr = FLOOR( MOD( prefYMDhm, 1.0e12_dp ) / 1.0e8_dp )
799 0 : origMt = FLOOR( MOD( prefYMDhm, 1.0e8_dp ) / 1.0e6_dp )
800 0 : origDy = FLOOR( MOD( prefYMDhm, 1.0e6_dp ) / 1.0e4_dp )
801 0 : origHr = FLOOR( MOD( prefYMDhm, 1.0e4_dp ) / 1.0e2_dp )
802 0 : origMi = FLOOR( MOD( prefYMDhm, 1.0e2_dp ) )
803 :
804 : ! Extract new attribute from availYMDhm and insert into prefYMDhm. Pick
805 : ! closest available value.
806 0 : SELECT CASE ( level )
807 : ! --- Year
808 : CASE ( 1 )
809 : modVal = 1.0e12_dp
810 : div = 1.0e8_dp
811 0 : refAttr = origYr
812 :
813 : ! --- Month
814 : CASE ( 2 )
815 0 : modVal = 1.0e8_dp
816 0 : div = 1.0e6_dp
817 0 : refAttr = origMt
818 :
819 : ! --- Day
820 : CASE ( 3 )
821 0 : modVal = 1.0e6_dp
822 0 : div = 1.0e4_dp
823 0 : refAttr = origMt
824 :
825 : ! --- Hour
826 : CASE ( 4 )
827 0 : modval = 1.0e4_dp
828 0 : div = 1.0e2_dp
829 0 : refAttr = origHr
830 :
831 : ! --- Minute
832 : CASE ( 5 )
833 0 : modVal = 1.0e2_dp
834 0 : div = 1.0_dp
835 0 : refAttr = origMi
836 :
837 : CASE DEFAULT
838 0 : RETURN
839 : END SELECT
840 :
841 : ! Maximum loop number:
842 : ! If tidx1 is already set, only search values in the past.
843 0 : IF ( tidx1 > 0 ) THEN
844 : IMIN = 1
845 : IMAX = tidx1
846 :
847 : ! If tidx1 is not yet set, prefYMDhm must be outside the range of
848 : ! availYMDhm. Pick only the closest available time stamp.
849 : ELSE
850 0 : IF ( prefYMDhm > availYMDhm(1) ) THEN
851 0 : IMIN = N
852 0 : IMAX = N
853 : ELSE
854 : IMIN = 1
855 : IMAX = 1
856 : ENDIF
857 : ENDIF
858 :
859 : ! Select current minimum value
860 0 : minDiff = 10000000000000000.0_dp
861 0 : newAttr = -1d0
862 0 : DO I = IMIN, IMAX
863 0 : tmpAttr = FLOOR( MOD(availYMDhm(I),modVal) / div )
864 0 : iDiff = ABS( tmpAttr - refAttr )
865 0 : IF ( iDiff < minDiff ) THEN
866 0 : newAttr = tmpAttr
867 0 : minDiff = iDiff
868 : ENDIF
869 : ENDDO
870 :
871 : ! Just reuse current value if no better value could be found
872 0 : IF ( newAttr < 0 ) THEN
873 0 : newAttr = refAttr
874 : ENDIF
875 :
876 : ! Update variable
877 : ! --- Year
878 0 : IF ( level == 1 ) THEN
879 : prefYMDhm = ( newAttr * 1.0e8_dp ) + &
880 : ( origMt * 1.0e6_dp ) + &
881 : ( origDy * 1.0e4_dp ) + &
882 : ( origHr * 1.0e2_dp ) + &
883 0 : ( origMi )
884 :
885 : ! --- Month
886 0 : ELSEIF ( level == 2 ) THEN
887 : prefYMDhm = ( origYr * 1.0e8_dp ) + &
888 : ( newAttr * 1.0e6_dp ) + &
889 : ( origDy * 1.0e4_dp ) + &
890 : ( origHr * 1.0e2_dp ) + &
891 0 : ( origMi )
892 :
893 : ! --- Day
894 0 : ELSEIF ( level == 3 ) THEN
895 : prefYMDhm = ( origYr * 1.0e8_dp ) + &
896 : ( origMt * 1.0e6_dp ) + &
897 : ( newAttr * 1.0e4_dp ) + &
898 : ( origHr * 1.0e2_dp ) + &
899 0 : ( origMi )
900 :
901 : ! --- Hour
902 0 : ELSEIF ( level == 4 ) THEN
903 : prefYMDhm = ( origYr * 1.0e8_dp ) + &
904 : ( origMt * 1.0e6_dp ) + &
905 : ( origDy * 1.0e4_dp ) + &
906 : ( newAttr * 1.0e2_dp ) + &
907 0 : ( origMi )
908 : ! --- Minute
909 0 : ELSEIF ( level == 5 ) THEN
910 : prefYMDhm = ( origYr * 1.0e8_dp ) + &
911 : ( origMt * 1.0e6_dp ) + &
912 : ( origDy * 1.0e4_dp ) + &
913 : ( origHr * 1.0e2_dp ) + &
914 0 : ( newAttr )
915 :
916 : ENDIF
917 :
918 : END SUBROUTINE prefYMDhm_Adjust
919 : !EOC
920 : !------------------------------------------------------------------------------
921 : ! Harmonized Emissions Component (HEMCO) !
922 : !------------------------------------------------------------------------------
923 : !BOP
924 : !
925 : ! !IROUTINE: Set_tIdx2
926 : !
927 : ! !DESCRIPTION: sets the upper time slice index by selecting the range
928 : ! of all elements in availYMDhm with the same date (year,month,day) as
929 : ! availYMDh(tidx1).
930 : !\\
931 : !\\
932 : ! !INTERFACE:
933 : !
934 0 : SUBROUTINE Set_tIdx2( N, availYMDhm, tidx1, tidx2 )
935 : !
936 : ! !INPUT PARAMETERS:
937 : !
938 : INTEGER, INTENT(IN) :: N ! Number of times
939 : REAL(dp), INTENT(IN) :: availYMDhm(N) ! Time stamp vector
940 : INTEGER, INTENT(IN) :: tidx1 ! Lower time slice index
941 : !
942 : ! !INPUT/OUTPUT PARAMETERS:
943 : !
944 : INTEGER, INTENT(OUT) :: tidx2 ! Upper time slice index
945 : !
946 : ! !REVISION HISTORY:
947 : ! 13 Mar 2013 - C. Keller - Initial version
948 : ! See https://github.com/geoschem/hemco for complete history
949 : !EOP
950 : !------------------------------------------------------------------------------
951 : !BOC
952 : !
953 : ! !LOCAL VARIABLES:
954 : !
955 : INTEGER :: YMD, I, IYMD
956 :
957 : !=================================================================
958 : ! SET_TIDX2 begins here!
959 : !=================================================================
960 :
961 : ! Init
962 0 : tidx2 = tidx1
963 :
964 : ! Sanity check
965 0 : IF ( tidx1 == N ) RETURN
966 :
967 : ! Get wanted YMD
968 0 : YMD = floor(availYMDhm(tidx1) / 1.0e4_dp)
969 :
970 : ! See how many more tile slices with the same YMD exist from index
971 : ! tidx1 onwards.
972 0 : DO I = tidx1, N
973 0 : iYMD = floor(availYMDhm(I) / 1.0e4_dp)
974 0 : IF ( iYMD == YMD ) THEN
975 0 : tidx2 = I
976 0 : ELSEIF ( iYMD > YMD ) THEN
977 : EXIT
978 : ENDIF
979 : ENDDO
980 :
981 : END SUBROUTINE Set_tIdx2
982 : !EOC
983 : !------------------------------------------------------------------------------
984 : ! Harmonized Emissions Component (HEMCO) !
985 : !------------------------------------------------------------------------------
986 : !BOP
987 : !
988 : ! !IROUTINE: IsClosest
989 : !
990 : ! !DESCRIPTION: function IsClosest returns true if the selected time index
991 : ! is the 'closest' one. It is defined as being closest if:
992 : ! (a) the currently selected index exactly matches the preferred one.
993 : ! (b) the time gap between the preferred time stamp and the currently selected
994 : ! index is at least as small as any other gap of consecutive prior time stamps.
995 : !\\
996 : !\\
997 : ! !INTERFACE:
998 : !
999 0 : FUNCTION IsClosest ( prefYMDhm, availYMDhm, nTime, ctidx1 ) RESULT ( Closest )
1000 : !
1001 : ! !INPUT PARAMETERS:
1002 : !
1003 : INTEGER, INTENT(IN) :: nTime
1004 : REAL(dp), INTENT(IN) :: prefYMDhm
1005 : REAL(dp), INTENT(IN) :: availYMDhm(nTime)
1006 : INTEGER, INTENT(IN) :: ctidx1
1007 : !
1008 : ! !OUTPUT PARAMETERS:
1009 : !
1010 : LOGICAL :: Closest
1011 : !
1012 : ! !REVISION HISTORY:
1013 : ! 03 Mar 2015 - C. Keller - Initial version
1014 : ! See https://github.com/geoschem/hemco for complete history
1015 : !EOP
1016 : !------------------------------------------------------------------------------
1017 : !BOC
1018 : !
1019 : ! !LOCAL VARIABLES:
1020 : !
1021 : INTEGER :: N
1022 : INTEGER :: diff, idiff
1023 :
1024 : !=================================================================
1025 : ! IsClosest begins here!
1026 : !=================================================================
1027 :
1028 : ! Init
1029 0 : Closest = .TRUE.
1030 :
1031 : ! It's not closest if index is not defined
1032 0 : IF ( ctidx1 <= 0 ) THEN
1033 0 : Closest = .FALSE.
1034 : RETURN
1035 : ENDIF
1036 :
1037 : ! It's closest if it is the first index
1038 0 : IF ( ctidx1 == 1 ) RETURN
1039 :
1040 : ! It's closest if it matches date exactly
1041 : ! NOTE: Epsilon test is more robust than an equality test
1042 : ! for double-precision variables (bmy, 4/11/17)
1043 0 : IF ( ABS( availYMDhm(ctidx1) - prefYMDhm ) < EPSILON ) RETURN
1044 :
1045 : ! It's closest if current select one is in the future
1046 0 : IF ( availYMDhm(ctidx1) > prefYMDhm ) RETURN
1047 :
1048 : ! Check if any of the time stamps in the past have closer intervals
1049 : ! than the current select time stamp to it's previous one
1050 0 : diff = prefYMDhm - availYMDhm(ctidx1)
1051 0 : DO N = 2, ctidx1
1052 0 : idiff = availYMDhm(N) - availYMDhm(N-1)
1053 0 : IF ( idiff < diff ) THEN
1054 0 : Closest = .FALSE.
1055 : RETURN
1056 : ENDIF
1057 : ENDDO
1058 :
1059 : END FUNCTION IsClosest
1060 : !EOC
1061 : !------------------------------------------------------------------------------
1062 : ! Harmonized Emissions Component (HEMCO) !
1063 : !------------------------------------------------------------------------------
1064 : !BOP
1065 : !
1066 : ! !IROUTINE: GetIndex2Interp
1067 : !
1068 : ! !DESCRIPTION: GetIndex2Interp
1069 : !\\
1070 : !\\
1071 : ! !INTERFACE:
1072 : !
1073 0 : SUBROUTINE GetIndex2Interp ( HcoState, Lct, &
1074 0 : nTime, availYMDhm, &
1075 : prefYMDhm, origYMDhm, tidx1, &
1076 : tidx2, wgt1, wgt2, RC )
1077 : !
1078 : ! !INPUT PARAMETERS:
1079 : !
1080 : TYPE(HCO_State), POINTER :: HcoState
1081 : TYPE(ListCont), POINTER :: Lct
1082 : INTEGER, INTENT(IN) :: nTime
1083 : REAL(dp), INTENT(IN) :: availYMDhm(nTime)
1084 : REAL(dp), INTENT(IN) :: prefYMDhm
1085 : REAL(dp), INTENT(IN) :: origYMDhm
1086 : INTEGER, INTENT(IN) :: tidx1
1087 : !
1088 : ! !OUTPUT PARAMETERS:
1089 : !
1090 : INTEGER, INTENT(OUT) :: tidx2
1091 : !
1092 : ! !INPUT/OUTPUT PARAMETERS:
1093 : !
1094 : REAL(sp), INTENT(INOUT) :: wgt1
1095 : REAL(sp), INTENT(INOUT) :: wgt2
1096 : INTEGER, INTENT(INOUT) :: RC
1097 : !
1098 : ! !REVISION HISTORY:
1099 : ! 02 Mar 2015 - C. Keller - Initial version
1100 : ! See https://github.com/geoschem/hemco for complete history
1101 : !EOP
1102 : !------------------------------------------------------------------------------
1103 : !BOC
1104 : !
1105 : ! !LOCAL VARIABLES:
1106 : !
1107 : ! Scalars
1108 : INTEGER :: I
1109 : REAL(dp) :: tmpYMDhm
1110 : LOGICAL :: verb
1111 :
1112 : ! Strings
1113 : CHARACTER(LEN=255) :: MSG
1114 : CHARACTER(LEN=255) :: LOC = 'GetIndex2Interp (hcoio_util_mod.F90)'
1115 :
1116 : !=================================================================
1117 : ! GetIndex2Interp begins here
1118 : !=================================================================
1119 :
1120 : ! Verbose mode?
1121 : verb = HCO_IsVerb( HcoState%Config%Err )
1122 :
1123 : ! If the originally wanted datetime was below the available data
1124 : ! range, set all weights to the first index.
1125 0 : IF ( origYMDhm <= availYMDhm(1) ) THEN
1126 0 : tidx2 = tidx1
1127 0 : wgt1 = 1.0_sp
1128 0 : wgt2 = 0.0_sp
1129 :
1130 : ! If the originally wanted datetime is beyond the available data
1131 : ! range, set tidx2 to tidx1 but leave weights in their original
1132 : ! values (-1.0). The reason is that we will attempt to interpolate
1133 : ! between a second file, which is only done if the weights are
1134 : ! negative.
1135 0 : ELSEIF ( origYMDhm >= availYMDhm(nTime) ) THEN
1136 0 : tidx2 = tidx1
1137 :
1138 : ! No interpolation needed if there is a time slices that exactly
1139 : ! matches the (originally) preferred datetime.
1140 : ! NOTE: An Epsilon test is more robust than an equality test
1141 : ! for double-precision variables (bmy, 4/11/17)
1142 0 : ELSEIF ( ABS( origYMDhm - availYMDhm(tidx1) ) < EPSILON ) THEN
1143 0 : tidx2 = tidx1
1144 0 : wgt1 = 1.0_sp
1145 0 : wgt2 = 0.0_sp
1146 :
1147 : ! If we are inside the data range but none of the time slices
1148 : ! matches the preferred datetime, get the second time slices that
1149 : ! shall be used for data interpolation. This not necessarily needs
1150 : ! to be the consecutive time slice. For instance, imagine a data
1151 : ! set that contains montlhly data for years 2005 and 2010. For
1152 : ! Feb 2007, we would want to interpolate between Feb 2005 and Feb
1153 : ! 2010 data. The index tidx1 already points to Feb 2005, but the
1154 : ! upper index tidx2 needs to be set accordingly.
1155 : ELSE
1156 :
1157 : ! Init
1158 0 : tidx2 = -1
1159 :
1160 : ! Search for a time slice in the future that has the same
1161 : ! month/day/hour as currently selected time slice.
1162 : tmpYMDhm = availYMDhm(tidx1)
1163 : DO
1164 : ! Increase by one year
1165 0 : tmpYMDhm = tmpYMDhm + 1.0e8_dp
1166 :
1167 : ! Exit if we are beyond available dates
1168 0 : IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
1169 :
1170 : ! Check if there is a time slice with that date
1171 0 : DO I = tidx1,nTime
1172 0 : IF ( tmpYMDhm == availYMDhm(I) ) THEN
1173 0 : tidx2 = I
1174 0 : EXIT
1175 : ENDIF
1176 : ENDDO
1177 0 : IF ( tidx2 > 0 ) EXIT
1178 : ENDDO
1179 :
1180 : ! Repeat above but now only modify day
1181 0 : IF ( tidx2 < 0 ) THEN
1182 : tmpYMDhm = availYMDhm(tidx1)
1183 : DO
1184 : ! Increase by one day
1185 0 : tmpYMDhm = tmpYMDhm + 1.0e4_dp
1186 :
1187 : ! Exit if we are beyond available dates
1188 0 : IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
1189 :
1190 : ! Check if there is a time slice with that date
1191 0 : DO I = tidx1,nTime
1192 0 : IF ( tmpYMDhm == availYMDhm(I) ) THEN
1193 0 : tidx2 = I
1194 0 : EXIT
1195 : ENDIF
1196 : ENDDO
1197 0 : IF ( tidx2 > 0 ) EXIT
1198 : ENDDO
1199 : ENDIF
1200 :
1201 : ! If all of those tests failed, simply get the next time
1202 : ! slice.
1203 0 : IF ( tidx2 < 0 ) THEN
1204 0 : tidx2 = tidx1 + 1
1205 :
1206 : ! Make sure that tidx2 does not exceed nTime, which is
1207 : ! the number of time slices in the file. This can cause
1208 : ! an out-of-bounds error. (bmy, 3/7/19)
1209 0 : IF ( tidx2 > nTime ) tidx2 = nTime
1210 :
1211 : ! Prompt warning
1212 0 : WRITE(MSG,*) 'Having problems in finding the next time slice ', &
1213 0 : 'to interpolate from, just take the next available ', &
1214 0 : 'slice. Interpolation will be performed from ', &
1215 0 : availYMDhm(tidx1), ' to ', availYMDhm(tidx2), '. Data ', &
1216 0 : 'container: ', TRIM(Lct%Dct%cName)
1217 0 : CALL HCO_WARNING(HcoState%Config%Err, MSG, RC, THISLOC=LOC)
1218 : ENDIF
1219 :
1220 : ! Calculate weights wgt1 and wgt2 to be given to slice 1 and
1221 : ! slice2, respectively.
1222 0 : CALL GetWeights ( availYMDhm(tidx1), availYMDhm(tidx2), origYMDhm, &
1223 0 : wgt1, wgt2 )
1224 :
1225 : ENDIF
1226 :
1227 : ! Return w/ success
1228 0 : RC = HCO_SUCCESS
1229 :
1230 0 : END SUBROUTINE GetIndex2Interp
1231 : !EOC
1232 : !------------------------------------------------------------------------------
1233 : ! Harmonized Emissions Component (HEMCO) !
1234 : !------------------------------------------------------------------------------
1235 : !BOP
1236 : !
1237 : ! !IROUTINE: GetWeights
1238 : !
1239 : ! !DESCRIPTION: Helper function to get the interpolation weights between
1240 : ! two datetime intervals (int1, int2) and for a given time cur.
1241 : !\\
1242 : !\\
1243 : ! !INTERFACE:
1244 : !
1245 0 : SUBROUTINE GetWeights ( int1, int2, cur, wgt1, wgt2 )
1246 : !
1247 : ! !INPUT PARAMETERS:
1248 : !
1249 : REAL(dp), INTENT(IN ) :: int1, int2, cur
1250 : !
1251 : ! !INPUT/OUTPUT PARAMETERS:
1252 : !
1253 : REAL(sp), INTENT( OUT) :: wgt1, wgt2
1254 : !
1255 : ! !REVISION HISTORY:
1256 : ! 04 Mar 2015 - C. Keller - Initial version
1257 : ! See https://github.com/geoschem/hemco for complete history
1258 : !EOP
1259 : !------------------------------------------------------------------------------
1260 : !BOC
1261 : !
1262 : ! !LOCAL VARIABLES:
1263 : !
1264 : REAL(dp) :: diff1, diff2
1265 : REAL(dp) :: jdc, jd1, jd2
1266 :
1267 : !=================================================================
1268 : ! GetWeights begins here!
1269 : !=================================================================
1270 :
1271 : ! Convert dates to Julian dates
1272 0 : jdc = YMDhm2jd ( cur )
1273 0 : jd1 = YMDhm2jd ( int1 )
1274 0 : jd2 = YMDhm2jd ( int2 )
1275 :
1276 : ! Check if outside of range
1277 0 : IF ( jdc <= jd1 ) THEN
1278 0 : wgt1 = 1.0_sp
1279 0 : ELSEIF ( jdc >= jd2 ) THEN
1280 0 : wgt1 = 0.0_sp
1281 : ELSE
1282 0 : diff1 = jd2 - jdc
1283 0 : diff2 = jd2 - jd1
1284 0 : wgt1 = diff1 / diff2
1285 : ENDIF
1286 :
1287 : ! second weight is just complement of wgt1
1288 0 : wgt2 = 1.0_sp - wgt1
1289 :
1290 0 : END SUBROUTINE GetWeights
1291 : !EOC
1292 : !------------------------------------------------------------------------------
1293 : ! Harmonized Emissions Component (HEMCO) !
1294 : !------------------------------------------------------------------------------
1295 : !BOP
1296 : !
1297 : ! !IROUTINE: YMDhm2jd
1298 : !
1299 : ! !DESCRIPTION: returns the julian date of element YMDhm.
1300 : !\\
1301 : !\\
1302 : ! !INTERFACE:
1303 : !
1304 0 : FUNCTION YMDhm2jd ( YMDhm ) RESULT ( jd )
1305 : !
1306 : ! !USES:
1307 : !
1308 : USE HCO_Julday_Mod
1309 : !
1310 : ! !INPUT PARAMETERS:
1311 : !
1312 : REAL(dp), INTENT(IN) :: YMDhm
1313 : !
1314 : ! !INPUT/OUTPUT PARAMETERS:
1315 : !
1316 : REAL(hp) :: jd
1317 : !
1318 : ! !REVISION HISTORY:
1319 : ! 24 Feb 2019 - C. Keller - Initial version
1320 : ! See https://github.com/geoschem/hemco for complete history
1321 : !EOP
1322 : !------------------------------------------------------------------------------
1323 : !BOC
1324 : !
1325 : ! !LOCAL VARIABLES:
1326 : !
1327 : INTEGER :: yr, mt, dy, hr, mn
1328 : REAL(dp) :: utc, day
1329 :
1330 : !=================================================================
1331 : ! YMDh2jd begins here!
1332 : !=================================================================
1333 0 : yr = FLOOR( MOD( YMDhm, 1.0e12_dp ) / 1.0e8_dp )
1334 0 : mt = FLOOR( MOD( YMDhm, 1.0e8_dp ) / 1.0e6_dp )
1335 0 : dy = FLOOR( MOD( YMDhm, 1.0e6_dp ) / 1.0e4_dp )
1336 0 : hr = FLOOR( MOD( YMDhm, 1.0e4_dp ) / 1.0e2_dp )
1337 0 : mn = FLOOR( MOD( YMDhm, 1.0e2_dp ) )
1338 : utc = ( REAL(hr,dp) / 24.0_dp ) + &
1339 : ( REAL(mn,dp) / 1440.0_dp ) + &
1340 0 : ( REAL(0 ,dp) / 86400.0_dp )
1341 0 : day = REAL(dy,dp) + utc
1342 0 : jd = JULDAY( yr, mt, day )
1343 :
1344 0 : END FUNCTION YMDhm2jd
1345 : !EOC
1346 : !------------------------------------------------------------------------------
1347 : ! Harmonized Emissions Component (HEMCO) !
1348 : !------------------------------------------------------------------------------
1349 : !BOP
1350 : !
1351 : ! !IROUTINE: YMDhm2hrs
1352 : !
1353 : ! !DESCRIPTION: returns the hours of element YMDhm. For simplicity, 30 days are
1354 : ! assigned to every month. At the moment, this routine is only called to
1355 : ! determine the time interval between two emission time slices (DeltaT) and
1356 : ! this approximation is good enough.
1357 : !\\
1358 : !\\
1359 : ! !INTERFACE:
1360 : !
1361 0 : FUNCTION YMDhm2hrs ( YMDhm ) RESULT ( hrs )
1362 : !
1363 : ! !INPUT PARAMETERS:
1364 : !
1365 : REAL(dp), INTENT(IN) :: YMDhm
1366 : !
1367 : ! !INPUT/OUTPUT PARAMETERS:
1368 : !
1369 : INTEGER :: hrs
1370 : !
1371 : ! !REVISION HISTORY:
1372 : ! 26 Jan 2015 - C. Keller - Initial version
1373 : ! See https://github.com/geoschem/hemco for complete history
1374 : !EOP
1375 : !------------------------------------------------------------------------------
1376 : !BOC
1377 :
1378 : !=================================================================
1379 : ! YMDh2hrs begins here!
1380 : !=================================================================
1381 : hrs = FLOOR( MOD( YMDhm, 1.0e12_dp ) / 1.0e8_dp ) * 8760 + &
1382 : FLOOR( MOD( YMDhm, 1.0e8_dp ) / 1.0e6_dp ) * 720 + &
1383 : FLOOR( MOD( YMDhm, 1.0e6_dp ) / 1.0e4_dp ) * 24 + &
1384 0 : FLOOR( MOD( YMDhm, 1.0e4_dp ) / 1.0e2_dp )
1385 :
1386 0 : END FUNCTION YMDhm2hrs
1387 : !EOC
1388 : !------------------------------------------------------------------------------
1389 : ! Harmonized Emissions Component (HEMCO) !
1390 : !------------------------------------------------------------------------------
1391 : !BOP
1392 : !
1393 : ! !IROUTINE: Normalize_Area
1394 : !
1395 : ! !DESCRIPTION: Subroutine Normalize\_Area normalizes the given array
1396 : ! by the surface area calculated from the given netCDF file.
1397 : !\\
1398 : !\\
1399 : ! !INTERFACE:
1400 : !
1401 0 : SUBROUTINE Normalize_Area( HcoState, Array, nlon, LatEdge, FN, RC )
1402 : !
1403 : ! !INPUT PARAMETERS:
1404 : !
1405 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1406 : INTEGER, INTENT(IN ) :: nlon ! # of lon midpoints
1407 : REAL(hp), POINTER :: LatEdge(:) ! lat edges
1408 : CHARACTER(LEN=*), INTENT(IN ) :: FN ! filename
1409 : !
1410 : ! !INPUT/OUTPUT PARAMETERS:
1411 : !
1412 : REAL(sp), POINTER :: Array(:,:,:,:) ! Data
1413 : INTEGER, INTENT(INOUT) :: RC ! Return code
1414 : !
1415 : ! !REVISION HISTORY:
1416 : ! 13 Mar 2013 - C. Keller - Initial version
1417 : ! See https://github.com/geoschem/hemco for complete history
1418 : !EOP
1419 : !------------------------------------------------------------------------------
1420 : !BOC
1421 : !
1422 : ! !LOCAL VARIABLES:
1423 : !
1424 : REAL(hp) :: DLAT, AREA
1425 : INTEGER :: NLAT, J
1426 : CHARACTER(LEN=255) :: MSG, LOC
1427 :
1428 : !=================================================================
1429 : ! NORNALIZE_AREA begins here!
1430 : !=================================================================
1431 :
1432 : ! Initialize
1433 0 : LOC = 'NORMALIZE_AREA (hcoio_util_mod.F90 )'
1434 :
1435 : ! Check array size
1436 0 : NLAT = SIZE(LatEdge,1) - 1
1437 :
1438 0 : IF ( SIZE(Array,1) /= nlon ) THEN
1439 0 : MSG = 'Array size does not agree with nlon: ' // TRIM(FN)
1440 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1441 0 : RETURN
1442 : ENDIF
1443 0 : IF ( SIZE(Array,2) /= NLAT ) THEN
1444 0 : MSG = 'Array size does not agree with nlat: ' // TRIM(FN)
1445 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1446 0 : RETURN
1447 : ENDIF
1448 :
1449 : ! Loop over all latitudes
1450 0 : DO J = 1, NLAT
1451 : ! get grid box area in m2 for grid box with lower and upper latitude
1452 : ! llat/ulat: Area = 2 * PI * Re^2 * DLAT / nlon,
1453 : ! where DLAT = abs( sin(ulat) - sin(llat) )
1454 0 : DLAT = ABS( SIN(LatEdge(J+1)*HcoState%Phys%PI_180) &
1455 0 : - SIN(LatEdge(J)*HcoState%Phys%PI_180) )
1456 : AREA = ( 2_hp * HcoState%Phys%PI * DLAT * HcoState%Phys%Re**2 ) &
1457 0 : / REAL(nlon,hp)
1458 :
1459 : ! convert array data to m-2
1460 0 : ARRAY(:,J,:,:) = ARRAY(:,J,:,:) / AREA
1461 : ENDDO
1462 :
1463 : ! Prompt a warning
1464 0 : WRITE(MSG,*) 'No area unit found in ' // TRIM(FN) // ' - convert to m-2!'
1465 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
1466 :
1467 : ! Leave w/ success
1468 0 : RC = HCO_SUCCESS
1469 :
1470 : END SUBROUTINE Normalize_Area
1471 : !EOC
1472 : !------------------------------------------------------------------------------
1473 : ! Harmonized Emissions Component (HEMCO) !
1474 : !------------------------------------------------------------------------------
1475 : !BOP
1476 : !
1477 : ! !IROUTINE: SrcFile_Parse
1478 : !
1479 : ! !DESCRIPTION: Routine SrcFile\_Parse parses the source file name ('ncFile')
1480 : ! of the provided list container Lct. In particular, it searches for tokens
1481 : ! such as $ROOT, $YYYY, etc., within the file name and replaces those values
1482 : ! with the intendend characters. The parsed file name is returned in string
1483 : ! srcFile, while the original file name is retained in Lct.
1484 : !\\
1485 : !\\
1486 : ! It now also checks if the file exists. If the file does not exist and the
1487 : ! file name contains date tokens, it tries to adjust the file name to the
1488 : ! closest available date in the past. The optional flag FUTURE can be used
1489 : ! to denote that the next available file in the future shall be selected,
1490 : ! even if there is a file that exactly matches the preferred date time. This
1491 : ! is useful for interpolation between fields.
1492 : !\\
1493 : !\\
1494 : ! !INTERFACE:
1495 : !
1496 0 : SUBROUTINE SrcFile_Parse ( HcoState, Lct, srcFile, FOUND, RC, &
1497 : Direction, Year )
1498 : !
1499 : ! !USES:
1500 : !
1501 : USE HCO_TIDX_MOD, ONLY : HCO_GetPrefTimeAttr
1502 : USE HCO_TIDX_MOD, ONLY : tIDx_IsInRange
1503 : USE HCO_CLOCK_MOD, ONLY : HcoClock_Get
1504 : USE HCO_CLOCK_MOD, ONLY : Get_LastDayOfMonth
1505 : !
1506 : ! !INPUT PARAMETERS:
1507 : !
1508 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1509 : TYPE(ListCont), POINTER :: Lct ! HEMCO list
1510 : INTEGER, INTENT(IN ), OPTIONAL :: Direction ! Look for file in
1511 : ! future (+1) or
1512 : ! past (-1)
1513 : INTEGER, INTENT(IN ), OPTIONAL :: Year ! To use fixed year
1514 : !
1515 : ! !OUTPUT PARAMETERS:
1516 : !
1517 : CHARACTER(LEN=*), INTENT( OUT) :: srcFile ! output string
1518 : LOGICAL, INTENT( OUT) :: FOUND ! Does file exist?
1519 : !
1520 : ! !INPUT/OUTPUT PARAMETERS:
1521 : !
1522 : INTEGER, INTENT(INOUT) :: RC ! return code
1523 : !
1524 : ! !REVISION HISTORY:
1525 : ! 01 Oct 2014 - C. Keller - Initial version
1526 : ! See https://github.com/geoschem/hemco for complete history
1527 : !EOP
1528 : !------------------------------------------------------------------------------
1529 : !BOC
1530 : !
1531 : ! !LOCAL VARIABLES:
1532 : !
1533 : INTEGER :: INC, CNT, TYPCNT, TYP, NEWTYP
1534 : INTEGER :: prefYr, prefMt, prefDy, prefHr, prefMn
1535 : INTEGER :: origYr, origMt, origDy, origHr
1536 : LOGICAL :: hasFile, hasYr, hasMt, hasDy, hasHr
1537 : LOGICAL :: nextTyp
1538 : CHARACTER(LEN=1023) :: MSG, LOC
1539 : CHARACTER(LEN=1023) :: srcFileOrig
1540 :
1541 : ! maximum # of iterations for file search
1542 : INTEGER, PARAMETER :: MAXIT = 10000
1543 :
1544 : !=================================================================
1545 : ! SrcFile_Parse
1546 : !=================================================================
1547 :
1548 : ! Initialize
1549 0 : LOC = 'SrcFile_Parse (HCOIO_UTIL_MOD.F90)'
1550 0 : RC = HCO_SUCCESS
1551 0 : found = .FALSE.
1552 0 : srcFile = Lct%Dct%Dta%ncFile
1553 :
1554 : ! verbose mode
1555 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1556 0 : WRITE(MSG,*) 'Parsing source file and replacing tokens'
1557 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1558 : ENDIF
1559 :
1560 : ! Get preferred dates (to be passed to parser)
1561 : CALL HCO_GetPrefTimeAttr ( HcoState, Lct, &
1562 0 : prefYr, prefMt, prefDy, prefHr, prefMn, RC )
1563 0 : IF ( RC /= HCO_SUCCESS ) THEN
1564 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
1565 0 : RETURN
1566 : ENDIF
1567 :
1568 : ! Make sure dates are not negative
1569 0 : IF ( prefYr <= 0 ) THEN
1570 0 : CALL HcoClock_Get( HcoState%Clock, cYYYY = prefYr, RC = RC )
1571 0 : IF ( RC /= HCO_SUCCESS ) THEN
1572 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
1573 0 : RETURN
1574 : ENDIF
1575 : ENDIF
1576 0 : IF ( prefMt <= 0 ) THEN
1577 0 : CALL HcoClock_Get( HcoState%Clock, cMM = prefMt, RC = RC )
1578 0 : IF ( RC /= HCO_SUCCESS ) THEN
1579 0 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
1580 0 : RETURN
1581 : ENDIF
1582 : ENDIF
1583 0 : IF ( prefDy <= 0 ) THEN
1584 0 : CALL HcoClock_Get( HcoState%Clock, cDD = prefDy, RC = RC )
1585 0 : IF ( RC /= HCO_SUCCESS ) THEN
1586 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
1587 0 : RETURN
1588 : ENDIF
1589 : ENDIF
1590 0 : IF ( prefHr < 0 ) THEN
1591 0 : CALL HcoClock_Get( HcoState%Clock, cH = prefHr, RC = RC )
1592 0 : IF ( RC /= HCO_SUCCESS ) THEN
1593 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
1594 0 : RETURN
1595 : ENDIF
1596 : ENDIF
1597 :
1598 : ! Eventually replace default preferred year with specified one
1599 0 : IF ( PRESENT(Year) ) prefYr = Year
1600 :
1601 : ! Call the parser
1602 : CALL HCO_CharParse ( HcoState%Config, srcFile, prefYr, prefMt, &
1603 0 : prefDy, prefHr, prefMn, RC )
1604 0 : IF ( RC /= HCO_SUCCESS ) THEN
1605 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
1606 0 : RETURN
1607 : ENDIF
1608 0 : srcFileOrig = TRIM(srcFile)
1609 :
1610 : ! Check if file exists
1611 0 : INQUIRE( FILE=TRIM(srcFile), EXIST=HasFile )
1612 :
1613 : ! If the direction flag is on, force HasFile to be false.
1614 0 : IF ( PRESENT(Direction) ) THEN
1615 0 : IF ( Direction /= 0 ) HasFile = .FALSE.
1616 : ENDIF
1617 :
1618 : !-----------------------------------------------------------------------
1619 : ! If this is a HEMCO dry-run simulation, then do not enter the loop
1620 : ! where we will attempt to go back in time until a file is found.
1621 : ! For the dry-run we need to report all files, even missing.
1622 : ! This fixes Github issue geoschem/geos-chem #312. (bmy, 6/9/20)
1623 : !-----------------------------------------------------------------------
1624 0 : IF ( HcoState%Options%isDryRun ) THEN
1625 :
1626 : ! Make sure that the year is not 1, this indicates that the
1627 : ! preferred year is outside of the years specified in the
1628 : ! time range settings in the configuration file, and will
1629 : ! lead to files with a year of "0001" in the path.
1630 : ! (bmy, 6/9/20)
1631 0 : IF ( prefyr == 1 ) THEN
1632 : MSG = 'Cannot find file for current simulation time: ' // &
1633 : TRIM(srcFile) // ' - Cannot get field ' // &
1634 : TRIM(Lct%Dct%cName) // '. Please check file name ' // &
1635 0 : 'and time (incl. time range flag) in the config. file'
1636 0 : CALL HCO_ERROR( MSG, RC )
1637 0 : RETURN
1638 : ENDIF
1639 :
1640 : ! Otherwise return with success
1641 0 : RC = HCO_SUCCESS
1642 0 : Found = HasFile
1643 0 : RETURN
1644 : ENDIF
1645 :
1646 : ! If file does not exist, check if we can adjust prefYr, prefMt, etc.
1647 0 : IF ( .NOT. HasFile .AND. Lct%Dct%DctType /= HCO_CFLAG_EXACT ) THEN
1648 :
1649 : ! Check if any token exist
1650 0 : HasYr = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'YYYY') > 0 )
1651 0 : HasMt = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'MM' ) > 0 )
1652 0 : HasDy = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'DD' ) > 0 )
1653 0 : HasHr = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'HH' ) > 0 )
1654 :
1655 : ! Search for file
1656 0 : IF ( HasYr .OR. HasMt .OR. HasDy .OR. HasHr ) THEN
1657 :
1658 : ! Date increments
1659 0 : INC = -1
1660 0 : IF ( PRESENT(Direction) ) THEN
1661 0 : INC = Direction
1662 : ENDIF
1663 :
1664 : ! Initialize counters
1665 0 : CNT = 0
1666 :
1667 : ! Type is the update type (see below)
1668 0 : TYP = 0
1669 :
1670 : ! Mirror preferred variables
1671 0 : origYr = prefYr
1672 0 : origMt = prefMt
1673 0 : origDy = prefDy
1674 0 : origHr = prefHr
1675 :
1676 : ! Do until file is found or counter exceeds threshold
1677 0 : DO WHILE ( .NOT. HasFile )
1678 :
1679 : ! Inrease counter
1680 0 : CNT = CNT + 1
1681 0 : IF ( CNT > MAXIT ) EXIT
1682 :
1683 : ! Increase update type if needed:
1684 0 : nextTyp = .FALSE.
1685 :
1686 : ! Type 0: Initialization
1687 0 : IF ( TYP == 0 ) THEN
1688 : nextTyp = .TRUE.
1689 : ! Type 1: update hour only
1690 0 : ELSEIF ( TYP == 1 .AND. TYPCNT > 24 ) THEN
1691 : nextTyp = .TRUE.
1692 : ! Type 2: update day only
1693 0 : ELSEIF ( TYP == 2 .AND. TYPCNT > 31 ) THEN
1694 : nextTyp = .TRUE.
1695 : ! Type 3: update month only
1696 0 : ELSEIF ( TYP == 3 .AND. TYPCNT > 12 ) THEN
1697 : nextTyp = .TRUE.
1698 : ! Type 4: update year only
1699 0 : ELSEIF ( TYP == 4 .AND. TYPCNT > 300 ) THEN
1700 : nextTyp = .TRUE.
1701 : ! Type 5: update hour and day
1702 0 : ELSEIF ( TYP == 5 .AND. TYPCNT > 744 ) THEN
1703 : nextTyp = .TRUE.
1704 : ! Type 6: update day and month
1705 0 : ELSEIF ( TYP == 6 .AND. TYPCNT > 372 ) THEN
1706 : nextTyp = .TRUE.
1707 : ! Type 7: update month and year
1708 0 : ELSEIF ( TYP == 7 .AND. TYPCNT > 3600 ) THEN
1709 : EXIT
1710 : ENDIF
1711 :
1712 : ! Get next type
1713 : IF ( nextTyp ) THEN
1714 0 : NEWTYP = -1
1715 0 : IF ( hasHr .AND. TYP < 1 ) THEN
1716 : NEWTYP = 1
1717 0 : ELSEIF ( hasDy .AND. TYP < 2 ) THEN
1718 : NEWTYP = 2
1719 0 : ELSEIF ( hasMt .AND. TYP < 3 ) THEN
1720 : NEWTYP = 3
1721 0 : ELSEIF ( hasYr .AND. TYP < 4 ) THEN
1722 : NEWTYP = 4
1723 : ELSEIF ( hasDy .AND. TYP < 2 ) THEN
1724 : NEWTYP = 5
1725 : ELSEIF ( hasDy .AND. TYP < 2 ) THEN
1726 : NEWTYP = 6
1727 : ELSEIF ( hasDy .AND. TYP < 2 ) THEN
1728 : NEWTYP = 7
1729 : ENDIF
1730 :
1731 : ! Exit if no other type found
1732 : IF ( NEWTYP < 0 ) EXIT
1733 :
1734 : ! This is the new type, reset type counter
1735 0 : TYP = NEWTYP
1736 0 : TYPCNT = 0
1737 :
1738 : ! Make sure we reset all values
1739 0 : prefYr = origYr
1740 0 : prefMt = origMt
1741 0 : prefDy = origDy
1742 0 : prefHr = origHr
1743 :
1744 : ENDIF
1745 :
1746 : ! Update preferred datetimes
1747 0 : SELECT CASE ( TYP )
1748 : ! Adjust hour only
1749 : CASE ( 1 )
1750 0 : prefHr = prefHr + INC
1751 : ! Adjust day only
1752 : CASE ( 2 )
1753 0 : prefDy = prefDy + INC
1754 : ! Adjust month only
1755 : CASE ( 3 )
1756 0 : prefMt = prefMt + INC
1757 : ! Adjust year only
1758 : CASE ( 4 )
1759 0 : prefYr = prefYr + INC
1760 : ! Adjust hour and day
1761 : CASE ( 5 )
1762 0 : prefHr = prefHr + INC
1763 0 : IF ( MOD(TYPCNT,24) == 0 ) prefDy = prefDy + INC
1764 : ! Adjust day and month
1765 : CASE ( 6 )
1766 0 : prefDy = prefDy + INC
1767 0 : IF ( MOD(TYPCNT,31) == 0 ) prefMt = prefMt + INC
1768 : ! Adjust month and year
1769 : CASE ( 7 )
1770 0 : prefMt = prefMt + INC
1771 0 : IF ( MOD(TYPCNT,12) == 0 ) prefYr = prefYr + INC
1772 : CASE DEFAULT
1773 0 : EXIT
1774 : END SELECT
1775 :
1776 : ! Check if we need to adjust a year/month/day/hour
1777 0 : IF ( prefHr < 0 ) THEN
1778 0 : prefHr = 23
1779 0 : prefDy = prefDy - 1
1780 : ENDIF
1781 0 : IF ( prefHr > 23 ) THEN
1782 0 : prefHr = 0
1783 0 : prefDy = prefDy + 1
1784 : ENDIF
1785 0 : IF ( prefDy < 1 ) THEN
1786 0 : prefDy = 31
1787 0 : prefMt = prefMt - 1
1788 : ENDIF
1789 0 : IF ( prefDy > 31 ) THEN
1790 0 : prefDy = 1
1791 0 : prefMt = prefMt + 1
1792 : ENDIF
1793 0 : IF ( prefMt < 1 ) THEN
1794 0 : prefMt = 12
1795 0 : prefYr = prefYr - 1
1796 : ENDIF
1797 0 : IF ( prefMt > 12 ) THEN
1798 0 : prefMt = 1
1799 0 : prefYr = prefYr + 1
1800 : ENDIF
1801 :
1802 : ! Make sure day does not exceed max. number of days in this month
1803 0 : prefDy = MIN( prefDy, Get_LastDayOfMonth( prefMt, prefYr ) )
1804 :
1805 : ! Mirror original file
1806 0 : srcFile = Lct%Dct%Dta%ncFile
1807 :
1808 : ! Call the parser with adjusted values
1809 : CALL HCO_CharParse ( HcoState%Config, srcFile, prefYr, &
1810 0 : prefMt, prefDy, prefHr, prefMn, RC )
1811 0 : IF ( RC /= HCO_SUCCESS ) THEN
1812 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
1813 0 : RETURN
1814 : ENDIF
1815 :
1816 : ! Check if this file exists
1817 0 : INQUIRE( FILE=TRIM(srcFile), EXIST=HasFile )
1818 :
1819 : ! Update counter
1820 0 : TYPCNT = TYPCNT + 1
1821 : ENDDO
1822 : ENDIF
1823 : ENDIF
1824 :
1825 : ! Additional check for data with a given range: make sure that the selected
1826 : ! field is not outside of the given range
1827 0 : IF ( HasFile .AND. ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) ) THEN
1828 0 : HasFile = TIDX_IsInRange ( Lct, prefYr, prefMt, prefDy, prefHr )
1829 : ENDIF
1830 :
1831 : ! Restore original source file name and date to avoid confusion in log file
1832 0 : IF ( .not. HasFile ) THEN
1833 0 : srcFile = Trim(srcFileOrig)
1834 : ENDIF
1835 :
1836 : ! Return variable
1837 0 : FOUND = HasFile
1838 :
1839 : ! Return w/ success
1840 0 : RC = HCO_SUCCESS
1841 :
1842 0 : END SUBROUTINE SrcFile_Parse
1843 : !EOC
1844 : !------------------------------------------------------------------------------
1845 : ! Harmonized Emissions Component (HEMCO) !
1846 : !------------------------------------------------------------------------------
1847 : !BOP
1848 : !
1849 : ! !IROUTINE: SigmaMidToEdges
1850 : !
1851 : ! !DESCRIPTION: Helper routine to interpolate sigma mid point values to edges.
1852 : ! A simple linear interpolation is performed.
1853 : !\\
1854 : !\\
1855 : ! !INTERFACE:
1856 : !
1857 0 : SUBROUTINE SigmaMidToEdges ( HcoState, SigMid, SigEdge, RC )
1858 : !
1859 : ! !INPUT PARAMETERS:
1860 : !
1861 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
1862 : REAL(hp), POINTER :: SigMid(:,:,:) ! sigma levels
1863 : !
1864 : ! !OUTPUT PARAMETERS:
1865 : !
1866 : REAL(hp), POINTER :: SigEdge(:,:,:) ! sigma edges
1867 : INTEGER, INTENT( OUT) :: RC ! return code
1868 : !
1869 : ! !REVISION HISTORY:
1870 : ! 03 Oct 2013 - C. Keller - Initial version
1871 : ! See https://github.com/geoschem/hemco for complete history
1872 : !EOP
1873 : !------------------------------------------------------------------------------
1874 : !BOC
1875 : !
1876 : ! !LOCAL VARIABLES:
1877 : !
1878 : INTEGER :: L, AS
1879 : INTEGER :: nx, ny, nz
1880 : CHARACTER(LEN=255) :: MSG
1881 : CHARACTER(LEN=255) :: LOC = 'SigmaMidToEdges (hcoio_util_mod.F90)'
1882 :
1883 : !=================================================================
1884 : ! SigmaMidToEdges begins here!
1885 : !=================================================================
1886 :
1887 : ! Allocate space as required
1888 0 : nx = SIZE(SigMid,1)
1889 0 : ny = SIZE(SigMid,2)
1890 0 : nz = SIZE(SigMid,3)
1891 0 : IF ( ASSOCIATED(SigEdge) ) DEALLOCATE(SigEdge)
1892 0 : ALLOCATE(SigEdge(nx,ny,nz+1),STAT=AS)
1893 0 : IF ( AS/=0 ) THEN
1894 : CALL HCO_ERROR( 'Allocate SigEdge', RC, &
1895 0 : THISLOC=LOC )
1896 0 : RETURN
1897 : ENDIF
1898 0 : SigEdge = 0.0_hp
1899 :
1900 : ! Calculate sigma edges by linear interpolation (symmetric mid-points)
1901 0 : DO L = 1, nz-1
1902 0 : SigEdge(:,:,L+1) = ( SigMid(:,:,L) + SigMid(:,:,L+1) ) / 2.0_hp
1903 : ENDDO
1904 :
1905 : ! Get outermost values:
1906 0 : SigEdge(:,:,1 ) = SigMid(:,:,1 ) - ( SigEdge(:,:,2) - SigMid(:,:,1) )
1907 0 : SigEdge(:,:,nz+1) = SigMid(:,:,nz) + ( SigMid(:,:,nz) - SigEdge(:,:,nz) )
1908 :
1909 : ! Return w/ success
1910 0 : RC = HCO_SUCCESS
1911 :
1912 : END SUBROUTINE SigmaMidToEdges
1913 : !EOC
1914 : !------------------------------------------------------------------------------
1915 : ! Harmonized Emissions Component (HEMCO) !
1916 : !------------------------------------------------------------------------------
1917 : !BOP
1918 : !
1919 : ! !IROUTINE: CheckMissVal
1920 : !
1921 : ! !DESCRIPTION: Checks for missing values in the passed array. Missing values
1922 : ! of base emissions and masks are set to 0, missing values of scale factors
1923 : ! are set to 1.
1924 : !\\
1925 : ! !INTERFACE:
1926 : !
1927 0 : SUBROUTINE CheckMissVal ( Lct, Arr )
1928 : !
1929 : ! !INPUT PARAMETERS:
1930 : !
1931 : TYPE(ListCont), POINTER :: Lct
1932 : REAL(sp), POINTER :: Arr(:,:,:,:)
1933 : !
1934 : ! !REVISION HISTORY:
1935 : ! 04 Mar 2015 - C. Keller - Initial version
1936 : ! See https://github.com/geoschem/hemco for complete history
1937 : !EOP
1938 : !------------------------------------------------------------------------------
1939 : !BOC
1940 : !
1941 : ! !LOCAL VARIABLES:
1942 : !
1943 : !=================================================================
1944 : ! CheckMissVal begins here!
1945 : !=================================================================
1946 :
1947 : ! Error trap
1948 0 : IF ( .NOT. ASSOCIATED(Arr) ) RETURN
1949 :
1950 0 : IF ( ANY(Arr == HCO_MISSVAL) ) THEN
1951 : ! Base emissions
1952 0 : IF ( Lct%Dct%DctType == HCO_DCTTYPE_BASE ) THEN
1953 0 : WHERE(Arr == HCO_MISSVAL) Arr = 0.0_sp
1954 : ! Scale factor
1955 0 : ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_SCAL ) THEN
1956 0 : WHERE(Arr == HCO_MISSVAL) Arr = 1.0_sp
1957 : ! Mask
1958 0 : ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
1959 0 : WHERE(Arr == HCO_MISSVAL) Arr = 0.0_sp
1960 : ENDIF
1961 : ENDIF
1962 :
1963 : END SUBROUTINE CheckMissVal
1964 : !EOC
1965 : !------------------------------------------------------------------------------
1966 : ! Harmonized Emissions Component (HEMCO) !
1967 : !------------------------------------------------------------------------------
1968 : !BOP
1969 : !
1970 : ! !IROUTINE: GetArbDimIndex
1971 : !
1972 : ! !DESCRIPTION: Subroutine GetArbDimIndex returns the index of the arbitrary
1973 : ! file dimension. -1 if no such dimension is defined.
1974 : !\\
1975 : ! !INTERFACE:
1976 : !
1977 0 : SUBROUTINE GetArbDimIndex( HcoState, Lun, Lct, ArbIdx, RC )
1978 : !
1979 : ! !USES:
1980 : !
1981 : USE HCO_m_netcdf_io_checks
1982 : USE HCO_m_netcdf_io_get_dimlen
1983 : USE HCO_ExtList_Mod, ONLY : GetExtOpt
1984 : !
1985 : ! !INPUT PARAMETERS:
1986 : !
1987 : TYPE(HCO_State), POINTER :: HcoState
1988 : INTEGER, INTENT(IN ) :: Lun
1989 : TYPE(ListCont), POINTER :: Lct
1990 : !
1991 : ! !OUTPUT PARAMETERS:
1992 : !
1993 : INTEGER, INTENT( OUT) :: ArbIdx
1994 : INTEGER, INTENT( OUT) :: RC
1995 : !
1996 : ! !REVISION HISTORY:
1997 : ! 22 Sep 2015 - C. Keller - Initial version
1998 : ! See https://github.com/geoschem/hemco for complete history
1999 : !EOP
2000 : !------------------------------------------------------------------------------
2001 : !BOC
2002 : !
2003 : ! !LOCAL VARIABLES:
2004 : !
2005 : INTEGER :: TargetVal, nVal
2006 : LOGICAL :: Found
2007 : CHARACTER(LEN=255) :: ArbDimVal
2008 : CHARACTER(LEN=511) :: MSG
2009 : CHARACTER(LEN=255) :: LOC = 'GetArbDimIndex (hcoio_util_mod.F90)'
2010 :
2011 : !=================================================================
2012 : ! GetArbDimIndex
2013 : !=================================================================
2014 :
2015 : ! Assume success until otherwise
2016 0 : RC = HCO_SUCCESS
2017 :
2018 : ! Init
2019 0 : ArbIdx = -1
2020 0 : IF ( TRIM(Lct%Dct%Dta%ArbDimName) == 'none' ) RETURN
2021 :
2022 : ! Check if variable exists
2023 0 : Found = Ncdoes_Dim_Exist ( Lun, TRIM(Lct%Dct%Dta%ArbDimName) )
2024 0 : IF ( .NOT. Found ) THEN
2025 : MSG = 'Cannot read dimension ' // TRIM(Lct%Dct%Dta%ArbDimName) &
2026 : // ' from file ' // &
2027 0 : TRIM(Lct%Dct%Dta%ncFile)
2028 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2029 0 : RETURN
2030 : ENDIF
2031 :
2032 : ! Get dimension length
2033 0 : CALL Ncget_Dimlen ( Lun, TRIM(Lct%Dct%Dta%ArbDimName), nVal )
2034 :
2035 : ! Get value to look for. This is archived in variable ArbDimVal.
2036 : ! Eventually need to extract value from HEMCO settings
2037 0 : ArbDimVal = TRIM(Lct%Dct%Dta%ArbDimVal)
2038 :
2039 : ! If string starts with a number, evaluate value directly
2040 : IF ( ArbDimVal(1:1) == '0' .OR. &
2041 : ArbDimVal(1:1) == '1' .OR. &
2042 : ArbDimVal(1:1) == '2' .OR. &
2043 : ArbDimVal(1:1) == '3' .OR. &
2044 : ArbDimVal(1:1) == '4' .OR. &
2045 : ArbDimVal(1:1) == '5' .OR. &
2046 : ArbDimVal(1:1) == '6' .OR. &
2047 : ArbDimVal(1:1) == '7' .OR. &
2048 0 : ArbDimVal(1:1) == '8' .OR. &
2049 : ArbDimVal(1:1) == '9' ) THEN
2050 0 : READ(ArbDimVal,*) TargetVal
2051 :
2052 : ! Otherwise, assume this is a HEMCO option (including a token)
2053 : ELSE
2054 0 : IF ( ArbDimVal(1:1) == '$' ) ArbDimVal = ArbDimVal(2:LEN(ArbDimVal))
2055 : CALL GetExtOpt ( HcoState%Config, ExtNr=-999, &
2056 : OptName=TRIM(ArbDimVal), &
2057 0 : OptValInt=TargetVal, FOUND=Found, RC=RC )
2058 0 : IF ( RC /= HCO_SUCCESS ) THEN
2059 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
2060 0 : RETURN
2061 : ENDIF
2062 0 : IF ( .NOT. Found ) THEN
2063 0 : WRITE(MSG,*) 'Cannot evaluate additional dimension value ', &
2064 0 : TRIM(ArbDimVal), '. This does not seem to be a number nor ', &
2065 0 : 'a HEMCO token/setting. This error happened when evaluating ', &
2066 0 : 'dimension ', TRIM(Lct%Dct%Dta%ArbDimName), ' belonging to ', &
2067 0 : 'file ', TRIM(Lct%Dct%Dta%ncFile)
2068 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2069 0 : RETURN
2070 : ENDIF
2071 : ENDIF
2072 :
2073 0 : IF ( TargetVal > nVal ) THEN
2074 0 : WRITE(MSG,*) 'Desired dimension value ', TargetVal, &
2075 0 : ' exceeds corresponding dimension length on that file: ', nVal, &
2076 0 : 'This error happened when evaluating ', &
2077 0 : 'dimension ', TRIM(Lct%Dct%Dta%ArbDimName), ' belonging to ', &
2078 0 : 'file ', TRIM(Lct%Dct%Dta%ncFile)
2079 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2080 0 : RETURN
2081 :
2082 : ELSE
2083 0 : ArbIdx = TargetVal
2084 : ENDIF
2085 :
2086 : ! Verbose
2087 0 : IF ( HcoState%amIRoot .AND. HCO_IsVerb( HcoState%Config%Err ) ) THEN
2088 0 : WRITE(MSG,*) 'Additional dimension ', TRIM(Lct%Dct%Dta%ArbDimName), &
2089 0 : ' in ', TRIM(Lct%Dct%Dta%ncFile), ': use index ', &
2090 0 : ArbIdx, ' (set: ', Lct%Dct%Dta%ArbDimVal, ')'
2091 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2092 : ENDIF
2093 :
2094 : ! Return w/ success
2095 0 : RC = HCO_SUCCESS
2096 :
2097 : END SUBROUTINE GetArbDimIndex
2098 : !EOC
2099 : #endif
2100 : !------------------------------------------------------------------------------
2101 : ! Harmonized Emissions Component (HEMCO) !
2102 : !------------------------------------------------------------------------------
2103 : !BOP
2104 : !
2105 : ! !IROUTINE: HCOIO_ReadOther
2106 : !
2107 : ! !DESCRIPTION: Subroutine HCOIO\_ReadOther is a wrapper routine to
2108 : ! read data from sources other than netCDF.
2109 : !\\
2110 : !\\
2111 : ! If a file name is given (ending with '.txt'), the data are assumed
2112 : ! to hold country-specific values (e.g. diurnal scale factors). In all
2113 : ! other cases, the data is directly read from the configuration file
2114 : ! (scalars).
2115 : !\\
2116 : !\\
2117 : ! !INTERFACE:
2118 : !
2119 0 : SUBROUTINE HCOIO_ReadOther( HcoState, Lct, RC )
2120 : !
2121 : ! !USES:
2122 : !
2123 : !
2124 : ! !INPUT PARAMTERS:
2125 : !
2126 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
2127 : !
2128 : ! !INPUT/OUTPUT PARAMETERS:
2129 : !
2130 : TYPE(ListCont), POINTER :: Lct
2131 : INTEGER, INTENT(INOUT) :: RC
2132 : !
2133 : ! !REVISION HISTORY:
2134 : ! 22 Dec 2014 - C. Keller: Initial version
2135 : ! See https://github.com/geoschem/hemco for complete history
2136 : !EOP
2137 : !------------------------------------------------------------------------------
2138 : !BOC
2139 : !
2140 : ! !LOCAL VARIABLES:
2141 : !
2142 : CHARACTER(LEN=255) :: MSG, LOC
2143 :
2144 : !======================================================================
2145 : ! HCOIO_ReadOther begins here
2146 : !======================================================================
2147 0 : LOC = 'HCOIO_ReadOther (HCOIO_UTIL_MOD.F90)'
2148 :
2149 : ! Error check: data must be in local time
2150 0 : IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
2151 : MSG = 'Cannot read data from file that is not in local time: ' // &
2152 0 : TRIM(Lct%Dct%cName)
2153 0 : CALL HCO_ERROR( MSG, RC, THISLOC='HCOIO_ReadOther (hcoio_dataread_mod.F90)' )
2154 0 : RETURN
2155 : ENDIF
2156 :
2157 : ! Read an ASCII file as country values
2158 0 : IF ( INDEX( TRIM(Lct%Dct%Dta%ncFile), '.txt' ) > 0 ) THEN
2159 0 : CALL HCOIO_ReadCountryValues( HcoState, Lct, RC )
2160 0 : IF ( RC /= HCO_SUCCESS ) THEN
2161 0 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
2162 0 : RETURN
2163 : ENDIF
2164 :
2165 : ! Directly read from configuration file otherwise
2166 : ELSE
2167 0 : CALL HCOIO_ReadFromConfig( HcoState, Lct, RC )
2168 0 : IF ( RC /= HCO_SUCCESS ) THEN
2169 0 : CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
2170 0 : RETURN
2171 : ENDIF
2172 : ENDIF
2173 :
2174 : ! Return w/ success
2175 0 : RC = HCO_SUCCESS
2176 :
2177 : END SUBROUTINE HCOIO_ReadOther
2178 : !EOC
2179 : !------------------------------------------------------------------------------
2180 : ! Harmonized Emissions Component (HEMCO) !
2181 : !------------------------------------------------------------------------------
2182 : !BOP
2183 : !
2184 : ! !IROUTINE: HCOIO_ReadCountryValues
2185 : !
2186 : ! !DESCRIPTION: Subroutine HCOIO\_ReadCountryValues
2187 : !\\
2188 : !\\
2189 : ! !INTERFACE:
2190 : !
2191 0 : SUBROUTINE HCOIO_ReadCountryValues ( HcoState, Lct, RC )
2192 : !
2193 : ! !USES:
2194 : !
2195 : USE HCO_inquireMod, ONLY : findFreeLUN
2196 : USE HCO_CHARTOOLS_MOD, ONLY : HCO_CMT, HCO_SPC, NextCharPos
2197 : USE HCO_EmisList_Mod, ONLY : HCO_GetPtr
2198 : USE HCO_FileData_Mod, ONLY : FileData_ArrCheck
2199 : !
2200 : ! !INPUT PARAMTERS:
2201 : !
2202 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
2203 : !
2204 : ! !INPUT/OUTPUT PARAMETERS:
2205 : !
2206 : TYPE(ListCont), POINTER :: Lct
2207 : INTEGER, INTENT(INOUT) :: RC
2208 : !
2209 : ! !REVISION HISTORY:
2210 : ! 22 Dec 2014 - C. Keller: Initial version
2211 : ! See https://github.com/geoschem/hemco for complete history
2212 : !EOP
2213 : !------------------------------------------------------------------------------
2214 : !BOC
2215 : !
2216 : ! !LOCAL VARIABLES:
2217 : !
2218 : INTEGER :: IUFILE, IOS
2219 : INTEGER :: ID1, ID2, I, NT, CID, NLINE
2220 0 : REAL(sp), POINTER :: CNTR(:,:)
2221 0 : INTEGER, ALLOCATABLE :: CIDS(:,:)
2222 0 : REAL(hp), POINTER :: Vals(:)
2223 : LOGICAL :: Verb
2224 : CHARACTER(LEN=2047) :: LINE
2225 : CHARACTER(LEN=255) :: MSG, DUM, CNT
2226 : CHARACTER(LEN=255) :: LOC = 'HCOIO_ReadCountryValues (hcoio_util_mod.F90)'
2227 :
2228 : !======================================================================
2229 : ! HCOIO_ReadCountryValues begins here
2230 : !======================================================================
2231 :
2232 : ! Init
2233 0 : CNTR => NULL()
2234 0 : Vals => NULL()
2235 :
2236 : ! verbose mode?
2237 0 : Verb = HCO_IsVerb( HcoState%Config%Err )
2238 :
2239 : ! Verbose
2240 0 : IF ( Verb ) THEN
2241 0 : MSG = 'Use country-specific values for ' // TRIM(Lct%Dct%cName)
2242 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2243 0 : MSG = '- Source file: ' // TRIM(Lct%Dct%Dta%ncFile)
2244 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2245 : ENDIF
2246 :
2247 : ! Open file
2248 0 : IUFILE = FindFreeLun()
2249 0 : OPEN ( IUFILE, FILE=TRIM( Lct%Dct%Dta%ncFile ), STATUS='OLD', IOSTAT=IOS )
2250 0 : IF ( IOS /= 0 ) THEN
2251 0 : MSG = 'Cannot open ' // TRIM(Lct%Dct%Dta%ncFile)
2252 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2253 0 : RETURN
2254 : ENDIF
2255 :
2256 : ! Repeat for every line
2257 : NLINE = 0
2258 : DO
2259 :
2260 : ! Read line
2261 0 : READ( IUFILE, '(a)', IOSTAT=IOS ) LINE
2262 :
2263 : ! End of file?
2264 0 : IF ( IOS < 0 ) EXIT
2265 :
2266 : ! Error?
2267 0 : IF ( IOS > 0 ) THEN
2268 0 : MSG = 'Error reading ' // TRIM(Lct%Dct%Dta%ncFile)
2269 0 : MSG = TRIM(MSG) // ' - last valid line: ' // TRIM(LINE)
2270 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2271 0 : RETURN
2272 : ENDIF
2273 :
2274 : ! Skip commented lines and/or empty lines
2275 0 : IF ( TRIM(LINE) == '' ) CYCLE
2276 0 : IF ( LINE(1:1) == HCO_CMT ) CYCLE
2277 :
2278 : ! First (valid) line holds the name of the mask container
2279 0 : IF ( NLINE == 0 ) THEN
2280 :
2281 : ! Get pointer to mask. Convert to integer
2282 0 : CALL HCO_GetPtr( HcoState, TRIM(LINE), CNTR, RC )
2283 0 : IF ( RC /= HCO_SUCCESS ) THEN
2284 0 : CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
2285 0 : RETURN
2286 : ENDIF
2287 0 : ALLOCATE( CIDS(HcoState%NX, HcoState%NY), STAT=IOS )
2288 0 : IF ( IOS /= 0 ) THEN
2289 0 : CALL HCO_ERROR( 'Cannot allocate CIDS', RC, THISLOC=LOC )
2290 0 : RETURN
2291 : ENDIF
2292 0 : CIDS = NINT(CNTR)
2293 :
2294 : ! Verbose
2295 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
2296 0 : MSG = '- Use ID mask ' // TRIM(LINE)
2297 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2298 : ENDIF
2299 :
2300 : ! Go to next line
2301 0 : NLINE = NLINE + 1
2302 0 : CYCLE
2303 : ENDIF
2304 :
2305 : ! Get first space character to skip country name.
2306 : ! We assume here that a country name is given right at the
2307 : ! beginning of the line, e.g. 'USA 744 1.05/1.02/...'
2308 0 : ID1 = NextCharPos( LINE, HCO_SPC )
2309 0 : CNT = LINE(1:ID1)
2310 :
2311 : ! Get country ID
2312 0 : DO I = ID1, LEN(LINE)
2313 0 : IF ( LINE(I:I) /= HCO_SPC ) EXIT
2314 : ENDDO
2315 0 : ID1 = I
2316 0 : ID2 = NextCharPos( LINE, HCO_SPC, START=ID1 )
2317 :
2318 0 : IF ( ID2 >= LEN(LINE) .OR. ID2 < 0 ) THEN
2319 0 : MSG = 'Cannot extract country ID from: ' // TRIM(LINE)
2320 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2321 0 : RETURN
2322 : ENDIF
2323 0 : DUM = LINE(ID1:ID2)
2324 0 : READ( DUM, * ) CID
2325 :
2326 : ! Extract data values
2327 0 : ID1 = ID2+1
2328 0 : ID2 = LEN(LINE)
2329 0 : LINE = LINE(ID1:ID2)
2330 0 : CALL GetDataVals( HcoState, Lct, LINE, Vals, RC )
2331 0 : IF ( RC /= HCO_SUCCESS ) THEN
2332 0 : CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
2333 0 : RETURN
2334 : ENDIF
2335 :
2336 : ! Check data / array dimensions
2337 0 : NT = SIZE(Vals,1)
2338 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
2339 0 : HcoState%NX, HcoState%NY, NT, RC )
2340 0 : IF ( RC /= HCO_SUCCESS ) THEN
2341 0 : CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
2342 0 : RETURN
2343 : ENDIF
2344 :
2345 : ! Pass to data array. If the country ID is larger than zero, fill
2346 : ! only those grid boxes. Otherwise, fill all grid boxes that have
2347 : ! not yet been filled.
2348 0 : DO I = 1, NT
2349 0 : IF ( CID == 0 ) THEN
2350 0 : WHERE ( Lct%Dct%Dta%V2(I)%Val <= 0.0_sp )
2351 0 : Lct%Dct%Dta%V2(I)%Val = Vals(I)
2352 : ENDWHERE
2353 : ELSE
2354 0 : WHERE ( CIDS == CID )
2355 0 : Lct%Dct%Dta%V2(I)%Val = Vals(I)
2356 : ENDWHERE
2357 : ENDIF
2358 : ENDDO
2359 :
2360 : ! Verbose
2361 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
2362 0 : WRITE(MSG,*) '- Obtained values for ',TRIM(CNT),' ==> ID:', CID
2363 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2364 : ENDIF
2365 :
2366 : ! Cleanup
2367 0 : IF ( ASSOCIATED(Vals) ) DEALLOCATE( Vals )
2368 0 : Vals => NULL()
2369 :
2370 : ! Update # of read lines
2371 0 : NLINE = NLINE + 1
2372 : ENDDO
2373 :
2374 : ! Close file
2375 0 : CLOSE ( IUFILE )
2376 :
2377 : ! Data is 2D
2378 0 : Lct%Dct%Dta%SpaceDim = 2
2379 :
2380 : ! Make sure data is in local time
2381 0 : IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
2382 0 : Lct%Dct%Dta%IsLocTime = .TRUE.
2383 : MSG = 'Data assigned to mask regions will be treated in local time: '//&
2384 0 : TRIM(Lct%Dct%cName)
2385 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2386 : ENDIF
2387 :
2388 : ! Cleanup
2389 0 : Cntr => NULL()
2390 0 : IF ( ALLOCATED(CIDS) ) DEALLOCATE ( CIDS )
2391 :
2392 : ! Return w/ success
2393 0 : RC = HCO_SUCCESS
2394 :
2395 0 : END SUBROUTINE HCOIO_ReadCountryValues
2396 : !EOC
2397 : !------------------------------------------------------------------------------
2398 : ! Harmonized Emissions Component (HEMCO) !
2399 : !------------------------------------------------------------------------------
2400 : !BOP
2401 : !
2402 : ! !IROUTINE: HCOIO_ReadFromConfig
2403 : !
2404 : ! !DESCRIPTION: Subroutine HCOIO\_ReadFromConfig reads data directly from
2405 : ! the configuration file (instead of reading it from a netCDF file).
2406 : ! These data is always assumed to be spatially uniform, but it is possible
2407 : ! to specify multiple time slices by separating the individual time slice
2408 : ! values by the HEMCO separator sign ('/' by default). The time dimension
2409 : ! of these data is either determined from the srcTime attribute or estimated
2410 : ! from the number of time slices provided. For example, if no srcTime is
2411 : ! specified and 24 time slices are provided, data is assumed to represent
2412 : ! hourly data. Similarly, data is assumed to represent weekdaily or monthly
2413 : ! data for 7 or 12 time slices, respectively.
2414 : !\\
2415 : !\\
2416 : ! If the srcTime attribute is defined, the time slices are determined from
2417 : ! this attribute. Only one time dimension (year, month, day, or hour) can
2418 : ! be defined for scalar fields!
2419 : !\\
2420 : !\\
2421 : ! !INTERFACE:
2422 : !
2423 0 : SUBROUTINE HCOIO_ReadFromConfig( HcoState, Lct, RC )
2424 : !
2425 : ! !USES:
2426 : !
2427 : USE HCO_FILEDATA_MOD, ONLY : FileData_ArrCheck
2428 : !
2429 : ! !INPUT PARAMTERS:
2430 : !
2431 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
2432 : !
2433 : ! !INPUT/OUTPUT PARAMETERS:
2434 : !
2435 : TYPE(ListCont), POINTER :: Lct
2436 : INTEGER, INTENT(INOUT) :: RC
2437 : !
2438 : ! !REVISION HISTORY:
2439 : ! 24 Jul 2014 - C. Keller: Initial version
2440 : ! See https://github.com/geoschem/hemco for complete history
2441 : !EOP
2442 : !------------------------------------------------------------------------------
2443 : !BOC
2444 : !
2445 : ! !LOCAL VARIABLES:
2446 : !
2447 : INTEGER :: I, NT
2448 0 : REAL(hp), POINTER :: Vals(:)
2449 : CHARACTER(LEN=255) :: MSG
2450 : CHARACTER(LEN=255) :: LOC = 'HCOIO_ReadFromConfig (hcoio_util_mod.F90)'
2451 :
2452 : !======================================================================
2453 : ! HCOIO_ReadFromConfig begins here
2454 : !======================================================================
2455 :
2456 : ! Init
2457 0 : Vals => NULL()
2458 :
2459 : ! Verbose
2460 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
2461 0 : WRITE(MSG, *) 'Read from config file: ', TRIM(Lct%Dct%cName)
2462 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2463 : ENDIF
2464 :
2465 : !-------------------------------------------------------------------
2466 : ! Get data values for this time step.
2467 : !-------------------------------------------------------------------
2468 0 : CALL GetDataVals( HcoState, Lct, Lct%Dct%Dta%ncFile, Vals, RC )
2469 0 : IF ( RC /= HCO_SUCCESS ) THEN
2470 0 : CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
2471 0 : RETURN
2472 : ENDIF
2473 :
2474 : !-------------------------------------------------------------------
2475 : ! Copy data into array.
2476 : !-------------------------------------------------------------------
2477 :
2478 : ! Number of values
2479 0 : NT = SIZE(Vals,1)
2480 :
2481 : ! For masks, interpret data as mask corners (lon1/lat1/lon2/lat2)
2482 : ! with no time dimension
2483 0 : IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
2484 :
2485 : ! Make sure data is allocated
2486 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
2487 0 : HcoState%NX, HcoState%NY, 1, RC )
2488 0 : IF ( RC /= HCO_SUCCESS ) THEN
2489 0 : CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
2490 0 : RETURN
2491 : ENDIF
2492 :
2493 : ! Fill array: 1.0 within grid box, 0.0 outside.
2494 0 : CALL FillMaskBox( HcoState, Lct, Vals, RC )
2495 0 : IF ( RC /= HCO_SUCCESS ) THEN
2496 0 : CALL HCO_ERROR( 'ERROR 16', RC, THISLOC=LOC )
2497 0 : RETURN
2498 : ENDIF
2499 :
2500 : ! Data is 2D
2501 0 : Lct%Dct%Dta%SpaceDim = 2
2502 :
2503 : ! For base emissions and scale factors, interpret data as scalar
2504 : ! values with a time dimension.
2505 : ELSE
2506 :
2507 0 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, 1, 1, NT, RC )
2508 0 : IF ( RC /= HCO_SUCCESS ) THEN
2509 0 : CALL HCO_ERROR( 'ERROR 17', RC, THISLOC=LOC )
2510 0 : RETURN
2511 : ENDIF
2512 0 : DO I = 1, NT
2513 0 : Lct%Dct%Dta%V2(I)%Val(1,1) = Vals(I)
2514 : !==============================================================================
2515 : ! KLUDGE BY BOB YANTOSCA (05 Jan 2016)
2516 : !
2517 : ! This WRITE statement avoids a seg fault in some Intel Fortran Compiler
2518 : ! versions, such as ifort 12 and ifort 13. The ADVANCE="no" prevents
2519 : ! carriage returns from being added to the log file, and the '' character
2520 : ! will prevent text from creeping across the screen.
2521 : !
2522 : ! NOTE: This section only gets executed during the initialization phase,
2523 : ! when we save data not read from netCDF files into the HEMCO data structure.
2524 : ! This type of data includes scale factors and mask data specified as vectors
2525 : ! in the HEMCO configuration file. Therefore, this section will only get
2526 : ! executed at startup, so the WRITE statment should not add significant
2527 : ! overhead to the simulation.
2528 : !
2529 : ! The root issue seems to be an optimization bug in the compiler.
2530 : !==============================================================================
2531 : #if defined( LINUX_IFORT )
2532 : WRITE( 6, '(a)', ADVANCE='no' ) ''
2533 : #endif
2534 :
2535 : ENDDO
2536 :
2537 : ! Data is 1D
2538 0 : Lct%Dct%Dta%SpaceDim = 1
2539 :
2540 : ! Make sure data is in local time
2541 0 : IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
2542 0 : Lct%Dct%Dta%IsLocTime = .TRUE.
2543 : MSG = 'Scale factors read from file are treated as local time: '// &
2544 0 : TRIM(Lct%Dct%cName)
2545 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2546 : ENDIF
2547 :
2548 : ENDIF
2549 :
2550 : ! Cleanup
2551 0 : IF ( ASSOCIATED(Vals) ) DEALLOCATE(Vals)
2552 :
2553 : ! Return w/ success
2554 0 : RC = HCO_SUCCESS
2555 :
2556 0 : END SUBROUTINE HCOIO_ReadFromConfig
2557 : !EOC
2558 : !------------------------------------------------------------------------------
2559 : ! Harmonized Emissions Component (HEMCO) !
2560 : !------------------------------------------------------------------------------
2561 : !BOP
2562 : !
2563 : ! !IROUTINE: GetSliceIdx
2564 : !
2565 : ! !DESCRIPTION: gets the time slice index to be used for data directly
2566 : ! read from the HEMCO configuration file. prefDt denotes the preferred
2567 : ! time attribute (year, month, or day). DtType is used to identify the
2568 : ! time attribute type (1=year, 2=month, 3=day). The time slice index will
2569 : ! be selected based upon those two variables. IDX is the selected time
2570 : ! slice index. It will be set to -1 if the current simulation date
2571 : ! is outside of the specified time range and the time cycle attribute is
2572 : ! not enabled for this field.
2573 : !\\
2574 : !\\
2575 : ! !INTERFACE:
2576 : !
2577 0 : SUBROUTINE GetSliceIdx ( HcoState, Lct, DtType, prefDt, IDX, RC )
2578 : !
2579 : ! !INPUT PARAMETERS:
2580 : !
2581 : TYPE(HCO_State), POINTER :: HcoState
2582 : TYPE(ListCont), POINTER :: Lct
2583 : INTEGER, INTENT(IN ) :: DtType
2584 : INTEGER, INTENT(IN ) :: prefDt
2585 : !
2586 : ! !INPUT/OUTPUT PARAMETERS:
2587 : !
2588 : INTEGER, INTENT(INOUT) :: IDX
2589 : INTEGER, INTENT(INOUT) :: RC
2590 : !
2591 : ! !REVISION HISTORY:
2592 : ! 13 Mar 2013 - C. Keller - Initial version
2593 : ! See https://github.com/geoschem/hemco for complete history
2594 : !EOP
2595 : !------------------------------------------------------------------------------
2596 : !BOC
2597 : !
2598 : ! !LOCAL VARIABLES:
2599 : !
2600 : INTEGER :: lowDt, uppDt
2601 : CHARACTER(LEN=255) :: MSG
2602 : CHARACTER(LEN=255) :: LOC = 'GetSliceIdx (hcoio_util_mod.F90)'
2603 :
2604 : !=================================================================
2605 : ! GetSliceIdx begins here!
2606 : !=================================================================
2607 :
2608 : ! Init
2609 0 : RC = HCO_SUCCESS
2610 :
2611 : ! Get upper and lower time range
2612 0 : IF ( DtType == 1 ) THEN
2613 0 : lowDt = Lct%Dct%Dta%ncYrs(1)
2614 0 : uppDt = Lct%Dct%Dta%ncYrs(2)
2615 0 : ELSEIF ( DtType == 2 ) THEN
2616 0 : lowDt = Lct%Dct%Dta%ncMts(1)
2617 0 : uppDt = Lct%Dct%Dta%ncMts(2)
2618 0 : ELSEIF ( DtType == 3 ) THEN
2619 0 : lowDt = Lct%Dct%Dta%ncDys(1)
2620 0 : uppDt = Lct%Dct%Dta%ncDys(2)
2621 : ELSE
2622 0 : WRITE(MSG,*) "DtType must be one of 1, 2, 3: ", DtType
2623 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2624 0 : RETURN
2625 : ENDIF
2626 :
2627 : ! Check for cycle flags:
2628 :
2629 : ! Data cycle set to range or exact date: in these cases, the
2630 : ! the preferred date will be equal to the current date, so
2631 : ! check if the preferred date is indeed within the available
2632 : ! range (lowDt, uppDt).
2633 : ! For data only to be used within the specified range, set
2634 : ! index to -1. This will force the scale factors to be set to
2635 : ! zero!
2636 0 : IF ( prefDt < lowDt .OR. prefDt > uppDt ) THEN
2637 0 : IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) .OR. &
2638 : ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) ) THEN
2639 0 : IDX = -1
2640 0 : RETURN
2641 : ELSE
2642 : ! this here should never happen, since for a cycle flag of 1,
2643 : ! the preferred date should always be restricted to the range
2644 : ! of available time stamps.
2645 0 : MSG = 'preferred date is outside of range: ' // TRIM(Lct%Dct%cName)
2646 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2647 0 : RETURN
2648 : ENDIF
2649 : ENDIF
2650 :
2651 : ! If the code makes it to here, prefDt is within the available data range
2652 : ! and we simply get the wanted index from the current index and the lowest
2653 : ! available index.
2654 0 : IDX = prefDt - lowDt + 1
2655 :
2656 : END SUBROUTINE GetSliceIdx
2657 : !EOC
2658 : !------------------------------------------------------------------------------
2659 : ! Harmonized Emissions Component (HEMCO) !
2660 : !------------------------------------------------------------------------------
2661 : !BOP
2662 : !
2663 : ! !IROUTINE: GetDataVals
2664 : !
2665 : ! !DESCRIPTION: Subroutine GetDataVals extracts the data values from ValStr
2666 : ! and writes them into vector Vals. ValStr is typically a character string
2667 : ! read from an external ASCII file or directly from the HEMCO configuration
2668 : ! file. Depending on the time specifications provided in the configuration
2669 : ! file, Vals will be filled with only a subset of the values of ValStr.
2670 : !\\
2671 : !\\
2672 : ! !INTERFACE:
2673 : !
2674 0 : SUBROUTINE GetDataVals ( HcoState, Lct, ValStr, Vals, RC )
2675 : !
2676 : ! !USES:
2677 : !
2678 : USE HCO_CHARTOOLS_MOD, ONLY : HCO_CharSplit
2679 : USE HCO_EXTLIST_MOD, ONLY : HCO_GetOpt
2680 : USE HCO_UNIT_MOD, ONLY : HCO_Unit_Change
2681 : USE HCO_tIdx_Mod, ONLY : HCO_GetPrefTimeAttr
2682 : USE HCO_CLOCK_MOD, ONLY : HcoClock_Get
2683 : !
2684 : ! !INPUT PARAMTERS:
2685 : !
2686 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
2687 : CHARACTER(LEN=*), INTENT(IN ) :: ValStr
2688 : !
2689 : ! !INPUT/OUTPUT PARAMETERS:
2690 : !
2691 : TYPE(ListCont), POINTER :: Lct
2692 : INTEGER, INTENT(INOUT) :: RC
2693 : !
2694 : ! !OUTPUT PARAMETERS:
2695 : !
2696 : REAL(hp), POINTER :: Vals(:)
2697 : !
2698 : ! !REVISION HISTORY:
2699 : ! 22 Dec 2014 - C. Keller: Initial version
2700 : ! See https://github.com/geoschem/hemco for complete history
2701 : !EOP
2702 : !------------------------------------------------------------------------------
2703 : !BOC
2704 : !
2705 : ! !LOCAL VARIABLES:
2706 : !
2707 : INTEGER :: HcoID
2708 : INTEGER :: I, N, NUSE, AS
2709 : INTEGER :: IDX1, IDX2
2710 : INTEGER :: AreaFlag, TimeFlag, Check
2711 : INTEGER :: prefYr, prefMt, prefDy, prefHr, prefMn
2712 : INTEGER :: cYr, cMt, cDy, cHr
2713 : REAL(hp) :: MW_g
2714 : REAL(hp) :: UnitFactor
2715 : REAL(hp) :: FileVals(100)
2716 0 : REAL(hp), POINTER :: FileArr(:,:,:,:)
2717 : LOGICAL :: IsPerArea
2718 : LOGICAL :: IsMath
2719 : CHARACTER(LEN=255) :: MSG
2720 : CHARACTER(LEN=255) :: LOC = 'GetDataVals (hcoio_util_mod.F90)'
2721 :
2722 : !======================================================================
2723 : ! GetDataVals begins here
2724 : !======================================================================
2725 :
2726 : ! Initialize
2727 0 : FileArr => NULL()
2728 :
2729 : ! Shadow species properties needed for unit conversion
2730 0 : HcoID = Lct%Dct%HcoID
2731 0 : IF ( HcoID > 0 ) THEN
2732 0 : MW_g = HcoState%Spc(HcoID)%MW_g
2733 : ELSE
2734 0 : MW_g = -999.0_hp
2735 : ENDIF
2736 :
2737 : ! Is this a math expression?
2738 0 : IsMath = .FALSE.
2739 0 : IF ( LEN(ValStr) > 5 ) THEN
2740 0 : IF ( ValStr(1:5)=='MATH:' ) IsMath = .TRUE.
2741 : ENDIF
2742 :
2743 : ! Evaluate math expression if string starts with 'MATH:'
2744 : IF ( IsMath ) THEN
2745 0 : CALL ReadMath ( HcoState, Lct, ValStr, FileVals, N, RC )
2746 0 : IF ( RC /= HCO_SUCCESS ) THEN
2747 0 : CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
2748 0 : RETURN
2749 : ENDIF
2750 :
2751 : ! Use regular string parser otherwise
2752 : ELSE
2753 : CALL HCO_CharSplit ( ValStr, &
2754 : HCO_GetOpt(HcoState%Config%ExtList,'Separator'), &
2755 : HCO_GetOpt(HcoState%Config%ExtList,'Wildcard'), &
2756 0 : FileVals, N, RC )
2757 0 : IF ( RC /= HCO_SUCCESS ) THEN
2758 0 : CALL HCO_ERROR( 'ERROR 19', RC, THISLOC=LOC )
2759 0 : RETURN
2760 : ENDIF
2761 : ENDIF
2762 :
2763 : ! Return w/ error if no scale factor defined
2764 0 : IF ( N == 0 ) THEN
2765 : MSG = 'Cannot read data: ' // TRIM(Lct%Dct%cName) // &
2766 0 : ': ' // TRIM(ValStr)
2767 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC)
2768 0 : RETURN
2769 : ENDIF
2770 :
2771 : ! Get the preferred times, i.e. the preferred year, month, day,
2772 : ! or hour (as specified in the configuration file).
2773 : CALL HCO_GetPrefTimeAttr( HcoState, Lct, &
2774 0 : prefYr, prefMt, prefDy, prefHr, prefMn, RC )
2775 0 : IF ( RC /= HCO_SUCCESS ) THEN
2776 0 : CALL HCO_ERROR( 'ERROR 20', RC, THISLOC=LOC )
2777 0 : RETURN
2778 : ENDIF
2779 :
2780 : ! ----------------------------------------------------------------
2781 : ! For masks, assume that values represent the corners of the mask
2782 : ! box, e.g. there must be four values. Masks are time-independent
2783 : ! and unitless
2784 : ! ----------------------------------------------------------------
2785 0 : IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
2786 :
2787 : ! There must be exactly four values
2788 0 : IF ( N /= 4 ) THEN
2789 : MSG = 'Mask values are not lon1/lat1/lon2/lat2: ' // &
2790 0 : TRIM(ValStr) // ' --> ' // TRIM(Lct%Dct%cName)
2791 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2792 0 : RETURN
2793 : ENDIF
2794 :
2795 : ! Pass to FileArr array (will be used below)
2796 0 : NUSE = 4
2797 0 : ALLOCATE( FileArr(1,1,1,NUSE), STAT=AS )
2798 : IF ( AS /= 0 ) THEN
2799 0 : MSG = 'Cannot allocate FileArr'
2800 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2801 0 : RETURN
2802 : ENDIF
2803 0 : FileArr(1,1,1,:) = FileVals(1:NUSE)
2804 :
2805 : ! ----------------------------------------------------------------
2806 : ! For non-masks, the data is interpreted as uniform values with
2807 : ! a time dimension. Need to select the time slices to be used at
2808 : ! this time (depending on the provided time attributes), as well
2809 : ! as to ensure that values are in the correct units.
2810 : ! Use all time slices unless a time interval is provided in
2811 : ! attribute srcTime of the configuration file.
2812 : ! ----------------------------------------------------------------
2813 : ELSE
2814 :
2815 : ! If there is only one value use this one and ignore any time
2816 : ! preferences.
2817 0 : IF ( N == 1 ) THEN
2818 0 : NUSE = 1
2819 0 : IDX1 = 1
2820 0 : IDX2 = 1
2821 :
2822 : ! If it's a math expression use all passed values
2823 0 : ELSEIF ( IsMath ) THEN
2824 0 : NUSE = N
2825 0 : IDX1 = 1
2826 0 : IDX2 = N
2827 :
2828 : ELSE
2829 : ! Currently, data read directly from the configuration file can only
2830 : ! represent one time dimension, i.e. it can only be yearly, monthly,
2831 : ! daily (or hourly data, but this is read all at the same time).
2832 :
2833 : ! Annual data
2834 0 : IF ( Lct%Dct%Dta%ncYrs(1) /= Lct%Dct%Dta%ncYrs(2) ) THEN
2835 : ! Error check
2836 : IF ( Lct%Dct%Dta%ncMts(1) /= Lct%Dct%Dta%ncMts(2) .OR. &
2837 0 : Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) .OR. &
2838 : Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2) ) THEN
2839 : MSG = 'Data must not have more than one time dimension: ' // &
2840 0 : TRIM(Lct%Dct%cName)
2841 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2842 0 : RETURN
2843 : ENDIF
2844 :
2845 0 : CALL GetSliceIdx ( HcoState, Lct, 1, prefYr, IDX1, RC )
2846 0 : IF ( RC /= HCO_SUCCESS ) THEN
2847 0 : CALL HCO_ERROR( 'ERROR 21', RC, THISLOC=LOC )
2848 0 : RETURN
2849 : ENDIF
2850 0 : IDX2 = IDX1
2851 0 : NUSE = 1
2852 :
2853 : ! Monthly data
2854 0 : ELSEIF ( Lct%Dct%Dta%ncMts(1) /= Lct%Dct%Dta%ncMts(2) ) THEN
2855 : ! Error check
2856 0 : IF ( Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) .OR. &
2857 : Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2) ) THEN
2858 : MSG = 'Data must only have one time dimension: ' // &
2859 0 : TRIM(Lct%Dct%cName)
2860 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2861 0 : RETURN
2862 : ENDIF
2863 :
2864 0 : CALL GetSliceIdx ( HcoState, Lct, 2, prefMt, IDX1, RC )
2865 0 : IF ( RC /= HCO_SUCCESS ) THEN
2866 0 : CALL HCO_ERROR( 'ERROR 22', RC, THISLOC=LOC )
2867 0 : RETURN
2868 : ENDIF
2869 0 : IDX2 = IDX1
2870 0 : NUSE = 1
2871 :
2872 : ! Daily data
2873 0 : ELSEIF ( Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) ) THEN
2874 : ! Error check
2875 0 : IF ( Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2) ) THEN
2876 : MSG = 'Data must only have one time dimension: ' // &
2877 0 : TRIM(Lct%Dct%cName)
2878 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2879 0 : RETURN
2880 : ENDIF
2881 :
2882 0 : CALL GetSliceIdx ( HcoState, Lct, 3, prefDy, IDX1, RC )
2883 0 : IF ( RC /= HCO_SUCCESS ) THEN
2884 0 : CALL HCO_ERROR( 'ERROR 23', RC, THISLOC=LOC )
2885 0 : RETURN
2886 : ENDIF
2887 0 : IDX2 = IDX1
2888 0 : NUSE = 1
2889 :
2890 : ! All other cases (incl. hourly data): read all time slices).
2891 : ELSE
2892 0 : IDX1 = 1
2893 0 : IDX2 = N
2894 0 : NUSE = N
2895 : ENDIF
2896 : ENDIF
2897 :
2898 : ! ----------------------------------------------------------------
2899 : ! Read selected time slice(s) into data array
2900 : ! ----------------------------------------------------------------
2901 0 : IF ( IDX2 > N ) THEN
2902 0 : WRITE(MSG,*) 'Index ', IDX2, ' is larger than number of ', &
2903 0 : 'values found: ', TRIM(Lct%Dct%cName)
2904 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2905 0 : RETURN
2906 : ENDIF
2907 :
2908 0 : ALLOCATE( FileArr(1,1,1,NUSE), STAT=AS )
2909 : IF ( AS /= 0 ) THEN
2910 0 : MSG = 'Cannot allocate FileArr'
2911 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
2912 0 : RETURN
2913 : ENDIF
2914 :
2915 : ! Check for range/exact flag
2916 : ! If range is given, the preferred Yr/Mt/Dy/Hr will be negative
2917 : ! if we are outside the desired range.
2918 0 : IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) THEN
2919 0 : IF ( prefYr == -1 .OR. prefMt == -1 .OR. prefDy == -1 ) IDX1 = -1
2920 0 : IF ( Lct%Dct%Dta%ncHrs(1) >= 0 .AND. prefHr == -1 ) IDX1 = -1
2921 :
2922 : ! If flag is exact, the preferred date must be equal to the current
2923 : ! simulation date.
2924 0 : ELSEIF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) THEN
2925 0 : IF ( Lct%Dct%Dta%ncYrs(1) > 0 ) THEN
2926 0 : IF ( prefYr < Lct%Dct%Dta%ncYrs(1) .OR. &
2927 0 : prefYr > Lct%Dct%Dta%ncYrs(2) ) IDX1 = -1
2928 : ENDIF
2929 0 : IF ( Lct%Dct%Dta%ncMts(1) > 0 ) THEN
2930 0 : IF ( prefMt < Lct%Dct%Dta%ncMts(1) .OR. &
2931 0 : prefMt > Lct%Dct%Dta%ncMts(2) ) IDX1 = -1
2932 : ENDIF
2933 0 : IF ( Lct%Dct%Dta%ncDys(1) > 0 ) THEN
2934 0 : IF ( prefDy < Lct%Dct%Dta%ncDys(1) .OR. &
2935 0 : prefDy > Lct%Dct%Dta%ncDys(2) ) IDX1 = -1
2936 : ENDIF
2937 0 : IF ( Lct%Dct%Dta%ncHrs(1) >= 0 ) THEN
2938 0 : IF ( prefHr < Lct%Dct%Dta%ncHrs(1) .OR. &
2939 0 : prefHr > Lct%Dct%Dta%ncHrs(2) ) IDX1 = -1
2940 : ENDIF
2941 : ENDIF
2942 :
2943 : ! IDX1 becomes -1 for data that is outside of the valid range
2944 : ! (and no time cycling enabled). In this case, make sure that
2945 : ! scale factor is set to zero.
2946 0 : IF ( IDX1 < 0 ) THEN
2947 0 : IF ( Lct%Dct%DctType == HCO_DCTTYPE_BASE ) THEN
2948 0 : FileArr(1,1,1,:) = 0.0_hp
2949 : MSG = 'Base field outside of range - set to zero: ' // &
2950 0 : TRIM(Lct%Dct%cName)
2951 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2952 : #if defined( MODEL_GEOS )
2953 : ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
2954 : FileArr(1,1,1,:) = 0.0_hp
2955 : MSG = 'Mask outside of range - set to zero: ' // &
2956 : TRIM(Lct%Dct%cName)
2957 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2958 : #endif
2959 : ELSE
2960 0 : FileArr(1,1,1,:) = 1.0_hp
2961 : MSG = 'Scale factor outside of range - set to one: ' // &
2962 0 : TRIM(Lct%Dct%cName)
2963 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2964 : ENDIF
2965 : ELSE
2966 0 : FileArr(1,1,1,:) = FileVals(IDX1:IDX2)
2967 : ENDIF
2968 :
2969 : ! ----------------------------------------------------------------
2970 : ! Convert data to HEMCO units
2971 : ! ----------------------------------------------------------------
2972 : CALL HCO_UNIT_CHANGE( HcoConfig = HcoState%Config, &
2973 : Array = FileArr, &
2974 : Units = TRIM(Lct%Dct%Dta%OrigUnit), &
2975 : MW = MW_g, &
2976 : YYYY = -999, &
2977 : MM = -999, &
2978 : AreaFlag = AreaFlag, &
2979 : TimeFlag = TimeFlag, &
2980 : FACT = UnitFactor, &
2981 0 : RC = RC )
2982 0 : IF ( RC /= HCO_SUCCESS ) THEN
2983 0 : CALL HCO_ERROR( 'ERROR 24', RC, THISLOC=LOC )
2984 0 : RETURN
2985 : ENDIF
2986 :
2987 : ! Verbose mode
2988 0 : IF ( UnitFactor /= 1.0_hp ) THEN
2989 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
2990 0 : WRITE(MSG,*) 'Data was in units of ', TRIM(Lct%Dct%Dta%OrigUnit), &
2991 0 : ' - converted to HEMCO units by applying ', &
2992 0 : 'scale factor ', UnitFactor
2993 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2994 : ENDIF
2995 : ELSE
2996 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
2997 0 : WRITE(MSG,*) 'Data was in units of ', TRIM(Lct%Dct%Dta%OrigUnit), &
2998 0 : ' - unit conversion factor is ', UnitFactor
2999 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3000 : ENDIF
3001 : ENDIF
3002 :
3003 : ! Data must be ...
3004 : ! ... concentration ...
3005 0 : IF ( AreaFlag == 3 .AND. TimeFlag == 0 ) THEN
3006 0 : Lct%Dct%Dta%IsConc = .TRUE.
3007 :
3008 0 : ELSEIF ( AreaFlag == 3 .AND. TimeFlag == 1 ) THEN
3009 0 : Lct%Dct%Dta%IsConc = .TRUE.
3010 0 : FileArr = FileArr * HcoState%TS_EMIS
3011 : MSG = 'Data converted from kg/m3/s to kg/m3: ' // &
3012 0 : TRIM(Lct%Dct%cName) // ': ' // TRIM(Lct%Dct%Dta%OrigUnit)
3013 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
3014 :
3015 : ! ... emissions or unitless ...
3016 0 : ELSEIF ( (AreaFlag == -1 .AND. TimeFlag == -1) .OR. &
3017 : (AreaFlag == 2 .AND. TimeFlag == 1) ) THEN
3018 0 : Lct%Dct%Dta%IsConc = .FALSE.
3019 :
3020 : ! ... invalid otherwise:
3021 : ELSE
3022 : MSG = 'Unit must be unitless, emission or concentration: ' // &
3023 0 : TRIM(Lct%Dct%cName) // ': ' // TRIM(Lct%Dct%Dta%OrigUnit)
3024 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3025 0 : RETURN
3026 : ENDIF
3027 :
3028 : ! Auto-detect delta t [in hours] between time slices.
3029 : ! Scale factors can be:
3030 : ! length 1 : constant
3031 : ! length 7 : weekday factors: Sun, Mon, ..., Sat
3032 : ! length 12: monthly factors: Jan, Feb, ..., Dec
3033 : ! length 24: hourly factors: 12am, 1am, ... 11pm
3034 0 : IF ( NUSE == 1 ) THEN
3035 0 : Lct%Dct%Dta%DeltaT = 0
3036 0 : ELSEIF ( NUSE == 7 ) THEN
3037 0 : Lct%Dct%Dta%DeltaT = 24
3038 0 : ELSEIF ( NUSE == 12 ) THEN
3039 0 : Lct%Dct%Dta%DeltaT = 720
3040 0 : ELSEIF ( NUSE == 24 ) THEN
3041 0 : Lct%Dct%Dta%DeltaT = 1
3042 : ELSE
3043 : MSG = 'Factor must be of length 1, 7, 12, or 24!' // &
3044 0 : TRIM(Lct%Dct%cName)
3045 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC)
3046 0 : RETURN
3047 : ENDIF
3048 :
3049 : ENDIF ! Masks vs. non-masks
3050 :
3051 : ! Copy data into output array.
3052 0 : IF ( ASSOCIATED(Vals) ) DEALLOCATE( Vals )
3053 0 : ALLOCATE( Vals(NUSE), STAT=AS )
3054 : IF ( AS /= 0 ) THEN
3055 0 : MSG = 'Cannot allocate Vals'
3056 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3057 0 : RETURN
3058 : ENDIF
3059 0 : Vals(:) = FileArr(1,1,1,:)
3060 :
3061 : ! Cleanup
3062 0 : IF ( ASSOCIATED(FileArr) ) DEALLOCATE(FileArr)
3063 0 : FileArr => NULL()
3064 :
3065 : ! Return w/ success
3066 0 : RC = HCO_SUCCESS
3067 :
3068 0 : END SUBROUTINE GetDataVals
3069 : !EOC
3070 : !------------------------------------------------------------------------------
3071 : ! Harmonized Emissions Component (HEMCO) !
3072 : !------------------------------------------------------------------------------
3073 : !BOP
3074 : !
3075 : ! !IROUTINE: FillMaskBox
3076 : !
3077 : ! !DESCRIPTION: Subroutine FillMaskBox fills the data array of the passed list
3078 : ! container Lct according to the mask region provided in Vals. Vals contains
3079 : ! the mask region of interest, denoted by the lower left and upper right grid
3080 : ! box corners: lon1, lat1, lon2, lat2. The data array of Lct is filled such
3081 : ! that all grid boxes are set to 1 whose mid-point is inside of the given box
3082 : ! range.
3083 : !\\
3084 : !\\
3085 : ! !INTERFACE:
3086 : !
3087 0 : SUBROUTINE FillMaskBox ( HcoState, Lct, Vals, RC )
3088 : !
3089 : ! !USES:
3090 : !
3091 : !
3092 : ! !INPUT PARAMTERS:
3093 : !
3094 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
3095 : REAL(hp) , POINTER :: Vals(:)
3096 : !
3097 : ! !INPUT/OUTPUT PARAMETERS:
3098 : !
3099 : TYPE(ListCont), POINTER :: Lct
3100 : INTEGER, INTENT(INOUT) :: RC
3101 : !
3102 : ! !REVISION HISTORY:
3103 : ! 29 Dec 2014 - C. Keller - Initial version
3104 : ! See https://github.com/geoschem/hemco for complete history
3105 : !EOP
3106 : !------------------------------------------------------------------------------
3107 : !BOC
3108 : !
3109 : ! !LOCAL VARIABLES:
3110 : !
3111 : LOGICAL :: GridPoint
3112 : INTEGER :: I, J
3113 : REAL(hp) :: LON1, LON2, LAT1, LAT2
3114 : REAL(hp) :: XDG1, XDG2, YDG1, YDG2
3115 : REAL(hp) :: ILON, ILAT
3116 : CHARACTER(LEN=255) :: MSG
3117 : CHARACTER(LEN=255) :: LOC = 'FillMaskBox (hcoio_util_mod.F90)'
3118 :
3119 : !=================================================================
3120 : ! FillMaskBox begins here!
3121 : !=================================================================
3122 :
3123 : ! Extract lon1, lon2, lat1, lat2
3124 0 : LON1 = VALS(1)
3125 0 : LAT1 = VALS(2)
3126 0 : LON2 = VALS(3)
3127 0 : LAT2 = VALS(4)
3128 :
3129 : ! Check if this is mask is a point. In this case, we need the grid
3130 : ! box edges being defined.
3131 0 : GridPoint = .FALSE.
3132 0 : IF ( ( LON1 == LON2 ) .AND. ( LAT1 == LAT2 ) ) THEN
3133 0 : IF ( .NOT. ASSOCIATED(HcoState%Grid%XEDGE%Val) .OR. &
3134 : .NOT. ASSOCIATED(HcoState%Grid%YEDGE%Val) ) THEN
3135 : MSG = 'Cannot evaluate grid point mask - need grid box ' // &
3136 : 'edges for this. This error occurs if a mask covers '// &
3137 : 'a fixed grid point (e.g. lon1=lon2 and lat1=lat2) ' // &
3138 0 : 'but HEMCO grid edges are not defined.'
3139 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3140 0 : RETURN
3141 : ENDIF
3142 : GridPoint = .TRUE.
3143 : ENDIF
3144 :
3145 : ! Check for every grid box if mid point is within mask region.
3146 : ! Set to 1.0 if this is the case.
3147 : !$OMP PARALLEL DO &
3148 : !$OMP DEFAULT( SHARED ) &
3149 : !$OMP PRIVATE( I, J, ILON, ILAT ) &
3150 : !$OMP PRIVATE( XDG1, XDG2, YDG1, YDG2 ) &
3151 : !$OMP SCHEDULE( DYNAMIC )
3152 0 : DO J = 1, HcoState%NY
3153 0 : DO I = 1, HcoState%NX
3154 :
3155 : ! If it's a grid point, check if it's within this
3156 : ! grid box
3157 0 : IF ( GridPoint ) THEN
3158 0 : XDG1 = HcoState%Grid%XEDGE%Val(I ,J )
3159 0 : XDG2 = HcoState%Grid%XEDGE%Val(I+1,J )
3160 0 : YDG1 = HcoState%Grid%YEDGE%Val(I ,J )
3161 0 : YDG2 = HcoState%Grid%YEDGE%Val(I ,J+1)
3162 0 : IF ( XDG1 >= 180.0_hp ) XDG1 = XDG1 - 360.0_hp
3163 0 : IF ( XDG2 >= 180.0_hp ) XDG2 = XDG2 - 360.0_hp
3164 :
3165 : IF ( LON1 >= XDG1 .AND. LON1 <= XDG2 .AND. &
3166 0 : LAT1 >= YDG1 .AND. LAT1 <= YDG2 ) THEN
3167 0 : Lct%Dct%Dta%V2(1)%Val(I,J) = 1.0_sp
3168 : ENDIF
3169 :
3170 : ! Check if mid point is within mask region
3171 : ELSE
3172 : ! Get longitude and latitude at this grid box
3173 0 : ILON = HcoState%Grid%XMID%Val(I,J)
3174 0 : ILAT = HcoState%Grid%YMID%Val(I,J)
3175 0 : IF ( ILON >= 180.0_hp ) ILON = ILON - 360.0_hp
3176 :
3177 : IF ( ILON >= LON1 .AND. ILON <= LON2 .AND. &
3178 0 : ILAT >= LAT1 .AND. ILAT <= LAT2 ) THEN
3179 0 : Lct%Dct%Dta%V2(1)%Val(I,J) = 1.0_sp
3180 : ENDIF
3181 : ENDIF
3182 :
3183 : ENDDO
3184 : ENDDO
3185 : !$OMP END PARALLEL DO
3186 :
3187 : ! Return w/ success
3188 0 : RC = HCO_SUCCESS
3189 :
3190 : END SUBROUTINE FillMaskBox
3191 : !EOC
3192 : !------------------------------------------------------------------------------
3193 : ! Harmonized Emissions Component (HEMCO) !
3194 : !------------------------------------------------------------------------------
3195 : !BOP
3196 : !
3197 : ! !IROUTINE: ReadMath
3198 : !
3199 : ! !DESCRIPTION: Subroutine ReadMath reads and evaluates a mathematical
3200 : ! expression. Mathematical expressions can combine time-stamps with
3201 : ! mathematical functions, e.g. to yield the sine of current simulation hour.
3202 : ! Mathematical expressions must start with the identifier 'MATH:', followed
3203 : ! by the actual expression. Each expression must include at least one
3204 : ! variable (evaluated at runtime). The following variables are currently
3205 : ! supported: YYYY (year), MM (month), DD (day), HH (hour), LH (local hour),
3206 : ! NN (minute), SS (second), WD (weekday), LWD (local weekday),
3207 : ! DOY (day of year), ELH (elapsed hours), ELS (elapsed seconds).
3208 : ! In addition, the following variables can be used: PI (3.141...), DOM
3209 : ! (\# of days of current month).
3210 : ! For example, the following expression would yield a continuous sine
3211 : ! curve as function of hour of day: 'MATH:sin(HH/24*PI*2)'.
3212 : !\\
3213 : !\\
3214 : ! For a full list of valid mathematical expressions, see module interpreter.F90.
3215 : !\\
3216 : !\\
3217 : ! !INTERFACE:
3218 : !
3219 0 : SUBROUTINE ReadMath( HcoState, Lct, ValStr, Vals, N, RC )
3220 : !
3221 : ! !USES:
3222 : !
3223 : USE HCO_CLOCK_MOD, ONLY : HcoClock_Get
3224 : USE HCO_tIdx_Mod, ONLY : HCO_GetPrefTimeAttr
3225 : USE INTERPRETER
3226 : !
3227 : ! !INPUT PARAMTERS:
3228 : !
3229 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
3230 : TYPE(ListCont), POINTER :: Lct
3231 : CHARACTER(LEN=*), INTENT(IN ) :: ValStr
3232 : !
3233 : ! !INPUT/OUTPUT PARAMETERS:
3234 : !
3235 : REAL(hp), INTENT(INOUT) :: Vals(:)
3236 : INTEGER, INTENT(INOUT) :: RC
3237 : !
3238 : ! !OUTPUT PARAMETERS:
3239 : !
3240 : INTEGER, INTENT( OUT) :: N
3241 : !
3242 : ! !REVISION HISTORY:
3243 : ! 11 May 2017 - C. Keller - Initial version
3244 : ! See https://github.com/geoschem/hemco for complete history
3245 : !EOP
3246 : !------------------------------------------------------------------------------
3247 : !BOC
3248 : !
3249 : ! !LOCAL VARIABLES:
3250 : !
3251 : LOGICAL :: EOS
3252 : INTEGER :: STRL
3253 : INTEGER :: I, NVAL, LHIDX, LWDIDX
3254 : INTEGER :: prefYr, prefMt, prefDy, prefHr, prefMn
3255 : INTEGER :: prefWD, prefDOY, prefS, LMD, cHr
3256 : INTEGER :: nSteps
3257 : REAL(hp) :: ELH, ELS
3258 : REAL(hp) :: Val
3259 : CHARACTER(LEN=255) :: MSG
3260 : CHARACTER(LEN=255) :: LOC = 'ReadMath (hcoio_util_mod.F90)'
3261 :
3262 : ! Variables used by the evaluator to build and to determine the value
3263 : ! of the expressions
3264 : character(len = 10) :: all_variables(12)
3265 : real(hp) :: all_variablesvalues(12)
3266 :
3267 : !String variable that will store the function that the evaluator will build
3268 : character (len = 275) :: func
3269 :
3270 : !String variable that will return the building of the expression result
3271 : !If everything was ok then statusflag = 'ok', otherwise statusflag = 'error'
3272 : character (len = 5) :: statusflag
3273 :
3274 : !======================================================================
3275 : ! ReadMath begins here
3276 : !======================================================================
3277 :
3278 : ! Substring (without flag 'MATH:')
3279 0 : STRL = LEN(ValStr)
3280 0 : IF ( STRL < 6 ) THEN
3281 : MSG = 'Math expression is too short - expected `MATH:<expr>`: ' &
3282 0 : //TRIM(ValStr)
3283 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3284 0 : RETURN
3285 : ENDIF
3286 0 : func = ValStr(6:STRL)
3287 :
3288 : ! Get preferred time stamps
3289 : CALL HCO_GetPrefTimeAttr( HcoState, Lct, &
3290 0 : prefYr, prefMt, prefDy, prefHr, prefMn, RC )
3291 0 : IF ( RC /= HCO_SUCCESS ) THEN
3292 0 : CALL HCO_ERROR( 'ERROR 25', RC, THISLOC=LOC )
3293 0 : RETURN
3294 : ENDIF
3295 :
3296 : ! Get some other current time stamps
3297 : CALL HcoClock_Get( HcoState%Clock, cS=prefS, cH=cHr, &
3298 : cWEEKDAY=prefWD, cDOY=prefDOY, LMD=LMD, &
3299 0 : nSteps=nSteps, RC=RC )
3300 0 : IF ( RC /= HCO_SUCCESS ) THEN
3301 0 : CALL HCO_ERROR( 'ERROR 26', RC, THISLOC=LOC )
3302 0 : RETURN
3303 : ENDIF
3304 :
3305 : ! GetPrefTimeAttr can return -999 for hour. In this case set to current
3306 : ! simulation hour
3307 0 : IF ( prefHr < 0 ) prefHr = cHr
3308 :
3309 : ! Parse function. This will replace any tokens in the function with the
3310 : ! actual token values. (ckeller, 7/7/17)
3311 : CALL HCO_CharParse ( HcoState%Config, func, &
3312 0 : prefYr, prefMt, prefDy, prefHr, prefMn, RC )
3313 0 : IF ( RC /= HCO_SUCCESS ) THEN
3314 0 : CALL HCO_ERROR( 'ERROR 27', RC, THISLOC=LOC )
3315 0 : RETURN
3316 : ENDIF
3317 :
3318 : ! Elapsed hours and seconds since start time
3319 0 : ELS = HcoState%TS_DYN * nSteps
3320 0 : ELH = ELS / 3600.0_hp
3321 :
3322 : ! Check which variables are in string.
3323 : ! Possible variables are YYYY, MM, DD, WD, HH, NN, SS, DOY, ELH, ELS
3324 0 : NVAL = 0
3325 0 : LHIDX = -1
3326 0 : LWDIDX = -1
3327 :
3328 0 : IF ( INDEX(func,'YYYY') > 0 ) THEN
3329 0 : NVAL = NVAL + 1
3330 0 : all_variables(NVAL) = 'yyyy'
3331 0 : all_variablesvalues(NVAL) = prefYr
3332 : ENDIF
3333 0 : IF ( INDEX(func,'MM') > 0 ) THEN
3334 0 : NVAL = NVAL + 1
3335 0 : all_variables(NVAL) = 'mm'
3336 0 : all_variablesvalues(NVAL) = prefMt
3337 : ENDIF
3338 0 : IF ( INDEX(func,'DD') > 0 ) THEN
3339 0 : NVAL = NVAL + 1
3340 0 : all_variables(NVAL) = 'dd'
3341 0 : all_variablesvalues(NVAL) = prefDy
3342 : ENDIF
3343 0 : IF ( INDEX(func,'WD') > 0 ) THEN
3344 0 : NVAL = NVAL + 1
3345 0 : all_variables(NVAL) = 'wd'
3346 0 : all_variablesvalues(NVAL) = prefWD
3347 : ENDIF
3348 0 : IF ( INDEX(func,'LWD') > 0 ) THEN
3349 0 : NVAL = NVAL + 1
3350 0 : all_variables(NVAL) = 'lwd'
3351 0 : all_variablesvalues(NVAL) = prefWD
3352 0 : LWDIDX = NVAL
3353 : ENDIF
3354 0 : IF ( INDEX(func,'HH') > 0 ) THEN
3355 0 : NVAL = NVAL + 1
3356 0 : all_variables(NVAL) = 'hh'
3357 0 : all_variablesvalues(NVAL) = prefHr
3358 : ENDIF
3359 0 : IF ( INDEX(func,'LH') > 0 ) THEN
3360 0 : NVAL = NVAL + 1
3361 0 : all_variables(NVAL) = 'lh'
3362 0 : all_variablesvalues(NVAL) = prefHr
3363 0 : LHIDX = NVAL
3364 : ENDIF
3365 0 : IF ( INDEX(func,'NN') > 0 ) THEN
3366 0 : NVAL = NVAL + 1
3367 0 : all_variables(NVAL) = 'nn'
3368 0 : all_variablesvalues(NVAL) = prefMn
3369 : ENDIF
3370 0 : IF ( INDEX(func,'SS') > 0 ) THEN
3371 0 : NVAL = NVAL + 1
3372 0 : all_variables(NVAL) = 'ss'
3373 0 : all_variablesvalues(NVAL) = prefS
3374 : ENDIF
3375 0 : IF ( INDEX(func,'DOY') > 0 ) THEN
3376 0 : NVAL = NVAL + 1
3377 0 : all_variables(NVAL) = 'doy'
3378 0 : all_variablesvalues(NVAL) = prefDOY
3379 : ENDIF
3380 0 : IF ( INDEX(func,'PI') > 0 ) THEN
3381 0 : NVAL = NVAL + 1
3382 0 : all_variables(NVAL) = 'pi'
3383 0 : all_variablesvalues(NVAL) = HcoState%Phys%PI
3384 : ENDIF
3385 0 : IF ( INDEX(func,'DOM') > 0 ) THEN
3386 0 : NVAL = NVAL + 1
3387 0 : all_variables(NVAL) = 'dom'
3388 0 : all_variablesvalues(NVAL) = LMD
3389 : ENDIF
3390 0 : IF ( INDEX(func,'ELH') > 0 ) THEN
3391 0 : NVAL = NVAL + 1
3392 0 : all_variables(NVAL) = 'elh'
3393 0 : all_variablesvalues(NVAL) = ELH
3394 : ENDIF
3395 0 : IF ( INDEX(func,'ELS') > 0 ) THEN
3396 0 : NVAL = NVAL + 1
3397 0 : all_variables(NVAL) = 'els'
3398 0 : all_variablesvalues(NVAL) = ELS
3399 : ENDIF
3400 :
3401 : ! Error trap: cannot have local hour and local weekday in
3402 : ! same expression
3403 0 : IF ( LHIDX > 0 .AND. LWDIDX > 0 ) THEN
3404 : MSG = 'Cannot have local hour and local weekday in '//&
3405 0 : 'same expression: '//TRIM(func)
3406 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3407 0 : RETURN
3408 : ENDIF
3409 :
3410 : ! N is the number of expressions.
3411 0 : Vals(:) = -999.0_hp
3412 0 : IF ( LHIDX > 0 ) THEN
3413 0 : N = 24
3414 0 : ELSEIF ( LWDIDX > 0 ) THEN
3415 0 : N = 7
3416 : ELSE
3417 0 : N = 1
3418 : ENDIF
3419 :
3420 : ! Evaluate expression
3421 : !Initialize function
3422 0 : call init (func, all_variables(1:NVAL), statusflag)
3423 0 : IF(statusflag == 'ok') THEN
3424 0 : DO I=1,N
3425 0 : IF ( LHIDX > 0 ) all_variablesvalues(LHIDX) = I-1
3426 0 : IF ( LWDIDX > 0 ) all_variablesvalues(LWDIDX) = I-1
3427 0 : Val = evaluate( all_variablesvalues(1:NVAL) )
3428 0 : Vals(I) = Val
3429 :
3430 : ! Verbose
3431 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
3432 0 : WRITE(MSG,*) 'Evaluated function: ',TRIM(func),' --> ', Val
3433 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3434 : ENDIF
3435 : ENDDO
3436 : ELSE
3437 0 : MSG = 'Error evaluation function: '//TRIM(func)
3438 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
3439 0 : RETURN
3440 : ENDIF
3441 0 : call destroyfunc()
3442 :
3443 : ! Return w/ success
3444 0 : RC = HCO_SUCCESS
3445 :
3446 : END SUBROUTINE ReadMath
3447 : !EOC
3448 : END MODULE HCOIO_Util_Mod
|
