Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hco_calc_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCO\_Calc\_Mod contains routines to calculate
9 : ! HEMCO core emissions based on the content of the HEMCO EmisList
10 : ! object. All emissions are in [kg/m2/s].
11 : !\\
12 : !\\
13 : ! Emissions for the current datetime are calculated by multiplying base
14 : ! emissions fields with the associated scale factors. Different
15 : ! inventories are merged/overlayed based upon the category and hierarchy
16 : ! attributes assigned to the individual base fields. Within the same
17 : ! category, fields of higher hierarchy overwrite lower-hierarchy fields.
18 : ! Emissions of different categories are always added.
19 : !\\
20 : !\\
21 : ! The assembled emission array is written into the corresponding emission
22 : ! rates array of the HEMCO state object: HcoState%Spc(HcoID)%Emis, where
23 : ! HcoID denotes the corresponding species ID. Emis covers dimension lon,
24 : ! lat, lev on the HEMCO grid, i.e. unlike the emission arrays in EmisList
25 : ! that only cover the levels defined in the source files, Emis extends
26 : ! over all vertical model levels.
27 : !\\
28 : !\\
29 : ! Negative emissions are not supported and are ignored. Surface
30 : ! deposition velocities are stored in HcoState%Spc(HcoID)%Depv and can
31 : ! be added therein.
32 : !\\
33 : !\\
34 : ! In addition to emissions and surface deposition rates, HEMCO also
35 : ! supports concentrations (kg/m3). Data is automatically written into
36 : ! the concentration array HcoState%Spc(HcoID)%Conc if the source data
37 : ! is marked as being concentration data (i.e. if Dta%IsConc is .TRUE.).
38 : ! The latter is automatically determined by HEMCO based upon the data
39 : ! units.
40 : !\\
41 : !\\
42 : ! All emission calculation settings are passed through the HcoState
43 : ! options object (HcoState%Options). These include:
44 : !
45 : ! \begin{itemize}
46 : ! \item ExtNr: extension number to be considered.
47 : ! \item SpcMin: lower species ID (HEMCO ID) to be considered.
48 : ! \item SpcMax: upper species ID (HEMCO ID) to be considered. If set
49 : ! to -1, all species above or equal to SpcMin are considered.
50 : ! \item CatMin: lower emission category to be considered.
51 : ! \item CatMax: upper emission category to be considered. If set to
52 : ! -1, all categories above or equal to CatMin are considered.
53 : ! \item FillBuffer: if set to TRUE, the emissions will be written into
54 : ! buffer array HcoState%Buffer3D instead of HcoState%Spc(ID)%Emis.
55 : ! If this option is enabled, only one species can be calculated at
56 : ! a time (by setting SpcMin/SpcMax, CatMin/CatMax and/or ExtNr
57 : ! accordingly). This option is useful for extensions, e.g. if
58 : ! additional scalings need to be done on some emission fields
59 : ! assembled by HEMCO (e.g. PARANOX extension).
60 : ! \end{itemize}
61 : !
62 : ! !INTERFACE:
63 : !
64 : MODULE HCO_Calc_Mod
65 : !
66 : ! !USES:
67 : !
68 : USE HCO_Diagn_Mod
69 : USE HCO_Error_Mod
70 : USE HCO_Types_Mod
71 : USE HCO_DataCont_Mod, ONLY : Pnt2DataCont
72 :
73 : IMPLICIT NONE
74 : PRIVATE
75 : !
76 : ! !PUBLIC MEMBER FUNCTIONS:
77 : !
78 : PUBLIC :: HCO_CalcEmis
79 : PUBLIC :: HCO_CheckDepv
80 : PUBLIC :: HCO_EvalFld
81 : PUBLIC :: HCO_MaskFld
82 : #ifdef ADJOINT
83 : PUBLIC :: GET_CURRENT_EMISSIONS_ADJ
84 : #endif
85 : !
86 : ! !PRIVATE MEMBER FUNCTIONS:
87 : !
88 : PRIVATE :: GET_CURRENT_EMISSIONS
89 : PRIVATE :: GetMaskVal
90 : PRIVATE :: GetDilFact
91 : PRIVATE :: GetVertIndx
92 : PRIVATE :: GetIdx
93 : !
94 : ! !PARAMETER
95 : !
96 : ! Mask threshold. All mask values below this value will be evaluated
97 : ! as zero (= outside of mask), and all values including and above this
98 : ! value as inside the mask. This is only of relevance if the MaskFractions
99 : ! option is false. If MaskFractions is true, the fractional mask values are
100 : ! considered, e.g. a grid box can contribute 40% to a mask region, etc.
101 : ! The MaskFractions toggle can be set in the settings section of the HEMCO
102 : ! configuration file (Use mask fractions: true/false). It defaults to false.
103 : REAL(sp), PARAMETER :: MASK_THRESHOLD = 0.5_sp
104 : !
105 : ! ============================================================================
106 : !
107 : ! !REVISION HISTORY:
108 : ! 25 Aug 2012 - C. Keller - Initial version.
109 : ! See https://github.com/geoschem/hemco for complete history
110 : !EOP
111 : !------------------------------------------------------------------------------
112 : !BOC
113 : INTERFACE HCO_EvalFld
114 : MODULE PROCEDURE HCO_EvalFld_2D
115 : MODULE PROCEDURE HCO_EvalFld_3D
116 : END INTERFACE HCO_EvalFld
117 :
118 : CONTAINS
119 : !EOC
120 : !------------------------------------------------------------------------------
121 : ! Harmonized Emissions Component (HEMCO) !
122 : !------------------------------------------------------------------------------
123 : !BOP
124 : !
125 : ! !IROUTINE: HCO_CalcEmis
126 : !
127 : ! !DESCRIPTION: Subroutine HCO\_CalcEmis calculates the 3D emission
128 : ! fields at current datetime for the specified species, categories, and
129 : ! extension number.
130 : !\\
131 : !\\
132 : ! !INTERFACE:
133 : !
134 0 : SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC )
135 : !
136 : ! !USES:
137 : !
138 : USE HCO_STATE_MOD, ONLY : HCO_State
139 : USE HCO_ARR_MOD, ONLY : HCO_ArrAssert
140 : USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont
141 : USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined
142 : USE HCO_Scale_Mod, ONLY : HCO_ScaleArr
143 : !
144 : ! !INPUT PARAMETERS:
145 : !
146 : LOGICAL, INTENT(IN ) :: UseConc ! Use concentration fields?
147 : !
148 : ! !INPUT/OUTPUT PARAMETERS:
149 : !
150 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
151 : INTEGER, INTENT(INOUT) :: RC ! Return code
152 : !
153 : ! !REVISION HISTORY:
154 : ! 25 Aug 2012 - C. Keller - Initial Version
155 : ! See https://github.com/geoschem/hemco for complete history
156 : !EOP
157 : !------------------------------------------------------------------------------
158 : !BOC
159 : !
160 : ! !LOCAL VARIABLES:
161 : !
162 : ! Working pointers: list and data container
163 : TYPE(ListCont), POINTER :: Lct
164 : TYPE(DataCont), POINTER :: Dct
165 :
166 : ! Temporary emission arrays
167 : REAL(hp), POINTER :: OutArr(:,:,:) => NULL()
168 : REAL(hp), TARGET :: SpcFlx( HcoState%NX, &
169 : HcoState%NY, &
170 0 : HcoState%NZ )
171 : REAL(hp), TARGET :: CatFlx( HcoState%NX, &
172 : HcoState%NY, &
173 0 : HcoState%NZ )
174 : REAL(hp), TARGET :: TmpFlx( HcoState%NX, &
175 : HcoState%NY, &
176 0 : HcoState%NZ )
177 : REAL(hp) :: Mask ( HcoState%NX, &
178 : HcoState%NY, &
179 0 : HcoState%NZ )
180 : REAL(hp) :: HirFlx( HcoState%NX, &
181 : HcoState%NY, &
182 0 : HcoState%NZ )
183 : REAL(hp) :: HirMsk( HcoState%NX, &
184 : HcoState%NY, &
185 0 : HcoState%NZ )
186 :
187 : ! Integers
188 : INTEGER :: ThisSpc, PrevSpc ! current and previous species ID
189 : INTEGER :: ThisCat, PrevCat ! current and previous category
190 : INTEGER :: ThisHir, PrevHir ! current and previous hierarchy
191 : INTEGER :: SpcMin, SpcMax ! range of species to be considered
192 : INTEGER :: CatMin, CatMax ! range of categories to be considered
193 : INTEGER :: ExtNr ! Extension Nr to be used
194 : INTEGER :: nI, nJ, nL
195 : INTEGER :: nnSpec, FLAG
196 :
197 : LOGICAL :: Found, DoDiagn, EOL, UpdateCat
198 :
199 : ! For error handling & verbose mode
200 : CHARACTER(LEN=255) :: MSG, LOC
201 :
202 : ! testing / debugging
203 : integer :: ix,iy
204 :
205 : !=================================================================
206 : ! HCO_CalcEmis begins here!
207 : !=================================================================
208 :
209 : ! testing only
210 0 : ix = 30
211 0 : iy = 34
212 :
213 : ! Initialize
214 0 : LOC = 'HCO_CalcEmis (HCO_CALC_MOD.F90)'
215 0 : Lct => NULL()
216 0 : Dct => NULL()
217 :
218 : ! Enter routine
219 0 : CALL HCO_ENTER (HcoState%Config%Err, LOC, RC )
220 0 : IF(RC /= HCO_SUCCESS) RETURN
221 :
222 : !-----------------------------------------------------------------
223 : ! Initialize variables
224 : !-----------------------------------------------------------------
225 :
226 : ! Initialize
227 0 : SpcFlx(:,:,:) = 0.0_hp
228 0 : CatFlx(:,:,:) = 0.0_hp
229 0 : HirFlx(:,:,:) = 0.0_hp
230 0 : HirMsk(:,:,:) = 0.0_hp
231 0 : PrevSpc = -1
232 0 : PrevHir = -1
233 0 : PrevCat = -1
234 0 : nnSpec = 0
235 :
236 : ! Pass emission grid dimensions
237 0 : nI = HcoState%NX
238 0 : nJ = HcoState%NY
239 0 : nL = HcoState%NZ
240 :
241 : ! Pass calculation options
242 0 : SpcMin = HcoState%Options%SpcMin !Lower species ID
243 0 : SpcMax = HcoState%Options%SpcMax !Upper species ID
244 0 : CatMin = HcoState%Options%CatMin !Lower emission category
245 0 : CatMax = HcoState%Options%CatMax !Upper emission category
246 0 : ExtNr = HcoState%Options%ExtNr !Extension number
247 0 : DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics?
248 :
249 : ! Verbose mode
250 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
251 0 : WRITE (MSG, *) 'Run HEMCO calculation w/ following options:'
252 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
253 0 : WRITE (MSG, "(A20,I5)") 'Extension number:', ExtNr
254 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
255 0 : WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax
256 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
257 0 : WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax
258 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
259 0 : WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn
260 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
261 : ENDIF
262 :
263 : !=================================================================
264 : ! Walk through all containers of EmisList and determine the
265 : ! emissions for all containers that qualify for calculation.
266 : ! The containers in EmisList are sorted by species, category and
267 : ! hierarchy. This enables a straightforward, piece-by-piece
268 : ! assembly of the final emission array (start with lowest
269 : ! hierarchy emissions, then overwrite piece-by-piece with higher
270 : ! hierarchy values).
271 : !=================================================================
272 :
273 : ! Point to the head of the emissions linked list
274 0 : EOL = .FALSE. ! End of list
275 0 : Lct => NULL()
276 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
277 :
278 : ! Do until end of EmisList (==> loop over all emission containers)
279 : DO
280 : ! Have we reached the end of the list?
281 0 : EOL = .FALSE.
282 0 : IF ( FLAG /= HCO_SUCCESS ) EOL = .TRUE.
283 :
284 : ! ------------------------------------------------------------
285 : ! Select container and update all working variables & arrays.
286 : ! ------------------------------------------------------------
287 : IF ( .NOT. EOL ) THEN
288 :
289 : ! Dct is the current data container
290 0 : Dct => Lct%Dct
291 :
292 : ! Check if this is a base field
293 0 : IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN
294 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
295 0 : CYCLE
296 : ENDIF
297 :
298 : ! Sanity check: Make sure this container holds data.
299 : ! 'Empty' containers are possible if the simulation time
300 : ! is outside of the specified data time range and time
301 : ! slice cycling is deactivated (CycleFlag > 1).
302 0 : IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN
303 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
304 0 : CYCLE
305 : ENDIF
306 :
307 : ! Check if this is the specified extension number
308 0 : IF ( Dct%ExtNr /= ExtNr ) THEN
309 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
310 0 : CYCLE
311 : ENDIF
312 :
313 : ! Advance to next container if the species ID is outside
314 : ! the specified species range (SpcMin - SpcMax). Consider
315 : ! all species above SpcMin if SpcMax is negative!
316 0 : IF( ( Dct%HcoID < SpcMin ) .OR. &
317 : ( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN
318 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
319 0 : CYCLE
320 : ENDIF
321 :
322 : ! Advance to next emission field if the emission category of
323 : ! the current container is outside of the specified species
324 : ! range (CatMin - CatMax). Consider all categories above CatMin
325 : ! if CatMax is negative!
326 0 : IF( ( Dct%Cat < CatMin ) .OR. &
327 : ( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN
328 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
329 0 : CYCLE
330 : ENDIF
331 :
332 : ! Check if this container holds data in the desired unit format,
333 : ! i.e. concentration data if UseConc is enabled, emission data
334 : ! otherwise.
335 0 : IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN
336 0 : CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
337 0 : CYCLE
338 : ENDIF
339 :
340 : ! Update working variables
341 0 : ThisSpc = Dct%HcoID
342 0 : ThisCat = Dct%Cat
343 0 : ThisHir = Dct%Hier
344 :
345 : ! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir
346 : ! to make sure that emissions are added to HEMCO in the section
347 : ! below!
348 : ELSE
349 0 : ThisSpc = -1
350 0 : ThisCat = -1
351 0 : ThisHir = -1
352 : ENDIF
353 :
354 : !--------------------------------------------------------------------
355 : ! Before computing emissions of current data container make sure that
356 : ! emissions of previous container are properly archived.
357 : !--------------------------------------------------------------------
358 :
359 : ! Add emissions on hierarchy level to the category flux array. Do
360 : ! this only if this is a new species, a new category or a new
361 : ! hierarchy level.
362 : ! Note: no need to add to diagnostics because hierarchy level
363 : ! diagnostics are filled right after computing the emissions of
364 : ! a given data container (towards the end of the DO loop).
365 : IF ( (ThisHir /= PrevHir) .OR. &
366 0 : (ThisSpc /= PrevSpc) .OR. &
367 : (ThisCat /= PrevCat) ) THEN
368 :
369 : ! Add hierarchy level emissions to category array over the
370 : ! covered regions.
371 0 : CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx
372 :
373 : ! Reset
374 0 : HirFlx = 0.0_hp
375 0 : HirMsk = 0.0_hp
376 : ENDIF
377 :
378 : !--------------------------------------------------------------------
379 : ! If this is a new species or category, pass the previously collected
380 : ! emissions to the species array. Update diagnostics at category level.
381 : ! Skip this step for first species, i.e. if PrevSpc is still -1.
382 : !--------------------------------------------------------------------
383 0 : UpdateCat = .FALSE.
384 0 : IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE.
385 0 : IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE.
386 0 : IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE.
387 0 : IF ( UpdateCat ) THEN
388 :
389 : ! CatFlx holds the emissions for this category. Pass this to
390 : ! the species array SpcFlx.
391 0 : SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:)
392 :
393 : ! verbose
394 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
395 0 : WRITE(MSG,*) 'Added category emissions to species array: '
396 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
397 0 : WRITE(MSG,*) 'Species : ', PrevSpc
398 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
399 0 : WRITE(MSG,*) 'Category : ', PrevCat
400 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
401 0 : WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx)
402 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
403 0 : WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx)
404 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
405 : ENDIF
406 :
407 : ! Add category emissions to diagnostics at category level
408 : ! (only if defined in the diagnostics list).
409 0 : IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN
410 : ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
411 : ! are updated, including those manually defined in other models
412 : ! (mps, 11/30/21)
413 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
414 : Cat=PrevCat, Hier=-1, HcoID=PrevSpc, &
415 0 : AutoFill=1, Array3D=CatFlx, COL=-1, RC=RC )
416 0 : IF ( RC /= HCO_SUCCESS ) THEN
417 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
418 0 : RETURN
419 : ENDIF
420 : #ifdef ADJOINT
421 : IF (HcoState%IsAdjoint) THEN
422 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
423 : Cat=PrevCat, Hier=-1, HcoID=PrevSpc, &
424 : AutoFill=1, Array3D=CatFlx, &
425 : COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
426 : IF ( RC /= HCO_SUCCESS ) THEN
427 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
428 : RETURN
429 : ENDIF
430 : ENDIF
431 : #endif
432 : ENDIF
433 :
434 : ! Reset CatFlx array and the previously used hierarchy
435 : ! ==> Emission hierarchies are only important within the
436 : ! same category, hence always start over at lowest hierarchy
437 : ! when entering a new category.
438 0 : CatFlx(:,:,:) = 0.0_hp
439 : PrevHir = -1
440 : ENDIF
441 :
442 : !--------------------------------------------------------------------
443 : ! If this is a new species, pass previously calculated emissions
444 : ! to the final emissions array in HcoState.
445 : ! Update diagnostics at extension number level.
446 : ! Don't do before first emission calculation, i.e. if PrevSpc
447 : ! is still the initialized value of -1!
448 : !--------------------------------------------------------------------
449 0 : IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN
450 :
451 : ! Add to OutArr
452 0 : OutArr(:,:,:) = OutArr(:,:,:) + SpcFlx(:,:,:)
453 :
454 : ! testing only
455 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
456 0 : WRITE(MSG,*) 'Added total emissions to output array: '
457 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
458 0 : WRITE(MSG,*) 'Species: ', PrevSpc
459 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
460 0 : WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx)
461 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
462 0 : WRITE(MSG,*) 'OutArr : ', SUM(OutArr)
463 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
464 : ENDIF
465 :
466 : ! Add to diagnostics at extension number level.
467 : ! The same diagnostics may be updated multiple times during
468 : ! the same time step, continuously adding emissions to it.
469 0 : IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN
470 : ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
471 : ! are updated, including those manually defined in other models
472 : ! (mps, 11/30/21)
473 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
474 : Cat=-1, Hier=-1, HcoID=PrevSpc, &
475 0 : AutoFill=1,Array3D=SpcFlx, COL=-1, RC=RC )
476 0 : IF ( RC /= HCO_SUCCESS ) THEN
477 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
478 0 : RETURN
479 : ENDIF
480 : #ifdef ADJOINT
481 : IF (HcoState%IsAdjoint) THEN
482 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
483 : Cat=-1, Hier=-1, HcoID=PrevSpc, &
484 : AutoFill=1,Array3D=SpcFlx, &
485 : COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
486 : IF ( RC /= HCO_SUCCESS ) THEN
487 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
488 : RETURN
489 : ENDIF
490 : ENDIF
491 : #endif
492 : ENDIF
493 :
494 : ! Reset arrays and previous hierarchy.
495 0 : SpcFlx(:,:,:) = 0.0_hp
496 0 : PrevCat = -1
497 0 : PrevHir = -1
498 0 : OutArr => NULL()
499 : ENDIF
500 :
501 : !--------------------------------------------------------------------
502 : ! Exit DO loop here if end of list
503 : !--------------------------------------------------------------------
504 0 : IF ( EOL ) EXIT
505 :
506 : !--------------------------------------------------------------------
507 : ! Update/archive information on species level if needed
508 : !--------------------------------------------------------------------
509 0 : IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN
510 :
511 : ! Update number of species for which emissions have been
512 : ! calculated.
513 0 : nnSpec = nnSpec + 1
514 :
515 : ! To write emissions into temporary array, make OutArr point
516 : ! to the buffer array HcoState%Buffer3D.
517 0 : IF ( HcoState%Options%FillBuffer ) THEN
518 :
519 : ! Cannot use temporary array for more than one species!
520 0 : IF ( nnSpec > 1 ) THEN
521 0 : MSG = 'Cannot fill buffer for more than one species!'
522 0 : CALL HCO_ERROR( MSG, RC )
523 0 : RETURN
524 : ENDIF
525 :
526 : ! Point to array and check allocation status as well as
527 : ! array size.
528 0 : OutArr => HcoState%Buffer3D%Val
529 0 : IF ( .NOT. ASSOCIATED( OutArr ) ) THEN
530 0 : MSG = 'Buffer array is not associated'
531 0 : CALL HCO_ERROR( MSG, RC )
532 0 : RETURN
533 : ENDIF
534 : IF ( (SIZE(OutArr,1) /= nI) .OR. &
535 0 : (SIZE(OutArr,2) /= nJ) .OR. &
536 : (SIZE(OutArr,3) /= nL) ) THEN
537 0 : MSG = 'Buffer array has wrong dimension!'
538 0 : CALL HCO_ERROR( MSG, RC )
539 0 : RETURN
540 : ENDIF
541 :
542 : ! To write emissions directly into HcoState, make OutArr
543 : ! point to current species' array in HcoState. Use emission
544 : ! array for emissions, and concentration array for concentrations.
545 : ELSE
546 :
547 : ! For concentrations:
548 0 : IF ( UseConc ) THEN
549 0 : CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Conc, &
550 0 : nI, nJ, nL, RC )
551 0 : IF ( RC /= HCO_SUCCESS ) THEN
552 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
553 0 : RETURN
554 : ENDIF
555 0 : OutArr => HcoState%Spc(ThisSpc)%Conc%Val
556 :
557 : ! For emissions:
558 : ELSE
559 0 : CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Emis, &
560 0 : nI, nJ, nL, RC )
561 0 : IF ( RC /= HCO_SUCCESS ) THEN
562 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
563 0 : RETURN
564 : ENDIF
565 0 : OutArr => HcoState%Spc(ThisSpc)%Emis%Val
566 : ENDIF
567 :
568 : ENDIF
569 :
570 : ! verbose mode
571 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
572 0 : WRITE(MSG,*) 'Calculating emissions for species ', &
573 0 : TRIM(HcoState%Spc(ThisSpc)%SpcName)
574 0 : CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' )
575 : ENDIF
576 : ENDIF
577 :
578 : !--------------------------------------------------------------------
579 : ! Get current emissions and write into TmpFlx array. The array Mask
580 : ! denotes all valid grid boxes for this inventory.
581 : !--------------------------------------------------------------------
582 0 : TmpFlx(:,:,:) = 0.0_hp
583 0 : CALL GET_CURRENT_EMISSIONS( HcoState, Dct, nI, nJ, nL, TmpFlx, Mask, RC )
584 0 : IF ( RC /= HCO_SUCCESS ) THEN
585 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
586 0 : RETURN
587 : ENDIF
588 :
589 : ! Eventually add universal scale factor
590 0 : CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC )
591 0 : IF ( RC /= HCO_SUCCESS ) THEN
592 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
593 0 : RETURN
594 : ENDIF
595 :
596 : ! Check for negative values according to the corresponding setting
597 : ! in the configuration file: 2 means allow negative values, 1 means
598 : ! set to zero and prompt a warning, else return with error.
599 0 : IF ( HcoState%Options%NegFlag /= 2 ) THEN
600 :
601 0 : IF ( ANY(TmpFlx < 0.0_hp) ) THEN
602 :
603 : ! Set to zero and prompt warning
604 0 : IF ( HcoState%Options%NegFlag == 1 ) THEN
605 0 : WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp
606 0 : MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName)
607 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
608 :
609 : ! Return with error
610 : ELSE
611 : MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // &
612 0 : 'To allow negatives, edit settings in the configuration file.'
613 0 : CALL HCO_ERROR( MSG, RC )
614 0 : RETURN
615 : ENDIF
616 : ENDIF
617 : ENDIF
618 :
619 : ! ------------------------------------------------------------
620 : ! Collect all emissions of the same category (and species) on
621 : ! the hierarchy level into array HirFlx. HirMsk contains the
622 : ! combined covered region. That is, if there are two regional
623 : ! inventories with the same hierarchy HirMsk will cover both
624 : ! of these regions.
625 : ! The specified field hierarchies determine whether the
626 : ! temporary emissions are added (if hierarchy is the same
627 : ! as the previously used hierarchy), or if they overwrite the
628 : ! previous values in HirFlx (if hierarchy is higher than the
629 : ! previous hierarchy).
630 : ! ------------------------------------------------------------
631 :
632 : ! Add emissions to the hierarchy array HirFlx if this hierarchy
633 : ! is the same as previous hierarchy
634 0 : IF ( ThisHir == PrevHir ) THEN
635 0 : HirFlx = HirFlx + TmpFlx
636 0 : HirMsk = HirMsk + Mask
637 :
638 : ! Make sure mask values do not exceed 1.0
639 0 : WHERE(HirMsk > 1.0 ) HirMsk = 1.0
640 :
641 : ! If hierarchy is larger than those of the previously used
642 : ! fields, overwrite HirFlx with new values.
643 : ELSE
644 :
645 0 : HirFlx = TmpFlx
646 0 : HirMsk = Mask
647 :
648 : ENDIF
649 :
650 : ! Update diagnostics at hierarchy level. Make sure that only
651 : ! positive values are used.
652 : ! The same diagnostics may be updated multiple times
653 : ! during the same time step, continuously adding
654 : ! emissions to it.
655 : ! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's
656 : ! ok to keep negative values (ckeller, 7/12/15).
657 0 : IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN
658 : ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
659 : ! are updated, including those manually defined in other models
660 : ! (mps, 11/30/21)
661 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
662 : Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, &
663 : AutoFill=1, Array3D=TmpFlx, &
664 0 : COL=-1, RC=RC )
665 0 : IF ( RC /= HCO_SUCCESS ) THEN
666 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
667 0 : RETURN
668 : ENDIF
669 : #ifdef ADJOINT
670 : IF (HcoState%IsAdjoint) THEN
671 : ! I don't know why I chose collection=-1 instead of
672 : ! collection=HcoState%Diagn%HcoDiagnIDAdjoint like in the other
673 : ! parts of the adjoint code here, but it's what worked in the
674 : ! old repo so I'm keeping it for now. May need to change
675 : CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
676 : Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, &
677 : AutoFill=1, Array3D=TmpFlx, &
678 : COL=-1, RC=RC )
679 : IF ( RC /= HCO_SUCCESS ) THEN
680 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
681 : RETURN
682 : ENDIF
683 : ENDIF
684 :
685 : #endif
686 : ENDIF
687 :
688 : ! Update previously used species, category and hierarchy
689 0 : PrevSpc = ThisSpc
690 0 : PrevCat = ThisCat
691 0 : PrevHir = ThisHir
692 :
693 : ! Advance to next emission container
694 0 : CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG )
695 :
696 : ENDDO ! Loop over EmisList
697 :
698 : ! Make sure internal pointers are nullified
699 0 : Lct => NULL()
700 0 : Dct => NULL()
701 0 : OutArr => NULL()
702 :
703 : ! verbose
704 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
705 0 : WRITE (MSG, *) 'HEMCO emissions successfully calculated!'
706 0 : CALL HCO_MSG ( HcoState%Config%Err, MSG )
707 : ENDIF
708 :
709 : ! Leave w/ success
710 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
711 :
712 : END SUBROUTINE HCO_CalcEmis
713 : !EOC
714 : !------------------------------------------------------------------------------
715 : ! Harmonized Emissions Component (HEMCO) !
716 : !------------------------------------------------------------------------------
717 : !BOP
718 : !
719 : ! !IROUTINE: HCO_CheckDepv
720 : !
721 : ! !DESCRIPTION: Subroutine HCO\_CheckDepv is a simple routine to check the
722 : ! dry deposition frequency value. This is to avoid unrealistically high
723 : ! deposition frequencies that may occur if grid box concentrations are very
724 : ! low. The deposition frequency is limited to a value that will make sure
725 : ! that the drydep exponent ( exp( -depfreq * dt ) ) is still small enough to
726 : ! remove all species mass. The maximum limit of depfreq * dt can be defined
727 : ! as a HEMCO option (MaxDepExp). Its default value is 20.0.
728 : !\\
729 : !\\
730 : ! !INTERFACE:
731 : !
732 0 : SUBROUTINE HCO_CheckDepv( HcoState, Depv, RC )
733 : !
734 : ! !USES:
735 : !
736 : USE HCO_STATE_MOD, ONLY : HCO_State
737 : !
738 : ! !INPUT/OUTPUT PARAMETERS:
739 : !
740 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
741 : REAL(hp), INTENT(INOUT) :: Depv ! Deposition velocity
742 : INTEGER, INTENT(INOUT) :: RC ! Return code
743 : !
744 : ! !REVISION HISTORY:
745 : ! 11 Mar 2015 - C. Keller - Initial Version
746 : ! See https://github.com/geoschem/hemco for complete history
747 : !EOP
748 : !------------------------------------------------------------------------------
749 : !BOC
750 : !
751 : ! !LOCAL VARIABLES:
752 : !
753 : REAL(hp) :: ExpVal
754 :
755 : !=================================================================
756 : ! HCO_CheckDepv begins here!
757 : !=================================================================
758 :
759 0 : ExpVal = Depv * HcoState%TS_EMIS
760 0 : IF ( ExpVal > HcoState%Options%MaxDepExp ) THEN
761 0 : Depv = HcoState%Options%MaxDepExp / HcoState%TS_EMIS
762 : ENDIF
763 :
764 0 : END SUBROUTINE HCO_CheckDepv
765 : !EOC
766 : !------------------------------------------------------------------------------
767 : ! Harmonized Emissions Component (HEMCO) !
768 : !------------------------------------------------------------------------------
769 : !BOP
770 : !
771 : ! !IROUTINE: Get_Current_Emissions
772 : !
773 : ! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current
774 : ! emissions for the specified emission container.
775 : ! This subroutine is only called by HCO\_CalcEmis and for base emission
776 : ! containers, i.e. containers of type 1.
777 : !\\
778 : !\\
779 : ! !INTERFACE:
780 : !
781 0 : SUBROUTINE Get_Current_Emissions( HcoState, BaseDct, nI, nJ, &
782 0 : nL, OUTARR_3D, MASK, RC, UseLL )
783 : !
784 : ! !USES:
785 : !
786 : USE HCO_State_Mod, ONLY : HCO_State
787 : USE HCO_tIdx_MOD, ONLY : tIDx_GetIndx
788 : USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined
789 : !
790 : ! !INPUT PARAMETERS:
791 : !
792 : INTEGER, INTENT(IN) :: nI ! # of lons
793 : INTEGER, INTENT(IN) :: nJ ! # of lats
794 : INTEGER, INTENT(IN) :: nL ! # of levs
795 : !
796 : ! !INPUT/OUTPUT PARAMETERS:
797 : !
798 :
799 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
800 : TYPE(DataCont), POINTER :: BaseDct ! base emission
801 : ! container
802 : REAL(hp), INTENT(INOUT) :: outArr_3D(nI,nJ,nL) ! output array
803 : REAL(hp), INTENT(INOUT) :: mask (nI,nJ,nL) ! mask array
804 : INTEGER, INTENT(INOUT) :: RC
805 : !
806 : ! !OUTPUT PARAMETERS:
807 : !
808 : INTEGER, INTENT( OUT), OPTIONAL :: UseLL
809 : !
810 : ! !REMARKS:
811 : ! This routine uses multiple loops over all grid boxes (base emissions
812 : ! and scale factors use separate loops). In an OMP environment, this approach
813 : ! seems to be faster than using only one single loop (but repeated calls to
814 : ! point to containers, etc.). The alternative approach is used in routine
815 : ! Get\_Current\_Emissions\_B at the end of this module and may be employed
816 : ! on request.
817 : !
818 : ! !REVISION HISTORY:
819 : ! 25 Aug 2012 - C. Keller - Initial Version
820 : ! See https://github.com/geoschem/hemco for complete history
821 : !EOP
822 : !------------------------------------------------------------------------------
823 : !BOC
824 : !
825 : ! !LOCAL VARIABLES:
826 : !
827 : ! Pointers
828 : TYPE(DataCont), POINTER :: ScalDct
829 : TYPE(DataCont), POINTER :: MaskDct
830 : TYPE(DataCont), POINTER :: LevDct1
831 : TYPE(DataCont), POINTER :: LevDct2
832 :
833 : ! Scalars
834 : REAL(sp) :: tmpVal, MaskScale
835 : REAL(hp) :: outData
836 : REAL(hp) :: ScalFact
837 : INTEGER :: tIDx, IDX
838 : INTEGER :: totLL, nnLL
839 : INTEGER :: I, J, L, N
840 : INTEGER :: LowLL, UppLL, ScalLL, TmpLL
841 : INTEGER :: EC, ERROR
842 : CHARACTER(LEN=255) :: MSG, thisLoc
843 : LOGICAL :: NegScalExist
844 : LOGICAL :: MaskFractions
845 : LOGICAL :: isLevDct1
846 : LOGICAL :: isLevDct2
847 : LOGICAL :: isMaskDct
848 : LOGICAL :: isPblHt
849 : LOGICAL :: isBoxHt
850 : INTEGER :: LevDct1_Unit
851 : INTEGER :: LevDct2_Unit
852 :
853 : ! testing only
854 : INTEGER, PARAMETER :: IX=25, IY=25
855 :
856 : !=================================================================
857 : ! GET_CURRENT_EMISSIONS begins here
858 : !=================================================================
859 :
860 : ! Initialize
861 0 : ScalDct => NULL()
862 0 : MaskDct => NULL()
863 0 : msg = ''
864 0 : thisLoc = 'GET_CURRENT_EMISSIONS (hco_calc_mod.F90)'
865 :
866 : ! Enter
867 0 : CALL HCO_Enter( HcoState%Config%Err, thisLoc, RC )
868 0 : IF ( RC /= HCO_SUCCESS ) RETURN
869 :
870 : ! Check if container contains data
871 0 : IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN
872 0 : msg = 'Array not defined: ' // TRIM(BaseDct%cName)
873 0 : CALL HCO_Error( msg, RC, thisLoc )
874 0 : RETURN
875 : ENDIF
876 :
877 : ! Initialize mask. By default, assume that we use all grid boxes.
878 0 : MASK = 1.0_hp
879 0 : MaskFractions = HcoState%Options%MaskFractions
880 :
881 : ! Verbose output
882 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
883 0 : WRITE(msg,*) 'Evaluate field ', TRIM(BaseDct%cName)
884 0 : CALL HCO_Msg( HcoState%Config%Err, msg, sep1=' ' )
885 : ENDIF
886 :
887 : #if !defined ( ESMF_ )
888 : ! Put check for PBLHEIGHT here, if not running in ESMF/MAPL
889 0 : IF ( .NOT. ASSOCIATED( HcoState%Grid%PBLHEIGHT%Val ) ) THEN
890 0 : msg = 'PBLHEIGHT (in meters) is missing in HEMCO state'
891 0 : CALL HCO_Error( msg, RC, thisLoc )
892 0 : RETURN
893 : ENDIF
894 : #endif
895 :
896 : !========================================================================
897 : ! Set base emissions
898 : !========================================================================
899 :
900 : ! Check for level index containers
901 : ! Move error checks here, outside of the parallel DO loop
902 : ! Remove ELSE blocks for efficiency
903 0 : LevDct1 => NULL()
904 0 : IF ( BaseDct%levScalID1 > 0 ) THEN
905 0 : CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC )
906 0 : IF ( RC /= HCO_SUCCESS ) THEN
907 0 : msg = 'Could not get a pointer to LevDct1 from HcoState!'
908 0 : CALL HCO_Error( msg, RC, thisLoc )
909 0 : RETURN
910 : ENDIF
911 : ENDIF
912 :
913 0 : LevDct2 => NULL()
914 0 : IF ( BaseDct%levScalID2 > 0 ) THEN
915 0 : CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC )
916 0 : IF ( RC /= HCO_SUCCESS ) THEN
917 0 : msg = 'Could not get a pointer to LevDct2 in HcoState!'
918 0 : CALL HCO_Error( msg, RC, thisLoc )
919 0 : RETURN
920 : ENDIF
921 : ENDIF
922 :
923 : ! Test whether LevDct1 and LevDct2 are associated
924 0 : isLevDct1 = ASSOCIATED( LevDct1 )
925 0 : isLevDct2 = ASSOCIATED( LevDct2 )
926 :
927 : ! Get the units of LevDct1 (if it exists)
928 0 : LevDct1_Unit = -1
929 0 : IF ( isLevDct1 ) THEN
930 0 : LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
931 0 : IF ( LevDct1_Unit < 0 ) THEN
932 0 : MSG = 'LevDct1 units are not defined!'
933 0 : CALL HCO_Error( msg, RC, thisLoc )
934 0 : RETURN
935 : ENDIF
936 : ENDIF
937 :
938 : ! Get the units of LevDct2 (if it exists)
939 0 : LevDct2_Unit = -1
940 0 : IF ( isLevDct2 ) THEN
941 0 : LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
942 0 : IF ( LevDct2_Unit < 0 ) THEN
943 0 : MSG = 'LevDct2_Units are not defined!'
944 0 : CALL HCO_Error( msg, RC, thisLoc )
945 0 : RETURN
946 : ENDIF
947 : ENDIF
948 :
949 : ! Throw an error if boxheight is missing and the units are in meters
950 0 : IF ( LevDct1_Unit == HCO_EMISL_M .or. &
951 : LevDct2_Unit == HCO_EMISL_M ) THEN
952 0 : IF ( .not. ASSOCIATED( HcoState%Grid%BXHEIGHT_M%Val ) ) THEN
953 0 : MSG = 'Boxheight (in meters) is missing in HEMCO state'
954 0 : CALL HCO_Error( msg, RC, thisLoc )
955 0 : RETURN
956 : ENDIF
957 : ENDIF
958 :
959 : ! Throw an error if boxheight is missing and the units are in PBL frac
960 0 : IF ( LevDct1_Unit == HCO_EMISL_PBL .or. &
961 : LevDct2_Unit == HCO_EMISL_PBL ) THEN
962 0 : IF ( .not. ASSOCIATED( HcoState%Grid%PBLHEIGHT%Val ) ) THEN
963 0 : MSG = 'Boundary layer height is missing in HEMCO state'
964 0 : CALL HCO_Error( msg, RC, thisLoc )
965 0 : RETURN
966 : ENDIF
967 : ENDIF
968 :
969 : ! Initialize non-private loop variables here
970 0 : error = 0
971 0 : totLL = 0.0
972 0 : nnLL = 0.0
973 :
974 : ! Loop over all latitudes and longitudes
975 : !$OMP PARALLEL DO &
976 : !$OMP DEFAULT( SHARED )&
977 : !$OMP PRIVATE( I, J, L, tIdx, tmpVal, LowLL, UppLL, EC )&
978 : !$OMP REDUCTION( +:totLL )&
979 : !$OMP REDUCTION( +:nnLL )&
980 : !$OMP COLLAPSE( 2 )&
981 : !$OMP SCHEDULE( DYNAMIC, 8 )
982 0 : DO J = 1, nJ
983 0 : DO I = 1, nI
984 :
985 : ! Continue to end of loop if an error has occurred
986 : ! (we cannot exit from a parallel loop)
987 0 : IF ( error > 0 ) CYCLE
988 :
989 : ! Zero private variables for safety's sake
990 0 : tmpVal = 0.0_hp
991 0 : lowLL = 0
992 0 : uppLL = 0
993 0 : EC = HCO_SUCCESS
994 :
995 : !---------------------------------------------------------------------
996 : ! Get current time index for this container and at this location
997 : !---------------------------------------------------------------------
998 0 : tIDx = tIDx_GetIndx( HcoState, BaseDct%Dta, I, J )
999 0 : IF ( tIDx < 1 ) THEN
1000 0 : WRITE( msg, * ) 'Cannot get time slice index at location ',I,J, &
1001 0 : ': ', TRIM(BaseDct%cName), tIDx
1002 0 : error = 1
1003 0 : CYCLE
1004 : ENDIF
1005 :
1006 : !---------------------------------------------------------------------
1007 : ! Get lower and upper vertical index
1008 : !---------------------------------------------------------------------
1009 : CALL GetVertIndx( HcoState, BaseDct, isLevDct1, LevDct1, &
1010 : LevDct1_Unit, isLevDct2, LevDct2, LevDct2_Unit, &
1011 : I, J, LowLL, UppLL, &
1012 0 : RC=EC )
1013 :
1014 : ! Trap error
1015 0 : IF ( EC /= HCO_SUCCESS ) THEN
1016 0 : WRITE(MSG,*) 'Error getting vertical index at location ',I,J,&
1017 0 : ': ', TRIM(BaseDct%cName)
1018 0 : error = 1
1019 0 : CYCLE
1020 : ENDIF
1021 :
1022 : !---------------------------------------------------------------------
1023 : ! Apply vertical dilution factor (if necessary)
1024 : !---------------------------------------------------------------------
1025 :
1026 : ! Update variables for computing the average level
1027 0 : totLL = totLL + UppLL
1028 0 : nnLL = nnLL + 1
1029 :
1030 : ! Loop over all levels
1031 0 : DO L = LowLL, UppLL
1032 :
1033 : ! Get base value. Use uniform value if scalar field.
1034 0 : tmpVal = Get_Value_From_DataCont( I, J, L, tIdx, BaseDct )
1035 :
1036 : ! If it's a missing value, mask box as unused
1037 : ! and set value to zero.
1038 0 : IF ( tmpVal == HCO_MISSVAL ) THEN
1039 0 : mask(I,J,:) = 0.0_hp
1040 0 : outArr_3D(I,J,L) = 0.0_hp
1041 0 : CYCLE
1042 : ENDIF
1043 :
1044 : ! Otherwise, apply the vertical dilution factor (if necessary)
1045 : CALL Apply_Dilution_Factor( I = I, &
1046 : J = J, &
1047 : L = L, &
1048 : LowLL = LowLL, &
1049 : UppLL = UppLL, &
1050 : HcoState = HcoState, &
1051 : BaseDct = BaseDct, &
1052 : inData = tmpVal, &
1053 0 : outData = outArr_3D(I,J,L), &
1054 0 : RC = EC )
1055 :
1056 : ! Trap errors
1057 0 : IF ( EC /= HCO_SUCCESS ) THEN
1058 0 : error = 1
1059 0 : msg = 'Error encountered in routine "Apply_Dilution_Factor"!'
1060 0 : CYCLE
1061 : ENDIF
1062 :
1063 : ENDDO !L
1064 :
1065 : ENDDO !I
1066 : ENDDO !J
1067 : !$OMP END PARALLEL DO
1068 :
1069 : !------------------------------------------------------------------------
1070 : ! Return if an error occurred in the parallel loop above
1071 : !------------------------------------------------------------------------
1072 0 : IF ( error == 1 ) THEN
1073 0 : CALL HCO_Error( msg, RC, thisLoc )
1074 0 : RETURN
1075 : ENDIF
1076 :
1077 : !========================================================================
1078 : ! Apply scale factors to base emissions
1079 : !
1080 : ! The container IDs of all scale factors associated with this base
1081 : ! container are stored in vector Scal_cID.
1082 : !========================================================================
1083 0 : IF ( BaseDct%nScalID > 0 ) THEN
1084 :
1085 : ! Loop over all scale factors for this base emissions
1086 0 : DO N = 1, BaseDct%nScalID
1087 :
1088 : !------------------------------------------------------------------
1089 : ! Get a pointer to the data container holding this scale factor
1090 : !------------------------------------------------------------------
1091 :
1092 : ! Get the Nth scale factor container ID
1093 0 : IDX = BaseDct%Scal_cID(N)
1094 :
1095 : ! Point to data container with the given container ID
1096 0 : CALL Pnt2DataCont( HcoState, IDX, ScalDct, RC )
1097 0 : IF ( RC /= HCO_SUCCESS ) THEN
1098 : msg = 'Could not get scale factor for base emissions field ' // &
1099 0 : TRIM( BaseDct%cName )
1100 0 : CALL HCO_Error( msg, RC, thisLoc )
1101 0 : RETURN
1102 : ENDIF
1103 :
1104 : ! Sanity check: scale field cannot be a base field
1105 0 : IF ( ScalDct%DctType == HCO_DCTTYPE_BASE ) THEN
1106 0 : msg = 'Wrong scale field type: ' // TRIM(ScalDct%cName)
1107 0 : CALL HCO_Error( msg, RC, thisLoc )
1108 0 : RETURN
1109 : ENDIF
1110 :
1111 : !------------------------------------------------------------------
1112 : ! Skip this scale factor if no data defined. This is possible
1113 : ! if scale factors are only defined for a given time range and
1114 : ! the simulation datetime is outside of this range.
1115 : !------------------------------------------------------------------
1116 0 : IF ( .not. FileData_ArrIsDefined( ScalDct%Dta ) ) THEN
1117 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1118 : msg = 'Skip scale factor ' // TRIM( ScalDct%cName ) // &
1119 0 : ' because it is not defined for this datetime.'
1120 0 : CALL HCO_MSG( HcoState%Config%Err, msg )
1121 : ENDIF
1122 : CYCLE
1123 : ENDIF
1124 :
1125 : ! Verbose printout
1126 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
1127 0 : MSG = 'Applying scale factor ' // TRIM(ScalDct%cName)
1128 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1129 : ENDIF
1130 :
1131 : !------------------------------------------------------------------
1132 : ! Get vertical extension of this scale factor array.
1133 : !------------------------------------------------------------------
1134 0 : ScalLL = 1
1135 0 : IF ( ScalDct%Dta%SpaceDim == 3 ) THEN
1136 0 : ScalLL = SIZE( ScalDct%Dta%V3(1)%Val, 3 )
1137 : ENDIF
1138 :
1139 : !------------------------------------------------------------------
1140 : ! Check if there is a mask field associated with this scale
1141 : ! factor. In this case, get a pointer to the corresponding
1142 : ! mask field and evaluate scale factors only inside the mask
1143 : ! region.
1144 : !------------------------------------------------------------------
1145 0 : IF ( ASSOCIATED(ScalDct%Scal_cID) ) THEN
1146 0 : CALL Pnt2DataCont( HcoState, ScalDct%Scal_cID(1), MaskDct, RC )
1147 0 : IF ( RC /= HCO_SUCCESS ) THEN
1148 : msg = 'Could not get mask field for scale factor: ' // &
1149 0 : TRIM( ScalDct%cName )
1150 0 : CALL HCO_Error( msg, RC, thisLoc )
1151 0 : RETURN
1152 : ENDIF
1153 :
1154 : ! Sanity check: The data container MaskDct must be a mask!
1155 0 : IF ( MaskDct%DctType /= HCO_DCTTYPE_MASK ) THEN
1156 : MSG = 'Invalid mask for scale factor: ' // &
1157 : TRIM( ScalDct%cName ) // '; mask: ' // &
1158 0 : TRIM( MaskDct%cName )
1159 0 : CALL HCO_Error( msg, RC, thisLoc )
1160 0 : RETURN
1161 : ENDIF
1162 : ENDIF
1163 :
1164 : ! Set a flag to denote whether MaskDct is associated
1165 : ! This can be done outside of the parallel loops below
1166 0 : isMaskDct = ASSOCIATED( MaskDct )
1167 :
1168 : !------------------------------------------------------------------
1169 : ! Apply the mask to the scale factor
1170 : !------------------------------------------------------------------
1171 :
1172 : ! Reinitialize error flag. Will be set to > 0 if error occurs,
1173 : ! and to -1 if negative scale factor is ignored.
1174 0 : error = 0
1175 :
1176 : ! Loop over all latitudes and longitudes
1177 : !$OMP PARALLEL DO &
1178 : !$OMP DEFAULT( SHARED )&
1179 : !$OMP PRIVATE( I, J, tIdx, tmpVal, L )&
1180 : !$OMP PRIVATE( LowLL, UppLL, tmpLL, MaskScale, EC )&
1181 : !$OMP COLLAPSE( 2 )&
1182 : !$OMP SCHEDULE( DYNAMIC, 8 )
1183 0 : DO J = 1, nJ
1184 0 : DO I = 1, nI
1185 :
1186 : ! Continue to end of loop if an error has occurred
1187 : ! (we cannot exit from a parallel loop)
1188 0 : IF ( error > 0 ) CYCLE
1189 :
1190 : !---------------------------------------------------------------
1191 : ! If there is a mask associated with this scale factors, check
1192 : ! if this grid box is within or outside of the mask region.
1193 : ! Values that partially fall into the mask region are either
1194 : ! treated as binary (100% inside or outside), or partially
1195 : ! (using the real grid area fractions), depending on the
1196 : ! HEMCO options.
1197 : !---------------------------------------------------------------
1198 :
1199 : ! Default mask scaling is 1.0 (no mask applied)
1200 0 : maskScale = 1.0_sp
1201 :
1202 : ! If there is a mask applied to this scale factor ...
1203 0 : IF ( isMaskDct ) THEN
1204 : CALL GetMaskVal ( MaskDct, I, J, &
1205 0 : MaskScale, MaskFractions, EC )
1206 0 : IF ( EC /= HCO_SUCCESS ) THEN
1207 : error = 4
1208 : CYCLE
1209 : ENDIF
1210 : ENDIF
1211 :
1212 : ! We continue an skip this grid box if mask is completely zero
1213 0 : IF ( maskScale <= 0.0_sp ) CYCLE
1214 :
1215 : ! Get current time index for this container and at this location
1216 0 : tIDx = tIDx_GetIndx( HcoState, ScalDct%Dta, I, J )
1217 0 : IF ( tIDx < 1 ) THEN
1218 0 : WRITE(*,*) 'Cannot get time slice index at location ',I,J, &
1219 0 : ': ', TRIM(ScalDct%cName), tIDx
1220 0 : error = 3
1221 0 : CYCLE
1222 : ENDIF
1223 :
1224 : !---------------------------------------------------------------
1225 : ! Check if this is a mask. If so, add mask values to the MASK
1226 : ! array. For now, we assume masks to be binary, i.e. 0 or 1.
1227 : ! We may want to change that in future to also support values
1228 : ! in between. This is especially important when regridding
1229 : ! high resolution masks onto coarser grids!
1230 : !---------------------------------------------------------------
1231 0 : IF ( ScalDct%DctType == HCO_DCTTYPE_MASK ) THEN
1232 :
1233 : ! Get mask value
1234 0 : CALL GetMaskVal( ScalDct, I, J, TMPVAL, MaskFractions, EC )
1235 0 : IF ( EC /= HCO_SUCCESS ) THEN
1236 : error = 4
1237 : CYCLE
1238 : ENDIF
1239 :
1240 : ! Pass to output mask
1241 0 : mask(I,J,:) = mask(I,J,:) * TMPVAL
1242 :
1243 : ! Verbose printout
1244 : IF ( HCO_IsVerb(HcoState%Config%Err) .and. &
1245 0 : I==1 .AND. J==1 ) THEN
1246 : msg = 'Mask field ' // TRIM( ScalDct%cName ) // &
1247 0 : ' found and added to temporary mask.'
1248 0 : CALL HCO_Msg( HcoState%Config%Err, msg )
1249 : ENDIF
1250 :
1251 : ! Advance to next grid box
1252 : CYCLE
1253 : ENDIF
1254 :
1255 : !---------------------------------------------------------------
1256 : ! For non-mask fields, apply scale factors to all levels
1257 : ! of the base field individually. If the scale factor
1258 : ! field has more than one vertical level, use the
1259 : ! vertical level closest to the corresponding vertical
1260 : ! level of the base emission field
1261 : !---------------------------------------------------------------
1262 :
1263 : ! Get lower and upper vertical index
1264 : CALL GetVertIndx( HcoState, BaseDct, isLevDct1, &
1265 : LevDct1, LevDct1_Unit, isLevDct2, &
1266 : LevDct2, LevDct2_Unit, I, &
1267 0 : J, LowLL, UppLL, EC )
1268 :
1269 : ! Trap errors
1270 0 : IF ( EC /= HCO_SUCCESS ) THEN
1271 : error = 1
1272 : CYCLE
1273 : ENDIF
1274 :
1275 : ! Loop over all vertical levels of the base field
1276 0 : DO L = LowLL,UppLL
1277 :
1278 : ! If the vertical level exceeds the number of available
1279 : ! scale factor levels, use the highest available level.
1280 : ! Otherwise use the same vertical level index.
1281 0 : TmpLL = L
1282 0 : IF ( L > ScalLL ) TmpLL = ScalLL
1283 :
1284 : !------------------------------------------------------------
1285 : ! Get scale factor for this grid box. Use same uniform
1286 : ! value if it's a scalar field.
1287 : !------------------------------------------------------------
1288 0 : tmpVal = Get_Value_From_DataCont( I, J, L, tIdX, ScalDct )
1289 :
1290 : ! Set missing value to one
1291 0 : IF ( tmpVal == HCO_MISSVAL ) tmpVal = 1.0_sp
1292 :
1293 : ! Eventually apply mask scaling
1294 0 : IF ( maskScale /= 1.0_sp ) THEN
1295 0 : tmpVal = tmpVal * MaskScale
1296 : ENDIF
1297 :
1298 : !------------------------------------------------------------
1299 : ! Negative scale factors: proceed according to the "negative
1300 : ! value" setting specified in the HEMCO configuration file
1301 : ! NegFlag = 2: Use this value (i.e. pass thru this IF stmt)
1302 : ! NegFlag = 1: Ignore ands how warning
1303 : ! NegFlag = 0: Return with error
1304 : !------------------------------------------------------------
1305 0 : IF ( tmpVal < 0.0_sp .and. HcoState%Options%NegFlag /= 2 ) THEN
1306 :
1307 : ! NegFlag = 1: ignore and show warning
1308 : ! Otherwise return with error
1309 0 : IF ( HcoState%Options%NegFlag == 1 ) THEN
1310 : error = -1
1311 : CYCLE
1312 : ELSE
1313 0 : WRITE( 6, * ) 'Negative scale factor at ', &
1314 0 : I, J, TmpLL, tidx, ': ', &
1315 0 : TRIM(ScalDct%cName), TMPVAL
1316 0 : error = 1
1317 0 : CYCLE
1318 : ENDIF
1319 : ENDIF
1320 :
1321 : !------------------------------------------------------------
1322 : ! Apply scale factor in accordance to field operator
1323 : !------------------------------------------------------------
1324 : CALL Apply_Scale_Factor( I = I, &
1325 : J = J, &
1326 : L = L, &
1327 : ScalDct = ScalDct, &
1328 : scalFac = tmpVal, &
1329 0 : dataVal = outArr_3D(I,J,L), &
1330 0 : RC = EC )
1331 :
1332 : ! Trap errors
1333 0 : IF ( EC /= HCO_SUCCESS ) THEN
1334 0 : error = 2
1335 0 : CYCLE
1336 : ENDIF
1337 :
1338 : ENDDO
1339 :
1340 : !---------------------------------------------------------------
1341 : ! Verbose printout
1342 : !----------------------------------------------------------------
1343 : IF ( HCO_IsVerb(HcoState%Config%Err) .and. &
1344 0 : I == ix .and. J == iy ) THEN
1345 0 : write(MSG,*) 'Scale field ', TRIM(ScalDct%cName)
1346 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1347 0 : write(MSG,*) 'Time slice: ', tIdx
1348 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1349 0 : write(MSG,*) 'IX, IY: ', IX, IY
1350 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1351 0 : write(MSG,*) 'Scale factor (IX,IY,L1): ', TMPVAL
1352 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1353 0 : write(MSG,*) 'Mathematical operation : ', ScalDct%Oper
1354 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1355 : ENDIF
1356 :
1357 : ENDDO !I
1358 : ENDDO !J
1359 : !$OMP END PARALLEL DO
1360 :
1361 : !------------------------------------------------------------------
1362 : ! Return with error if an error was encountered in the loop above
1363 : !------------------------------------------------------------------
1364 0 : IF ( error > 0 ) THEN
1365 :
1366 : ! Construct the appropriate error message
1367 0 : SELECT CASE( error )
1368 : CASE( 1 )
1369 0 : msg = 'Negative scale factor found (aborted): '
1370 : CASE( 2 )
1371 0 : msg = 'Illegal mathematical operator for scale factor: '
1372 : CASE( 3 )
1373 0 : msg = 'Encountered negative time index for scale factor: '
1374 : CASE( 4 )
1375 0 : msg = 'Error applying mask to scale factor: '
1376 : CASE DEFAULT
1377 0 : msg = 'Error when applying scale factor: '
1378 : END SELECT
1379 :
1380 : ! Append name of scale factor to error message
1381 0 : msg = TRIM( msg ) // TRIM( ScalDct%cName )
1382 :
1383 : ! Exit with error message
1384 0 : ScalDct => NULL()
1385 0 : CALL HCO_Error( msg, RC, thisLoc )
1386 0 : RETURN
1387 : ENDIF
1388 :
1389 : ! eventually prompt warning for negative values
1390 0 : IF ( ERROR == -1 ) THEN
1391 : msg = 'Negative scale factor found (ignored): ' // &
1392 0 : TRIM( ScalDct%cName )
1393 0 : CALL HCO_WARNING( HcoState%Config%Err, msg, RC )
1394 : ENDIF
1395 :
1396 : ! Free pointer
1397 0 : MaskDct => NULL()
1398 :
1399 : ENDDO ! N
1400 : ENDIF ! N > 0
1401 :
1402 : !========================================================================
1403 : ! Update optional variables
1404 : !========================================================================
1405 0 : IF ( PRESENT( UseLL) ) THEN
1406 0 : UseLL = 1
1407 0 : IF ( nnLL > 0 ) UseLL = NINT(REAL(TotLL,kind=sp)/REAL(nnLL,kind=sp))
1408 : ENDIF
1409 :
1410 : ! Weight output emissions by mask
1411 0 : outArr_3D = outArr_3D * mask
1412 :
1413 : ! Cleanup and leave w/ success
1414 0 : ScalDct => NULL()
1415 0 : CALL HCO_Leave( HcoState%Config%Err, RC )
1416 :
1417 : END SUBROUTINE Get_Current_Emissions
1418 : !EOC
1419 : !------------------------------------------------------------------------------
1420 : ! GEOS-Chem Global Chemical Transport Model !
1421 : !------------------------------------------------------------------------------
1422 : !BOP
1423 : !
1424 : ! !IROUTINE: Get_Value_From_DataCont
1425 : !
1426 : ! !DESCRIPTION: Returns a data value stored in a data container at a given
1427 : ! grid box and time.
1428 : !\\
1429 : !\\
1430 : ! !INTERFACE:
1431 : !
1432 0 : FUNCTION Get_Value_From_DataCont( I, J, L, tIdx, Dct ) RESULT( dataVal )
1433 : !
1434 : ! !INPUT PARAMETERS:
1435 : !
1436 : INTEGER, INTENT(IN) :: I ! Longitude (or X-dim) index
1437 : INTEGER, INTENT(IN) :: J ! Latitude (or Y-dim) index
1438 : INTEGER, INTENT(IN) :: L ! Vertical level index
1439 : INTEGER, INTENT(IN) :: tIdx ! Time index
1440 : TYPE(DataCont), POINTER :: Dct ! Data Container object
1441 : !
1442 : ! !RETURN VALUE:
1443 : !
1444 : REAL(sp) :: dataVal ! Data at this grid box and time
1445 : !
1446 : ! !REMARKS:
1447 : ! This code was abstracted out of Get_Current_Emisssions for clarity.
1448 : ! We have refactored the code to remove ELSE blocks for better efficiency.
1449 : !EOP
1450 : !------------------------------------------------------------------------------
1451 : !BOC
1452 :
1453 : ! Data is a 1-D scaler: Return a uniform value
1454 0 : IF ( Dct%Dta%SpaceDim == 1 ) THEN
1455 0 : dataVal = Dct%Dta%V2(tIDx)%Val(1,1)
1456 0 : RETURN
1457 : ENDIF
1458 :
1459 : ! Data is a 2-D array: Return value at (I,J)
1460 0 : IF ( Dct%Dta%SpaceDim == 2 ) THEN
1461 0 : dataVal = Dct%Dta%V2(tIDx)%Val(I,J)
1462 0 : RETURN
1463 : ENDIF
1464 :
1465 : ! Data is a 3-D array: Return value at (I,J,L)
1466 0 : dataVal = Dct%Dta%V3(tIDx)%Val(I,J,L)
1467 :
1468 0 : END FUNCTION Get_Value_From_DataCont
1469 : !EOC
1470 : !------------------------------------------------------------------------------
1471 : ! GEOS-Chem Global Chemical Transport Model !
1472 : !------------------------------------------------------------------------------
1473 : !BOP
1474 : !
1475 : ! !IROUTINE: Apply_Dilution_Factor
1476 : !
1477 : ! !DESCRIPTION: Applies the dilution factor to input data. This algorithm
1478 : ! has been abstracted out of Get_Current_Emissions for computational
1479 : ! efficiency.
1480 : !\\
1481 : !\\
1482 : ! !INTERFACE:
1483 : !
1484 0 : SUBROUTINE Apply_Dilution_Factor( I, J, L, LowLL, &
1485 : UppLL, HcoState, BaseDct, inData, &
1486 : outData, RC )
1487 : !
1488 : ! !USES:
1489 : !
1490 : USE HCO_Error_Mod
1491 : USE HCO_State_Mod, ONLY : HCO_State
1492 : !
1493 : ! !INPUT PARAMETERS:
1494 : !
1495 : INTEGER, INTENT(IN) :: I ! Longitude (or X-dim) index
1496 : INTEGER, INTENT(IN) :: J ! Latitude (or Y-dim) index
1497 : INTEGER, INTENT(IN) :: L ! Vertical level index
1498 : INTEGER, INTENT(IN) :: LowLL ! Lower level for emissions
1499 : INTEGER, INTENT(IN) :: UppLL ! Upper level for emissions
1500 : TYPE(DataCont), POINTER :: BaseDct ! Base data container
1501 : REAL(sp), INTENT(IN) :: inData ! Input data
1502 : !
1503 : ! !INPUT/OUTPUT PARAMETERS:
1504 : !
1505 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1506 : REAL(hp), INTENT(INOUT) :: outData ! Data w/ dil. factor applied
1507 : !
1508 : ! !OUTPUT PARAMETERS:
1509 : !
1510 : INTEGER, INTENT(OUT) :: RC ! Success or failure
1511 : !
1512 : !
1513 : ! !REMARKS:
1514 : ! This code was abstracted out of Get_Current_Emisssions for clarity. We
1515 : ! have also refactored the code to remove ELSE blocks for better efficiency.
1516 : !EOP
1517 : !------------------------------------------------------------------------------
1518 : !BOC
1519 : !
1520 : ! !LOCAL VARIABLES:
1521 : !
1522 : ! Scalars
1523 : REAL(hp) :: dilFact
1524 :
1525 : ! Strings
1526 : CHARACTER(LEN=255) :: errMsg, thisLoc
1527 :
1528 : !========================================================================
1529 : ! Apply_Dilution_Factor begins here!
1530 : !========================================================================
1531 :
1532 : ! Initialize
1533 0 : RC = HCO_SUCCESS
1534 0 : errMsg = ''
1535 0 : thisLoc = 'Apply_Dilution_Factor (in src/Core/hco_calc_mod.F90)'
1536 0 : dilFact = 0.0_hp
1537 :
1538 : !========================================================================
1539 : ! Get dilution factor. Never dilute 3D emissions.
1540 : !========================================================================
1541 0 : IF ( BaseDct%Dta%SpaceDim == 3 ) THEN
1542 0 : outData = inData
1543 0 : RETURN
1544 : ENDIF
1545 :
1546 : !========================================================================
1547 : ! If emission level mode is 2, copy emissions to all level
1548 : ! A separate scale factor should be used to distribute vertically
1549 : !========================================================================
1550 0 : IF ( BaseDct%Dta%EmisLmode == 2 ) THEN
1551 0 : outData = inData
1552 0 : RETURN
1553 : ENDIF
1554 :
1555 : !========================================================================
1556 : ! Otherwise, compute the vertical dilution factor
1557 : ! and apply it to the input data.
1558 : !========================================================================
1559 : CALL GetDilFact( I = I, &
1560 : J = J, &
1561 : L = L, &
1562 : LowLL = LowLL, &
1563 : UppLL = UppLL, &
1564 : HcoState = HcoState, &
1565 : EmisL1 = BaseDct%Dta%EmisL1, &
1566 : EmisL1Unit = BaseDct%Dta%EmisL1Unit, &
1567 : EmisL2 = BaseDct%Dta%EmisL2, &
1568 : EmisL2Unit = BaseDct%Dta%EmisL2Unit, &
1569 : DilFact = dilFact, &
1570 0 : RC = RC )
1571 :
1572 : ! Trap error
1573 0 : IF ( RC /= HCO_SUCCESS ) THEN
1574 0 : WRITE(*,*) 'Error getting dilution factor at ',I,J,&
1575 0 : ': ', TRIM(BaseDct%cName)
1576 0 : RC = 1
1577 0 : RETURN
1578 : ENDIF
1579 :
1580 : ! Scale base emission by dilution factor
1581 0 : outData = dilFact * inData
1582 :
1583 : END SUBROUTINE Apply_Dilution_Factor
1584 : !EOC
1585 : !------------------------------------------------------------------------------
1586 : ! GEOS-Chem Global Chemical Transport Model !
1587 : !------------------------------------------------------------------------------
1588 : !BOP
1589 : !
1590 : ! !IROUTINE: Apply_Scale_Factor
1591 : !
1592 : ! !DESCRIPTION: Applies scale factors to the input
1593 : !\\
1594 : !\\
1595 : ! !INTERFACE:
1596 : !
1597 0 : SUBROUTINE Apply_Scale_Factor( I, J, L, ScalDct, scalFac, dataVal, RC )
1598 : !
1599 : ! !USES:
1600 : !
1601 : USE HCO_Error_Mod
1602 : !
1603 : ! !INPUT PARAMETERS:
1604 : !
1605 : INTEGER, INTENT(IN) :: I ! Longitude (or X-dim) index
1606 : INTEGER, INTENT(IN) :: J ! Latitude (or Y-dim) index
1607 : INTEGER, INTENT(IN) :: L ! Vertical level index
1608 : TYPE(DataCont), POINTER :: ScalDct ! Scale Factor data container
1609 : REAL(sp), INTENT(IN) :: scalFac ! Scale factor
1610 : !
1611 : ! !INPUT/OUTPUT PARAMETERS:
1612 : !
1613 : REAL(hp), INTENT(INOUT) :: dataVal ! Data to be scaled
1614 : !
1615 : ! !OUTPUT PARAMETERS:
1616 : !
1617 : INTEGER, INTENT(OUT) :: RC ! Success or failure
1618 : !
1619 : ! !REMARKS:
1620 : ! This code was abstracted out of Get_Current_Emisssions for clarity. We
1621 : ! have also refactored the code to remove ELSE blocks for better efficiency.
1622 : !EOP
1623 : !------------------------------------------------------------------------------
1624 : !BOC
1625 : !
1626 : ! !LOCAL VARIABLES:
1627 : !
1628 : ! Strings
1629 : CHARACTER(LEN=255) :: errMsg, thisLoc
1630 :
1631 : !========================================================================
1632 : ! Apply_Scale_Factor begins here!
1633 : !========================================================================
1634 :
1635 : ! Initialize
1636 0 : RC = HCO_SUCCESS
1637 0 : errMsg = ''
1638 0 : thisLoc = 'Apply_Scale_Factor (src/Core/hco_calc_mod.F90)'
1639 :
1640 : !------------------------------------------------------------------------
1641 : ! Operation = 1: multiply
1642 : !------------------------------------------------------------------------
1643 0 : IF ( ScalDct%Oper == 1 ) THEN
1644 0 : dataVal = dataVal * scalFac
1645 0 : RETURN
1646 : ENDIF
1647 :
1648 : !------------------------------------------------------------------------
1649 : ! Operation = -1: divide
1650 : !------------------------------------------------------------------------
1651 0 : IF ( ScalDct%Oper == -1 ) THEN
1652 0 : IF ( scalFac /= 0.0_sp ) THEN
1653 0 : dataVal = dataVal / scalFac
1654 : ENDIF
1655 0 : RETURN
1656 : ENDIF
1657 :
1658 : !------------------------------------------------------------------------
1659 : ! Operation = 2: square
1660 : !------------------------------------------------------------------------
1661 0 : IF ( ScalDct%Oper == 2 ) THEN
1662 0 : dataVal = dataVal * scalFac * scalFac
1663 0 : RETURN
1664 : ENDIF
1665 :
1666 : !------------------------------------------------------------------------
1667 : ! Return w/ error otherwise (Oper 3 is only allowed for masks!)
1668 : !------------------------------------------------------------------------
1669 0 : WRITE( errMsg, * ) 'Illegal operator for: ', TRIM( ScalDct%cName ), &
1670 0 : ' operation: ', ScalDct%Oper
1671 0 : CALL HCO_Error( ErrMsg, RC, thisLoc )
1672 0 : RC = 2
1673 :
1674 : END SUBROUTINE Apply_Scale_Factor
1675 : !EOC
1676 : !------------------------------------------------------------------------------
1677 : ! Harmonized Emissions Component (HEMCO) !
1678 : !------------------------------------------------------------------------------
1679 : !BOP
1680 : !
1681 : ! !IROUTINE: HCO_EvalFld_3D
1682 : !
1683 : ! !DESCRIPTION: Subroutine HCO\_EvalFld\_3D returns the 3D data field belonging
1684 : ! to the emissions list data container with field name 'cName'. The returned
1685 : ! data field is the completely evaluated field, e.g. the base field multiplied
1686 : ! by all scale factors and with all masking being applied (as specified in the
1687 : ! HEMCO configuration file). This distinguished this routine from HCO\_GetPtr
1688 : ! in hco\_emislist\_mod.F90, which returns a reference to the unevaluated data
1689 : ! field.
1690 : !\\
1691 : !\\
1692 : ! !INTERFACE:
1693 : !
1694 0 : SUBROUTINE HCO_EvalFld_3D( HcoState, cName, Arr3D, RC, FOUND )
1695 : !
1696 : ! !USES:
1697 : !
1698 : USE HCO_STATE_MOD, ONLY : HCO_State
1699 : USE HCO_DATACONT_MOD, ONLY : ListCont_Find
1700 : !
1701 : ! !INPUT PARAMETERS:
1702 : !
1703 : CHARACTER(LEN=*), INTENT(IN ) :: cName
1704 : !
1705 : ! !INPUT/OUTPUT PARAMETERS:
1706 : !
1707 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1708 : REAL(hp), INTENT(INOUT) :: Arr3D(:,:,:) ! 3D array
1709 : INTEGER, INTENT(INOUT) :: RC ! Return code
1710 : !
1711 : ! !OUTPUT PARAMETERS:
1712 : !
1713 : LOGICAL, INTENT( OUT), OPTIONAL :: FOUND
1714 : !
1715 : ! !REVISION HISTORY:
1716 : ! 11 May 2015 - C. Keller - Initial Version
1717 : ! See https://github.com/geoschem/hemco for complete history
1718 : !EOP
1719 : !------------------------------------------------------------------------------
1720 : !BOC
1721 : !
1722 : ! !LOCAL VARIABLES:
1723 : !
1724 : ! Scalars
1725 : LOGICAL :: FND
1726 : INTEGER :: AS, nI, nJ, nL, FLAG
1727 :
1728 : ! Arrays
1729 0 : REAL(hp), ALLOCATABLE :: Mask(:,:,:)
1730 :
1731 : ! Working pointers: list and data container
1732 : TYPE(ListCont), POINTER :: Lct
1733 :
1734 : ! For error handling & verbose mode
1735 : CHARACTER(LEN=255) :: MSG
1736 : CHARACTER(LEN=255) :: LOC = "HCO_EvalFld_3d (HCO_calc_mod.F90)"
1737 :
1738 : !=================================================================
1739 : ! HCO_EvalFld_3D begins here!
1740 : !=================================================================
1741 :
1742 : ! Init
1743 0 : RC = HCO_SUCCESS
1744 0 : Lct => NULL()
1745 0 : IF ( PRESENT(FOUND) ) FOUND = .FALSE.
1746 :
1747 : ! Search for base container
1748 0 : CALL ListCont_Find ( HcoState%EmisList, TRIM(cName), FND, Lct )
1749 0 : IF ( PRESENT(FOUND) ) FOUND = FND
1750 :
1751 : ! If not found, return here
1752 0 : IF ( .NOT. FND ) THEN
1753 0 : IF ( PRESENT(FOUND) ) THEN
1754 0 : RETURN
1755 : ELSE
1756 0 : MSG = 'Cannot find in EmisList: ' // TRIM(cName)
1757 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1758 0 : RETURN
1759 : ENDIF
1760 : ENDIF
1761 :
1762 : ! Init
1763 0 : Arr3D = 0.0_hp
1764 :
1765 : ! Define output dimensions
1766 0 : nI = SIZE(Arr3D,1)
1767 0 : nJ = SIZE(Arr3D,2)
1768 0 : nL = SIZE(Arr3D,3)
1769 :
1770 : ! Sanity check: horizontal grid dimensions are expected to be on HEMCO grid
1771 0 : IF ( nI /= HcoState%NX .OR. nJ /= HcoState%nY ) THEN
1772 0 : WRITE(MSG,*) "Horizontal dimension error: ", TRIM(cName), nI, nJ
1773 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1774 0 : RETURN
1775 : ENDIF
1776 :
1777 : ! Make sure mask array is defined
1778 0 : ALLOCATE(MASK(nI,nJ,nL),STAT=AS)
1779 0 : IF ( AS /= 0 ) THEN
1780 0 : CALL HCO_ERROR( 'Cannot allocate MASK', RC, THISLOC=LOC )
1781 0 : RETURN
1782 : ENDIF
1783 0 : mask = 0.0_hp
1784 :
1785 : ! Calculate emissions for base container
1786 : CALL GET_CURRENT_EMISSIONS( HcoState, Lct%Dct, nI, nJ, &
1787 0 : nL, Arr3D, Mask, RC )
1788 0 : IF ( RC /= HCO_SUCCESS ) THEN
1789 0 : CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
1790 0 : RETURN
1791 : ENDIF
1792 :
1793 : ! All done
1794 0 : IF (ALLOCATED(MASK) ) DEALLOCATE(MASK)
1795 0 : Lct => NULL()
1796 :
1797 0 : END SUBROUTINE HCO_EvalFld_3D
1798 : !EOC
1799 : !------------------------------------------------------------------------------
1800 : ! Harmonized Emissions Component (HEMCO) !
1801 : !------------------------------------------------------------------------------
1802 : !BOP
1803 : !
1804 : ! !IROUTINE: HCO_EvalFld_2D
1805 : !
1806 : ! !DESCRIPTION: Subroutine HCO\_EvalFld\_2D returns the 2D data field belonging
1807 : ! to the emissions list data container with field name 'cName'. The returned
1808 : ! data field is the completely evaluated field, e.g. the base field multiplied
1809 : ! by all scale factors and with all masking being applied (as specified in the
1810 : ! HEMCO configuration file). This distinguished this routine from HCO\_GetPtr
1811 : ! in hco\_emislist\_mod.F90, which returns a reference to the unevaluated data
1812 : ! field.
1813 : !\\
1814 : !\\
1815 : !\\
1816 : ! !INTERFACE:
1817 : !
1818 0 : SUBROUTINE HCO_EvalFld_2D( HcoState, cName, Arr2D, RC, FOUND )
1819 : !
1820 : ! !USES:
1821 : !
1822 : USE HCO_STATE_MOD, ONLY : HCO_State
1823 : USE HCO_DATACONT_MOD, ONLY : ListCont_Find
1824 : !
1825 : ! !INPUT PARAMETERS:
1826 : !
1827 : CHARACTER(LEN=*), INTENT(IN ) :: cName
1828 : !
1829 : ! !INPUT/OUTPUT PARAMETERS:
1830 : !
1831 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1832 : REAL(hp), INTENT(INOUT) :: Arr2D(:,:) ! 2D array
1833 : INTEGER, INTENT(INOUT) :: RC ! Return code
1834 : !
1835 : ! !OUTPUT PARAMETERS:
1836 : !
1837 : LOGICAL, INTENT( OUT), OPTIONAL :: FOUND
1838 : !
1839 : ! !REVISION HISTORY:
1840 : ! 11 May 2015 - C. Keller - Initial Version
1841 : ! See https://github.com/geoschem/hemco for complete history
1842 : !EOP
1843 : !------------------------------------------------------------------------------
1844 : !BOC
1845 : !
1846 : ! !LOCAL VARIABLES:
1847 : !
1848 : ! Scalars
1849 : LOGICAL :: FND
1850 : INTEGER :: AS, nI, nJ, nL, UseLL, FLAG
1851 :
1852 : ! Arrays
1853 0 : REAL(hp), ALLOCATABLE :: Mask (:,:,:)
1854 0 : REAL(hp), ALLOCATABLE :: Arr3D(:,:,:)
1855 :
1856 : ! Working pointers: list and data container
1857 : TYPE(ListCont), POINTER :: Lct
1858 :
1859 : ! For error handling & verbose mode
1860 : CHARACTER(LEN=255) :: MSG
1861 : CHARACTER(LEN=255) :: LOC = "HCO_EvalFld_2d (HCO_calc_mod.F90)"
1862 :
1863 : !=================================================================
1864 : ! HCO_EvalFld_2D begins here!
1865 : !=================================================================
1866 :
1867 : ! Init
1868 0 : RC = HCO_SUCCESS
1869 0 : Lct => NULL()
1870 0 : IF ( PRESENT(FOUND) ) FOUND = .FALSE.
1871 :
1872 : ! Search for base container
1873 0 : CALL ListCont_Find ( HcoState%EmisList, TRIM(cName), FND, Lct )
1874 0 : IF ( PRESENT(FOUND) ) FOUND = FND
1875 :
1876 : ! If not found, return here
1877 0 : IF ( .NOT. FND ) THEN
1878 0 : IF ( PRESENT(FOUND) ) THEN
1879 0 : RETURN
1880 : ELSE
1881 0 : MSG = 'Cannot find in EmisList: ' // TRIM(cName)
1882 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1883 0 : RETURN
1884 : ENDIF
1885 : ENDIF
1886 :
1887 : ! Init Arr2D
1888 0 : Arr2D = 0.0_hp
1889 :
1890 : ! Define output dimensions
1891 0 : nI = SIZE(Arr2D,1)
1892 0 : nJ = SIZE(Arr2D,2)
1893 0 : nL = 1
1894 :
1895 : ! Sanity check: horizontal grid dimensions are expected to be on HEMCO grid
1896 0 : IF ( nI /= HcoState%NX .OR. nJ /= HcoState%nY ) THEN
1897 0 : WRITE(MSG,*) "Horizontal dimension error: ", TRIM(cName), nI, nJ
1898 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1899 0 : RETURN
1900 : ENDIF
1901 :
1902 : ! Make sure mask array is defined
1903 0 : ALLOCATE(MASK(nI,nJ,nL),Arr3D(nI,nJ,nL),STAT=AS)
1904 0 : IF ( AS /= 0 ) THEN
1905 0 : CALL HCO_ERROR( 'Cannot allocate MASK', RC, THISLOC=LOC )
1906 0 : RETURN
1907 : ENDIF
1908 0 : Arr3D = 0.0_hp
1909 0 : Mask = 0.0_hp
1910 :
1911 : ! Calculate emissions for base container
1912 : CALL GET_CURRENT_EMISSIONS( HcoState, Lct%Dct, nI, nJ, &
1913 0 : nL, Arr3D, Mask, RC, UseLL=UseLL )
1914 0 : IF ( RC /= HCO_SUCCESS ) THEN
1915 0 : CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
1916 0 : RETURN
1917 : ENDIF
1918 :
1919 : ! Place 3D array into 2D array. UseLL returns the vertical level into which
1920 : ! emissions have been added within GET_CURRENT_EMISSIONS. This should be
1921 : ! level 1 for most cases but it can be another level if specified so.
1922 : ! Return a warning if level is not 1 (ckeller, 11/1/16).
1923 0 : UseLL = MIN( MAX(useLL,1), SIZE(Arr3D,3) )
1924 0 : IF ( UseLL /= 1 ) THEN
1925 0 : WRITE(MSG,*) "2D data was emitted above surface - this information might be lost: " , TRIM(cName), UseLL
1926 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
1927 : ENDIF
1928 :
1929 : ! Pass 3D data to 2D array
1930 0 : Arr2D = Arr3D(:,:,UseLL)
1931 :
1932 : ! All done
1933 0 : IF (ALLOCATED(MASK ) ) DEALLOCATE(MASK )
1934 0 : IF (ALLOCATED(Arr3D) ) DEALLOCATE(Arr3D)
1935 0 : Lct => NULL()
1936 :
1937 0 : END SUBROUTINE HCO_EvalFld_2D
1938 : !EOC
1939 : !------------------------------------------------------------------------------
1940 : ! Harmonized Emissions Component (HEMCO) !
1941 : !------------------------------------------------------------------------------
1942 : !BOP
1943 : !
1944 : ! !IROUTINE: GetMaskVal
1945 : !
1946 : ! !DESCRIPTION: Subroutine GetMaskVal is a helper routine to get the mask
1947 : ! value at a given location.
1948 : !\\
1949 : !\\
1950 : ! !INTERFACE:
1951 : !
1952 0 : SUBROUTINE GetMaskVal ( Dct, I, J, MaskVal, Fractions, RC )
1953 : !
1954 : ! !USES:
1955 : !
1956 : !
1957 : ! !INPUT PARAMETERS:
1958 : !
1959 : INTEGER, INTENT(IN ) :: I ! # of lons
1960 : INTEGER, INTENT(IN ) :: J ! # of lats
1961 : LOGICAL, INTENT(IN ) :: Fractions ! Use fractions?
1962 : !
1963 : ! !INPUT/OUTPUT PARAMETERS:
1964 : !
1965 : TYPE(DataCont), POINTER :: Dct ! Mask container
1966 : REAL(sp), INTENT(INOUT) :: MaskVal
1967 : INTEGER, INTENT(INOUT) :: RC
1968 : !
1969 : ! !REVISION HISTORY:
1970 : ! 09 Apr 2015 - C. Keller - Initial Version
1971 : ! See https://github.com/geoschem/hemco for complete history
1972 : !EOP
1973 : !------------------------------------------------------------------------------
1974 : !BOC
1975 : !
1976 : ! !LOCAL VARIABLES:
1977 : !
1978 :
1979 : !=================================================================
1980 : ! GetMaskVal begins here
1981 : !=================================================================
1982 :
1983 : ! Mask value over this grid box
1984 0 : MaskVal = Dct%Dta%V2(1)%Val(I,J)
1985 :
1986 : ! Negative mask values are treated as zero (exclude).
1987 0 : IF ( (MaskVal <= 0.0_sp) .OR. (MaskVal == HCO_MISSVAL) ) THEN
1988 0 : MaskVal = 0.0_sp
1989 0 : ELSEIF ( MaskVal > 1.0_sp ) THEN
1990 0 : MaskVal = 1.0_sp
1991 : ENDIF
1992 :
1993 : ! For operator set to 3, mirror value
1994 : ! MaskVal=1 becomes 0 and MaskVal=0/missing becomes 1
1995 0 : IF ( Dct%Oper == 3 ) THEN
1996 0 : IF ( (MaskVal == 0.0_sp) .OR. (MaskVal == HCO_MISSVAL) ) THEN
1997 0 : MaskVal = 1.0_sp
1998 0 : ELSEIF ( MaskVal == 1.0_sp ) THEN
1999 0 : MaskVal = 1.0_sp - MaskVal
2000 : ENDIF
2001 : ENDIF
2002 :
2003 : ! Treat as binary?
2004 0 : IF ( .NOT. Fractions ) THEN
2005 0 : IF ( MaskVal < MASK_THRESHOLD ) THEN
2006 0 : MaskVal = 0.0_sp
2007 : ELSE
2008 0 : MaskVal = 1.0_sp
2009 : ENDIF
2010 : ENDIF
2011 :
2012 : ! Return w/ success
2013 0 : RC = HCO_SUCCESS
2014 :
2015 0 : END SUBROUTINE GetMaskVal
2016 : !EOC
2017 : !------------------------------------------------------------------------------
2018 : ! Harmonized Emissions Component (HEMCO) !
2019 : !------------------------------------------------------------------------------
2020 : !BOP
2021 : !
2022 : ! !IROUTINE: HCO_MaskFld
2023 : !
2024 : ! !DESCRIPTION: Subroutine HCO\_MaskFld is a helper routine to get the mask
2025 : ! field with the given name. The returned mask field is fully evaluated,
2026 : ! e.g. the data operation flag associated with this mask field is already
2027 : ! taken into account. For instance, if the data operator of a mask field is
2028 : ! set to 3, the returned array contains already the mirrored mask values.
2029 : !\\
2030 : !\\
2031 : ! !INTERFACE:
2032 : !
2033 0 : SUBROUTINE HCO_MaskFld ( HcoState, MaskName, Mask, RC, FOUND )
2034 : !
2035 : ! !USES:
2036 : !
2037 : USE HCO_STATE_MOD, ONLY : HCO_State
2038 : USE HCO_DATACONT_MOD, ONLY : ListCont_Find
2039 : !
2040 : ! !INPUT PARAMETERS:
2041 : !
2042 : TYPE(HCO_STATE), POINTER :: HcoState
2043 : CHARACTER(LEN=*),INTENT(IN ) :: MaskName
2044 : !
2045 : ! !INPUT/OUTPUT PARAMETERS:
2046 : !
2047 : REAL(sp), INTENT(INOUT) :: Mask(:,:)
2048 : INTEGER, INTENT(INOUT) :: RC
2049 : !
2050 : ! !OUTPUT PARAMETERS:
2051 : !
2052 : LOGICAL, INTENT( OUT), OPTIONAL :: FOUND
2053 : !
2054 : ! !REVISION HISTORY:
2055 : ! 11 Jun 2015 - C. Keller - Initial Version
2056 : ! See https://github.com/geoschem/hemco for complete history
2057 : !EOP
2058 : !------------------------------------------------------------------------------
2059 : !BOC
2060 : !
2061 : ! !LOCAL VARIABLES:
2062 : !
2063 : INTEGER :: I, J, FLAG
2064 :
2065 : LOGICAL :: FND, ERR
2066 : LOGICAL :: Fractions
2067 :
2068 : TYPE(ListCont), POINTER :: MaskLct
2069 :
2070 : CHARACTER(LEN=255) :: MSG
2071 : CHARACTER(LEN=255) :: LOC = 'HCO_MaskFld (hco_calc_mod.F90)'
2072 :
2073 : !=================================================================
2074 : ! HCO_MaskFld begins here
2075 : !=================================================================
2076 :
2077 : ! Nullify
2078 0 : MaskLct => NULL()
2079 :
2080 : ! Init: default is mask value of 1
2081 0 : MASK = 1.0_sp
2082 0 : ERR = .FALSE.
2083 0 : FND = .FALSE.
2084 :
2085 : ! Search for mask field within EmisList
2086 0 : CALL ListCont_Find ( HcoState%EmisList, TRIM(MaskName), FND, MaskLct )
2087 :
2088 0 : IF ( .NOT. FND .AND. .NOT. PRESENT(FOUND) ) THEN
2089 0 : MSG = 'Cannot find mask field ' // TRIM(MaskName)
2090 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1='!')
2091 : MSG = 'Make sure this field is listed in the mask section ' // &
2092 : 'of the HEMCO configuration file. You may also need to ' // &
2093 0 : 'set the optional attribute `ReadAlways` to `yes`, e.g.'
2094 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2095 0 : MSG = '5000 TESTMASK -140/10/-40/90 - - - xy 1 1 -140/10/-40/90 yes'
2096 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
2097 : CALL HCO_ERROR ( &
2098 0 : 'Error reading mask '//TRIM(MaskName), RC, THISLOC=LOC )
2099 0 : RETURN
2100 : ENDIF
2101 0 : IF ( PRESENT(FOUND) ) FOUND = FND
2102 :
2103 : ! Do only if found
2104 0 : IF ( FND ) THEN
2105 :
2106 : ! Use mask fractions?
2107 0 : Fractions = HcoState%Options%MaskFractions
2108 :
2109 : ! Make sure mask array has correct dimensions
2110 0 : IF ( SIZE(MASK,1) /= HcoState%NX .OR. SIZE(MASK,2) /= HcoState%NY ) THEN
2111 0 : WRITE(MSG,*) 'Input mask array has wrong dimensions. Must be ', &
2112 0 : HcoState%NX, HcoState%NY, ' but found ', SIZE(MASK,1), SIZE(MASK,2)
2113 0 : CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
2114 0 : RETURN
2115 : ENDIF
2116 :
2117 : ! Do for every grid box
2118 : !$OMP PARALLEL DO &
2119 : !$OMP DEFAULT( SHARED ) &
2120 : !$OMP PRIVATE( I, J )
2121 0 : DO J = 1, HcoState%NY
2122 0 : DO I = 1, HcoState%NX
2123 0 : CALL GetMaskVal( MaskLct%Dct, I, J, Mask(I,J), Fractions, RC )
2124 0 : IF ( RC /= HCO_SUCCESS ) THEN
2125 : ERR = .TRUE.
2126 : EXIT
2127 : ENDIF
2128 : ENDDO
2129 : ENDDO
2130 : !$OMP END PARALLEL DO
2131 :
2132 : ! Error check
2133 0 : IF ( ERR ) THEN
2134 0 : MSG = 'Error in GetMaskVal'
2135 0 : CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
2136 0 : RETURN
2137 : ENDIF
2138 :
2139 : ENDIF
2140 :
2141 : ! Free pointer
2142 0 : MaskLct => NULL()
2143 :
2144 : ! Return w/ success
2145 0 : RC = HCO_SUCCESS
2146 :
2147 : END SUBROUTINE HCO_MaskFld
2148 : !EOC
2149 : !------------------------------------------------------------------------------
2150 : ! Harmonized Emissions Component (HEMCO) !
2151 : !------------------------------------------------------------------------------
2152 : !BOP
2153 : !
2154 : ! !IROUTINE: GetVertIndx
2155 : !
2156 : ! !DESCRIPTION: Subroutine GetVertIndx is a helper routine to get the vertical
2157 : ! index range of the given data field.
2158 : !\\
2159 : !\\
2160 : ! !INTERFACE:
2161 : !
2162 0 : SUBROUTINE GetVertIndx( HcoState, Dct, isLevDct1, &
2163 : LevDct1, LevDct1_Unit, isLevDct2, &
2164 : LevDct2, LevDct2_Unit, I, &
2165 : J, LowLL, UppLL, RC )
2166 : !
2167 : ! !USES:
2168 : !
2169 : USE HCO_State_Mod, ONLY : HCO_State
2170 : !
2171 : ! !INPUT PARAMETERS:
2172 : !
2173 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2174 : LOGICAL, INTENT(IN) :: isLevDct1 ! Is LevDct1 not null?
2175 : TYPE(DataCont), POINTER :: LevDct1 ! Level index 1 container
2176 : INTEGER, INTENT(IN) :: LevDct1_Unit ! LevDct1 unit code
2177 : LOGICAL, INTENT(IN) :: isLevDct2 ! Is LevDct2 not null?
2178 : TYPE(DataCont), POINTER :: LevDct2 ! Level index 2 container
2179 : INTEGER, INTENT(IN) :: LevDct2_Unit ! LevDct2 unit code
2180 : INTEGER, INTENT(IN) :: I ! lon index
2181 : INTEGER, INTENT(IN) :: J ! lat index
2182 : !
2183 : ! !INPUT/OUTPUT PARAMETERS:
2184 : !
2185 : TYPE(DataCont), POINTER :: Dct ! Mask container
2186 : INTEGER, INTENT(INOUT) :: LowLL ! lower level index
2187 : INTEGER, INTENT(INOUT) :: UppLL ! upper level index
2188 : INTEGER, INTENT(INOUT) :: RC
2189 : !
2190 : ! !REVISION HISTORY:
2191 : ! 06 May 2016 - C. Keller - Initial Version
2192 : ! See https://github.com/geoschem/hemco for complete history
2193 : !EOP
2194 : !------------------------------------------------------------------------------
2195 : !BOC
2196 : !
2197 : ! !LOCAL VARIABLES:
2198 : !
2199 : INTEGER :: EmisLUnit
2200 : REAL(hp) :: EmisL
2201 : CHARACTER(LEN=255) :: errMsg, thisLoc
2202 :
2203 : !=======================================================================
2204 : ! GetVertIndx begins here
2205 : !=======================================================================
2206 :
2207 : ! Initialize
2208 0 : RC = HCO_SUCCESS
2209 0 : errMsg = ''
2210 0 : thisLoc = 'GetVertIndx (src/Core/hco_calc_mod.F90)'
2211 :
2212 : !-----------------------------------------------------------------------
2213 : ! Get vertical extension of base emission array.
2214 : !
2215 : ! Unlike the output array OUTARR_3D, the data containers do not
2216 : ! necessarily extent over the entire troposphere but only cover
2217 : ! the effectively filled vertical levels. For most inventories,
2218 : ! this is only the first model level.
2219 : !-----------------------------------------------------------------------
2220 0 : IF ( Dct%Dta%SpaceDim == 3 ) THEN
2221 0 : LowLL = 1
2222 0 : UppLL = SIZE( Dct%Dta%V3(1)%Val, 3 )
2223 0 : RC = HCO_SUCCESS
2224 0 : RETURN
2225 : ENDIF
2226 :
2227 : !-----------------------------------------------------------------------
2228 : ! For 2D field, check if it shall be spread out over multiple
2229 : ! levels. Possible to go from PBL to max. specified level.
2230 : !-----------------------------------------------------------------------
2231 :
2232 : ! Lower level
2233 : ! --> Check if scale factor is used to determine lower and/or
2234 : ! upper level
2235 : !
2236 : ! NOTE: Get rid of ELSE block for efficiency (bmy, 09 Mar 2023)
2237 0 : EmisL = Dct%Dta%EmisL1
2238 0 : EmisLUnit = Dct%Dta%EmisL1Unit
2239 :
2240 0 : IF ( isLevDct1 ) THEN
2241 0 : EmisL = GetEmisL( HcoState, LevDct1, I, J )
2242 0 : IF ( EmisL < 0.0_hp ) THEN
2243 0 : RC = HCO_FAIL
2244 0 : RETURN
2245 : ENDIF
2246 0 : EmisLUnit = LevDct1_Unit
2247 : ENDIF
2248 :
2249 0 : CALL GetIdx( HcoState, I, J, EmisL, EmisLUnit, LowLL, RC )
2250 0 : IF ( RC /= HCO_SUCCESS ) THEN
2251 0 : errMsg = 'Error encountered in rouitine "GetIdx" (for LevDct1)!'
2252 0 : CALL HCO_ERROR( errMsg, RC, thisLoc )
2253 0 : RETURN
2254 : ENDIF
2255 :
2256 : ! Upper level
2257 : ! NOTE: Get rid of ELSE block for efficiency (bmy, 09 Mar 2023)
2258 0 : EmisL = Dct%Dta%EmisL2
2259 0 : EmisLUnit = Dct%Dta%EmisL2Unit
2260 :
2261 0 : IF ( isLevDct2 ) THEN
2262 0 : EmisL = GetEmisL( HcoState, LevDct2, I, J )
2263 0 : IF ( EmisL < 0.0_hp ) THEN
2264 0 : RC = HCO_FAIL
2265 0 : RETURN
2266 : ENDIF
2267 0 : EmisLUnit = LevDct2_Unit
2268 : ENDIF
2269 :
2270 0 : CALL GetIdx( HcoState, I, J, EmisL, EmisLUnit, UppLL, RC )
2271 0 : IF ( RC /= HCO_SUCCESS ) THEN
2272 0 : errMsg = 'Error encountered in rouitine "GetIdx" (for LevDct2)!'
2273 0 : CALL HCO_ERROR( errMsg, RC, thisLoc )
2274 0 : RETURN
2275 : ENDIF
2276 :
2277 : ! Upper level must not be lower than lower level
2278 0 : UppLL = MAX(LowLL, UppLL)
2279 :
2280 : ! Return w/ success
2281 0 : RC = HCO_SUCCESS
2282 :
2283 : END SUBROUTINE GetVertIndx
2284 : !EOC
2285 : !------------------------------------------------------------------------------
2286 : ! Harmonized Emissions Component (HEMCO) !
2287 : !------------------------------------------------------------------------------
2288 : !BOP
2289 : !
2290 : ! !FUNCTION: GetEmisL
2291 : !
2292 : ! !DESCRIPTION: Returns the emission level read from a scale factor.
2293 : !\\
2294 : !\\
2295 : ! !INTERFACE:
2296 : !
2297 0 : FUNCTION GetEmisL( HcoState, LevDct, I, J ) RESULT ( EmisL )
2298 : !
2299 : ! !USES:
2300 : !
2301 : USE HCO_TYPES_MOD
2302 : USE HCO_STATE_MOD, ONLY : HCO_State
2303 : USE HCO_tIdx_MOD, ONLY : tIDx_GetIndx
2304 : !
2305 : ! !INPUT PARAMETERS:
2306 : !
2307 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2308 : TYPE(DataCont), POINTER :: LevDct ! Level index 1 container
2309 : INTEGER, INTENT(IN ) :: I, J ! horizontal index
2310 : !
2311 : ! !RETURN VALUE:
2312 : !
2313 : REAL(hp) :: EmisL
2314 : !
2315 : ! !REVISION HISTORY:
2316 : ! 26 Jan 2018 - C. Keller - Initial version
2317 : ! See https://github.com/geoschem/hemco for complete history
2318 : !EOP
2319 : !------------------------------------------------------------------------------
2320 : !BOC
2321 : !
2322 : ! !LOCAL VARIABLES:
2323 : !
2324 : INTEGER :: levtidx
2325 :
2326 : !=================================================================
2327 : ! GetEmisL begins here
2328 : !=================================================================
2329 0 : levtidx = tIDx_GetIndx( HcoState, LevDct%Dta, I, J )
2330 0 : IF ( levtidx <= 0 ) THEN
2331 : WRITE(*,*)' Cannot get time slice for field '//&
2332 0 : TRIM(LevDct%cName)//': GetEmisL (hco_calc_mod.F90)'
2333 0 : EmisL = -1.0
2334 0 : RETURN
2335 : ENDIF
2336 :
2337 0 : IF ( LevDct%Dta%SpaceDim == 1 ) THEN
2338 0 : EmisL = LevDct%Dta%V2(levtidx)%Val(1,1)
2339 0 : ELSEIF ( LevDct%Dta%SpaceDim == 2 ) THEN
2340 0 : EmisL = LevDct%Dta%V2(levtidx)%Val(I,J)
2341 0 : ELSEIF ( LevDct%Dta%SpaceDim == 3 ) THEN
2342 0 : EmisL = LevDct%Dta%V3(levtidx)%Val(I,J,1)
2343 : ENDIF
2344 :
2345 0 : IF ( EmisL == HCO_MISSVAL ) EmisL = 0.0_hp
2346 :
2347 : END FUNCTION GetEmisL
2348 : !EOC
2349 : !------------------------------------------------------------------------------
2350 : ! Harmonized Emissions Component (HEMCO) !
2351 : !------------------------------------------------------------------------------
2352 : !BOP
2353 : !
2354 : ! !FUNCTION: GetEmisLUnit
2355 : !
2356 : ! !DESCRIPTION: Returns the emission level unit read from a scale factor.
2357 : !\\
2358 : !\\
2359 : ! !INTERFACE:
2360 : !
2361 0 : FUNCTION GetEmisLUnit( HcoState, LevDct ) RESULT( EmisLUnit )
2362 : !
2363 : ! !USES:
2364 : !
2365 : USE HCO_TYPES_MOD
2366 : USE HCO_STATE_MOD, ONLY : HCO_State
2367 : !
2368 : ! !INPUT PARAMETERS:
2369 : !
2370 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2371 : TYPE(DataCont), POINTER :: LevDct ! Level index 1 container
2372 : !
2373 : ! !RETURN VALUE:
2374 : !
2375 : INTEGER :: EmisLUnit
2376 : !
2377 : ! !REVISION HISTORY:
2378 : ! 26 Jan 2018 - C. Keller - Initial version
2379 : ! See https://github.com/geoschem/hemco for complete history
2380 : !EOP
2381 : !------------------------------------------------------------------------------
2382 : !BOC
2383 : !
2384 : ! !LOCAL VARIABLES:
2385 : !
2386 : !=================================================================
2387 : ! GetEmisLUnit begins here
2388 : !=================================================================
2389 :
2390 : ! For now, only meters are supported
2391 0 : EmisLUnit = HCO_EMISL_M
2392 :
2393 : ! Dummy check that units on field are actually in meters
2394 0 : IF ( TRIM(LevDct%Dta%OrigUnit) /= 'm' .AND. &
2395 : TRIM(LevDct%Dta%OrigUnit) /= '1' ) THEN
2396 : WRITE(*,*) TRIM(LevDct%cName)// &
2397 : ' must have units of `m`, instead found '//&
2398 0 : TRIM(LevDct%Dta%OrigUnit)//': GetEmisLUnit (hco_calc_mod.F90)'
2399 0 : EmisLUnit = -1
2400 : ENDIF
2401 :
2402 0 : END FUNCTION GetEmisLUnit
2403 : !EOC
2404 : !------------------------------------------------------------------------------
2405 : ! Harmonized Emissions Component (HEMCO) !
2406 : !------------------------------------------------------------------------------
2407 : !BOP
2408 : !
2409 : ! !IROUTINE: GetIdx
2410 : !
2411 : ! !DESCRIPTION: Subroutine GetIdx is a helper routine to return the vertical
2412 : ! level index for a given altitude. The altitude can be provided in level
2413 : ! coordinates, in units of meters or as the 'PBL mixing height'.
2414 : !\\
2415 : !\\
2416 : ! !INTERFACE:
2417 : !
2418 0 : SUBROUTINE GetIdx( HcoState, I, J, alt, altu, lidx, RC )
2419 : !
2420 : ! !USES:
2421 : !
2422 : USE HCO_TYPES_MOD
2423 : USE HCO_STATE_MOD, ONLY : HCO_STATE
2424 : !
2425 : ! !INPUT PARAMETERS:
2426 : !
2427 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2428 : INTEGER, INTENT(IN ) :: I, J ! horizontal index
2429 : INTEGER, INTENT(IN ) :: altu ! altitude unit
2430 : !
2431 : ! !OUTPUT PARAMETERS:
2432 : !
2433 : INTEGER, INTENT( OUT) :: lidx ! level index
2434 : !
2435 : ! !INPUT/OUTPUT PARAMETERS:
2436 : !
2437 : REAL(hp), INTENT(INOUT) :: alt ! altitude
2438 : INTEGER, INTENT(INOUT) :: RC
2439 : !
2440 : ! !REMARKS:
2441 : ! The code was refactored to avoid ELSE blocks (which are computationally
2442 : ! expensive) following the "never-nesting" technique. (Bob Y., 10 Mar 2023)
2443 : !
2444 : ! !REVISION HISTORY:
2445 : ! 09 May 2016 - C. Keller - Initial version
2446 : ! See https://github.com/geoschem/hemco for complete history
2447 : !EOP
2448 : !------------------------------------------------------------------------------
2449 : !BOC
2450 : !
2451 : ! !LOCAL VARIABLES:
2452 : !
2453 : INTEGER :: L
2454 : REAL(hp) :: altb, altt
2455 : CHARACTER(LEN=255) :: MSG
2456 : CHARACTER(LEN=255) :: LOC = 'GetIdx (hco_calc_mod.F90)'
2457 :
2458 : !=================================================================
2459 : ! HCO_GetIdx begins here
2460 : !=================================================================
2461 :
2462 : ! Init
2463 0 : RC = HCO_SUCCESS
2464 :
2465 : ! Input data is in level coordinates;
2466 : ! Return the corresponding level index
2467 0 : IF ( altu == HCO_EMISL_LEV ) THEN
2468 0 : lidx = INT( alt )
2469 0 : RETURN
2470 : ENDIF
2471 :
2472 : ! Input specifies the top-of-atmosphere;
2473 : ! Return the top-of-atmosphere level index
2474 0 : IF ( altu == HCO_EMISL_TOP ) THEN
2475 0 : lidx = HCOState%NZ
2476 0 : RETURN
2477 : ENDIF
2478 :
2479 : ! Input data is in meters or specifies the PBL top;
2480 : ! Find the corresponding level index
2481 0 : IF ( altu == HCO_EMISL_M .OR. altu == HCO_EMISL_PBL ) THEN
2482 :
2483 : ! Eventually get altitude from PBL height
2484 0 : IF ( altu == HCO_EMISL_PBL ) THEN
2485 0 : alt = HcoState%Grid%PBLHEIGHT%Val(I,J)
2486 : ENDIF
2487 :
2488 : ! Special case of negative height
2489 0 : IF ( alt <= 0.0_hp ) THEN
2490 0 : lidx = 1
2491 0 : RETURN
2492 : ENDIF
2493 :
2494 : ! Loop over data until we are within desired level
2495 : ! NOTE: This can be rewritten more efficiently (bmy, 3/5/21)
2496 0 : altt = 0.0_hp
2497 0 : altb = 0.0_hp
2498 0 : lidx = -1
2499 0 : DO L = 1, HcoState%NZ
2500 0 : altt = altb + HcoState%Grid%BXHEIGHT_M%Val(I,J,L)
2501 0 : IF ( alt >= altb .AND. alt < altt ) THEN
2502 0 : lidx = L
2503 0 : RETURN
2504 : ENDIF
2505 0 : altb = altt
2506 : ENDDO
2507 :
2508 : ! If altitude is above maximum level
2509 0 : IF ( lidx == -1 .AND. alt >= altt ) THEN
2510 0 : lidx = HcoState%NZ
2511 : WRITE(MSG,*) &
2512 0 : 'Level is above max. grid box level - use top level ', alt
2513 0 : CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2514 0 : RETURN
2515 : ENDIF
2516 :
2517 : RETURN
2518 : ENDIF
2519 :
2520 : ! Error if we drop down to here
2521 0 : MSG = 'Illegal altitude unit'
2522 0 : CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
2523 :
2524 : END SUBROUTINE GetIdx
2525 : !EOC
2526 : !------------------------------------------------------------------------------
2527 : ! Harmonized Emissions Component (HEMCO) !
2528 : !------------------------------------------------------------------------------
2529 : !BOP
2530 : !
2531 : ! !IROUTINE: GetDilFact
2532 : !
2533 : ! !DESCRIPTION: Subroutine GetDilFact returns the vertical dilution factor,
2534 : ! that is the factor that is to be applied to distribute emissions into
2535 : ! multiple vertical levels. If grid box height information are available,
2536 : ! these are used to compute the distribution factor. Otherwise, equal weight
2537 : ! is given to all vertical levels.
2538 : !\\
2539 : !\\
2540 : ! !TODO: Dilution factors are currently only weighted by grid box heights
2541 : ! (if these information are available) but any pressure information are
2542 : ! ignored.
2543 : !\\
2544 : !\\
2545 : ! !INTERFACE:
2546 : !
2547 0 : SUBROUTINE GetDilFact( HcoState, EmisL1, EmisL1Unit, EmisL2, &
2548 : EmisL2Unit, I, J, L, &
2549 : LowLL, UppLL, DilFact, RC )
2550 : !
2551 : ! !USES:
2552 : !
2553 : USE HCO_STATE_MOD, ONLY : HCO_State
2554 : !
2555 : ! !INPUT PARAMETERS:
2556 : !
2557 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2558 : INTEGER, INTENT(IN) :: I ! lon index
2559 : INTEGER, INTENT(IN) :: J ! lat index
2560 : INTEGER, INTENT(IN) :: L ! lev index
2561 : INTEGER, INTENT(IN) :: LowLL ! lower level index
2562 : INTEGER, INTENT(IN) :: UppLL ! upper level index
2563 : !
2564 : ! !OUTPUT PARAMETERS:
2565 : !
2566 : REAL(hp), INTENT(OUT) :: DilFact ! Dilution factor
2567 : !
2568 : ! !INPUT/OUTPUT PARAMETERS:
2569 : !
2570 : REAL(hp), INTENT(INOUT) :: EmisL1
2571 : INTEGER, INTENT(INOUT) :: EmisL1Unit
2572 : REAL(hp), INTENT(INOUT) :: EmisL2
2573 : INTEGER, INTENT(INOUT) :: EmisL2Unit
2574 : INTEGER, INTENT(INOUT) :: RC
2575 : !
2576 : ! !REVISION HISTORY:
2577 : ! 06 May 2016 - C. Keller - Initial Version
2578 : ! See https://github.com/geoschem/hemco for complete history
2579 : !EOP
2580 : !------------------------------------------------------------------------------
2581 : !BOC
2582 : !
2583 : ! !LOCAL VARIABLES:
2584 : !
2585 : INTEGER :: L1
2586 : CHARACTER(LEN=255) :: MSG
2587 : CHARACTER(LEN=255) :: LOC = 'GetDilFact (hco_calc_mod.F90)'
2588 : REAL(hp) :: h1, h2, dh, dh1, dh2
2589 : REAL(hp) :: UppLLR, LowLLR
2590 :
2591 : !=================================================================
2592 : ! GetDilFact begins here
2593 : !=================================================================
2594 :
2595 : ! Init
2596 0 : DilFact = 1.0_hp
2597 0 : RC = HCO_SUCCESS
2598 :
2599 : ! Nothing to do if it's only one level
2600 0 : IF ( LowLL == UppLL ) RETURN
2601 :
2602 : ! Compute dilution factor based on boxheights if this information
2603 : ! is available
2604 0 : IF ( ASSOCIATED( HcoState%Grid%BXHEIGHT_M%Val ) ) THEN
2605 :
2606 : ! Get height of bottom level LowLL (in m)
2607 0 : IF ( EmisL1Unit == HCO_EMISL_M ) THEN
2608 0 : h1 = EmisL1
2609 0 : ELSEIF ( EmisL1Unit == HCO_EMISL_PBL ) THEN
2610 0 : h1 = HcoState%Grid%PBLHEIGHT%Val(I,J)
2611 : ELSE
2612 0 : IF ( LowLL > 1 ) THEN
2613 0 : h1 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:(LowLL-1)))
2614 : ELSE
2615 : h1 = 0.0_hp
2616 : ENDIF
2617 : ENDIF
2618 :
2619 : ! Get height of top level UppLL (in m)
2620 0 : IF ( EmisL2Unit == HCO_EMISL_M ) THEN
2621 0 : h2 = EmisL2
2622 0 : ELSEIF ( EmisL2Unit == HCO_EMISL_PBL ) THEN
2623 0 : h2 = HcoState%Grid%PBLHEIGHT%Val(I,J)
2624 : ELSE
2625 0 : h2 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:UppLL))
2626 : ENDIF
2627 :
2628 : ! If vertical weight option is enabled, calculate vertical
2629 : ! distribution factor relative to the grid cell heights. This
2630 : ! is the default (and recommended) option as this makes sure
2631 : ! that the same amount of mass is emitted into each layer.
2632 0 : IF ( HcoState%Options%VertWeight ) THEN
2633 :
2634 : ! Height of grid box of interest (in m)
2635 0 : dh = HcoState%Grid%BXHEIGHT_M%Val(I,J,L)
2636 :
2637 : ! Adjust dh if we are in lowest level
2638 0 : IF ( L == LowLL ) THEN
2639 0 : dh = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:LowLL)) - h1
2640 : ENDIF
2641 :
2642 : ! Adjust dh if we are in top level
2643 0 : IF ( L == UppLL ) THEN
2644 0 : dh = h2 - SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:(UppLL-1)))
2645 : ENDIF
2646 :
2647 : ! compute dilution factor: the new flux should emit the same mass per
2648 : ! volume, i.e. flux_total/column_total = flux_level/column_level
2649 : ! --> flux_level = fluxtotal * column_level / column_total.
2650 0 : IF ( h2 > h1 ) THEN
2651 0 : DilFact = dh / ( h2 - h1 )
2652 : ELSE
2653 0 : MSG = 'GetDilFact h2 not greater than h1'
2654 0 : CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2655 0 : RETURN
2656 : ENDIF
2657 :
2658 : ! If VertWeight option is turned off, emit the same flux in each layer.
2659 : ! Since model layers have different depths, this will result in differnt
2660 : ! total emissions per layer.
2661 : ELSE
2662 :
2663 : ! Get fractional layer indeces for lower and upper level. This makes
2664 : ! sure that only fractions of the lower and upper level are being
2665 : ! considered, so that double-counting is avoided if a model layer
2666 : ! serves both as the top layer and the bottom layer (e.g., wildfire
2667 : ! emissions emitted from bottom to the top of PBL, and from the top
2668 : ! of PBL to 5000m).
2669 0 : LowLLR = REAL(LowLL,hp) - 1.0_hp
2670 0 : UppLLR = REAL(UppLL,hp)
2671 0 : dh1 = 0.0_hp
2672 0 : DO L1 = 1, HcoState%NZ
2673 0 : dh2 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:L1))
2674 0 : IF ( h1 >= dh1 .AND. h1 < dh2 ) THEN
2675 0 : LowLLR = REAL(L1,hp) - ( (dh2-h1)/(dh2-dh1) )
2676 : ENDIF
2677 0 : IF ( h2 > dh1 .AND. h2 <= dh2 ) THEN
2678 0 : UppLLR = REAL(L1,hp) - ( (dh2-h2)/(dh2-dh1) )
2679 : ENDIF
2680 : ! top layer is bottom layer in next loop
2681 0 : dh1 = dh2
2682 : ENDDO
2683 :
2684 : ! Dilution factor using fractional levels
2685 0 : IF ( UppLLR <= LowLLR ) THEN
2686 0 : DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp)
2687 : ELSE
2688 0 : DilFact = 1.0_hp / (UppLLR-LowLLR)
2689 : ENDIF
2690 :
2691 : ENDIF
2692 :
2693 : ! Approximate dilution factor otherwise
2694 : ELSE
2695 :
2696 0 : DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp)
2697 : ENDIF
2698 :
2699 : ! Return w/ success
2700 0 : RC = HCO_SUCCESS
2701 :
2702 : END SUBROUTINE GetDilFact
2703 : #ifdef ADJOINT
2704 : !BOP
2705 : !
2706 : ! !IROUTINE: Get_Current_Emissions
2707 : !
2708 : ! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current
2709 : ! emissions for the specified emission container.
2710 : ! This subroutine is only called by HCO\_CalcEmis and for base emission
2711 : ! containers, i.e. containers of type 1.
2712 : !\\
2713 : !\\
2714 : ! !INTERFACE:
2715 : !
2716 : SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct, &
2717 : nI, nJ, nL, OUTARR_3D, MASK, RC, UseLL )
2718 : !
2719 : ! !USES:
2720 : !
2721 : USE HCO_State_Mod, ONLY : HCO_State
2722 : USE HCO_tIdx_MOD, ONLY : tIDx_GetIndx
2723 : USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined
2724 : !
2725 : ! !INPUT PARAMETERS:
2726 : !
2727 : INTEGER, INTENT(IN) :: nI ! # of lons
2728 : INTEGER, INTENT(IN) :: nJ ! # of lats
2729 : INTEGER, INTENT(IN) :: nL ! # of levs
2730 : !
2731 : ! !INPUT/OUTPUT PARAMETERS:
2732 : !
2733 :
2734 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
2735 : TYPE(DataCont), POINTER :: BaseDct ! base emission
2736 : ! container
2737 : REAL(hp), INTENT(INOUT) :: OUTARR_3D(nI,nJ,nL) ! output array
2738 : REAL(hp), INTENT(INOUT) :: MASK (nI,nJ,nL) ! mask array
2739 : INTEGER, INTENT(INOUT) :: RC
2740 : !
2741 : ! !OUTPUT PARAMETERS:
2742 : !
2743 : INTEGER, INTENT( OUT), OPTIONAL :: UseLL
2744 : !
2745 : ! !REMARKS:
2746 : ! This routine uses multiple loops over all grid boxes (base emissions
2747 : ! and scale factors use separate loops). In an OMP environment, this approach
2748 : ! seems to be faster than using only one single loop (but repeated calls to
2749 : ! point to containers, etc.). The alternative approach is used in routine
2750 : ! Get\_Current\_Emissions\_B at the end of this module and may be employed
2751 : ! on request.
2752 : !
2753 : ! !REVISION HISTORY:
2754 : ! 25 Aug 2012 - C. Keller - Initial Version
2755 : ! 09 Nov 2012 - C. Keller - MASK update. Masks are now treated
2756 : ! separately so that multiple masks can be
2757 : ! added.
2758 : ! 06 Jun 2014 - R. Yantosca - Cosmetic changes in ProTeX headers
2759 : ! 07 Sep 2014 - C. Keller - Mask update. Now set mask to zero as soon as
2760 : ! on of the applied masks is zero.
2761 : ! 03 Dec 2014 - C. Keller - Now calculate time slice index on-the-fly.
2762 : ! 29 Dec 2014 - C. Keller - Added scale factor masks.
2763 : ! 02 Mar 2015 - C. Keller - Now check for missing values. Missing values are
2764 : ! excluded from emission calculation.
2765 : ! 26 Oct 2016 - R. Yantosca - Don't nullify local ptrs in declaration stmts
2766 : ! 11 May 2017 - C. Keller - Added universal scaling
2767 : !EOP
2768 : !------------------------------------------------------------------------------
2769 : !BOC
2770 : !
2771 : ! !LOCAL VARIABLES:
2772 : !
2773 : ! Pointers
2774 : TYPE(DataCont), POINTER :: ScalDct
2775 : TYPE(DataCont), POINTER :: MaskDct
2776 : TYPE(DataCont), POINTER :: LevDct1
2777 : TYPE(DataCont), POINTER :: LevDct2
2778 :
2779 : ! Scalars
2780 : REAL(sp) :: TMPVAL, MaskScale
2781 : REAL(hp) :: DilFact
2782 : REAL(hp) :: ScalFact
2783 : INTEGER :: tIDx, IDX
2784 : INTEGER :: I, J, L, N
2785 : INTEGER :: LowLL, UppLL, ScalLL, TmpLL
2786 : INTEGER :: ERROR
2787 : INTEGER :: TotLL, nnLL
2788 : CHARACTER(LEN=255) :: MSG, LOC
2789 : LOGICAL :: NegScalExist
2790 : LOGICAL :: MaskFractions
2791 : LOGICAL :: isLevDct1
2792 : LOGICAL :: isLevDct2
2793 : LOGICAL :: isMaskDct
2794 : LOGICAL :: isPblHt
2795 : LOGICAL :: isBoxHt
2796 : INTEGER :: LevDct1_Unit
2797 : INTEGER :: LevDct2_Unit
2798 :
2799 : ! testing only
2800 : INTEGER :: IX, IY
2801 :
2802 : !=================================================================
2803 : ! GET_CURRENT_EMISSIONS begins here
2804 : !=================================================================
2805 :
2806 : ! Initialize
2807 : ScalDct => NULL()
2808 : MaskDct => NULL()
2809 : LOC = 'GET_CURRENT_EMISSIONS_ADJ (hco_calc_mod.F90)'
2810 :
2811 : ! Enter
2812 : CALL HCO_ENTER(HcoState%Config%Err, LOC, RC )
2813 : IF(RC /= HCO_SUCCESS) RETURN
2814 :
2815 : ! testing only:
2816 : IX = 3 !-1
2817 : IY = 8 !-1
2818 :
2819 : ! Check if container contains data
2820 : IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN
2821 : MSG = 'Array not defined: ' // TRIM(BaseDct%cName)
2822 : CALL HCO_ERROR( MSG, RC )
2823 : RETURN
2824 : ENDIF
2825 :
2826 : ! Initialize mask. By default, assume that we use all grid boxes.
2827 : MASK(:,:,:) = 1.0_hp
2828 : MaskFractions = HcoState%Options%MaskFractions
2829 :
2830 : ! Verbose
2831 : IF ( HCO_IsVerb(HcoState%Config%Err) ) THEN
2832 : WRITE(MSG,*) 'Evaluate field ', TRIM(BaseDct%cName)
2833 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
2834 : ENDIF
2835 :
2836 : ! ----------------------------------------------------------------
2837 : ! Set base emissions
2838 : ! ----------------------------------------------------------------
2839 :
2840 : ! Initialize ERROR. Will be set to 1 if error occurs below
2841 : ERROR = 0
2842 :
2843 : ! Initialize variables to compute average vertical level index
2844 : totLL = 0
2845 : nnLL = 0
2846 :
2847 : ! Check for level index containers
2848 : IF ( BaseDct%levScalID1 > 0 ) THEN
2849 : CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC )
2850 : IF ( RC /= HCO_SUCCESS ) THEN
2851 : CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
2852 : RETURN
2853 : ENDIF
2854 : ELSE
2855 : LevDct1 => NULL()
2856 : ENDIF
2857 : IF ( BaseDct%levScalID2 > 0 ) THEN
2858 : CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC )
2859 : IF ( RC /= HCO_SUCCESS ) THEN
2860 : CALL HCO_ERROR( 'ERROR 19', RC, THISLOC=LOC )
2861 : RETURN
2862 : ENDIF
2863 : ELSE
2864 : LevDct2 => NULL()
2865 : ENDIF
2866 :
2867 : ! Test whether LevDct1 and LevDct2 are associated
2868 : isLevDct1 = ASSOCIATED( LevDct1 )
2869 : isLevDct2 = ASSOCIATED( LevDct2 )
2870 :
2871 : ! Get the units of LevDct1 (if it exists)
2872 : IF ( isLevDct1 ) THEN
2873 : LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
2874 : IF ( LevDct1_Unit < 0 ) THEN
2875 : MSG = 'LevDct1 units are not defined!'
2876 : CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2877 : RC = HCO_FAIL
2878 : RETURN
2879 : ENDIF
2880 : ELSE
2881 : LevDct1_Unit = -1
2882 : ENDIF
2883 :
2884 : ! Get the units of LevDct2 (if it exists)
2885 : IF ( isLevDct2 ) THEN
2886 : LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
2887 : IF ( LevDct2_Unit < 0 ) THEN
2888 : MSG = 'LevDct2_Units are not defined!'
2889 : CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2890 : RETURN
2891 : ENDIF
2892 : ELSE
2893 : LevDct2_Unit = -1
2894 : ENDIF
2895 :
2896 : ! Throw an error if boxheight is missing and the units are in meters
2897 : IF ( LevDct1_Unit == HCO_EMISL_M .or. &
2898 : LevDct2_Unit == HCO_EMISL_M ) THEN
2899 : IF ( .NOT. ASSOCIATED(HcoState%Grid%BXHEIGHT_M%Val) ) THEN
2900 : MSG = 'Boxheight (in meters) is missing in HEMCO state'
2901 : CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2902 : RETURN
2903 : ENDIF
2904 : ENDIF
2905 :
2906 : ! Throw an error if boxheight is missing and the units are in PBL frac
2907 : IF ( LevDct1_Unit == HCO_EMISL_PBL .or. &
2908 : LevDct2_Unit == HCO_EMISL_PBL ) THEN
2909 : IF ( .NOT. ASSOCIATED(HcoState%Grid%PBLHEIGHT%Val) ) THEN
2910 : MSG = 'Boundary layer height is missing in HEMCO state'
2911 : CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
2912 : RETURN
2913 : ENDIF
2914 : ENDIF
2915 :
2916 : ! Loop over all latitudes and longitudes
2917 : !$OMP PARALLEL DO &
2918 : !$OMP DEFAULT( SHARED )&
2919 : !$OMP PRIVATE( I, J, L, tIdx, TMPVAL, DilFact, LowLL, UppLL )&
2920 : !$OMP COLLAPSE( 2 )&
2921 : !$OMP SCHEDULE( DYNAMIC, 4 )&
2922 : !$OMP REDUCTION( +:totLL )&
2923 : !$OMP REDUCTION( +:nnLL )
2924 : DO J = 1, nJ
2925 : DO I = 1, nI
2926 :
2927 : ! Zero for safety's sake
2928 : totLL = 0
2929 : nnLL = 0
2930 :
2931 : ! Get current time index for this container and at this location
2932 : tIDx = tIDx_GetIndx( HcoState, BaseDct%Dta, I, J )
2933 : IF ( tIDx < 1 ) THEN
2934 : WRITE(MSG,*) 'Cannot get time slice index at location ',I,J,&
2935 : ': ', TRIM(BaseDct%cName), tIDx
2936 : ERROR = 1
2937 : EXIT
2938 : ENDIF
2939 :
2940 : ! Get lower and upper vertical index
2941 : CALL GetVertIndx ( HcoState, BaseDct, isLevDct1, LevDct1, &
2942 : LevDct1_Unit, isLevDct2, LevDct2, LevDct2_Unit, &
2943 : I, J, LowLL, UppLL, &
2944 : RC )
2945 : IF ( RC /= HCO_SUCCESS ) THEN
2946 : WRITE(MSG,*) 'Error getting vertical index at location ',I,J,&
2947 : ': ', TRIM(BaseDct%cName)
2948 : ERROR = 1 ! Will cause error
2949 : EXIT
2950 : ENDIF
2951 :
2952 : ! average upper level
2953 : totLL = totLL + UppLL
2954 : nnLL = nnLL + 1
2955 :
2956 : ! Loop over all levels
2957 : DO L = LowLL, UppLL
2958 :
2959 : ! Get base value. Use uniform value if scalar field.
2960 : IF ( BaseDct%Dta%SpaceDim == 1 ) THEN
2961 : TMPVAL = BaseDct%Dta%V2(tIDx)%Val(1,1)
2962 : ELSEIF ( BaseDct%Dta%SpaceDim == 2 ) THEN
2963 : TMPVAL = BaseDct%Dta%V2(tIDx)%Val(I,J)
2964 : ELSE
2965 : TMPVAL = BaseDct%Dta%V3(tIDx)%Val(I,J,L)
2966 : ENDIF
2967 :
2968 : ! If it's a missing value, mask box as unused and set value to zero
2969 : IF ( TMPVAL == HCO_MISSVAL ) THEN
2970 : MASK(I,J,:) = 0.0_hp
2971 : OUTARR_3D(I,J,L) = 0.0_hp
2972 :
2973 : ! Pass base value to output array
2974 : ELSE
2975 :
2976 : ! Get dilution factor. Never dilute 3D emissions.
2977 : IF ( BaseDct%Dta%SpaceDim == 3 ) THEN
2978 : DilFact = 1.0_hp !1.0
2979 :
2980 : ! 2D dilution factor
2981 : ELSE
2982 : CALL GetDilFact ( HcoState, BaseDct%Dta%EmisL1, &
2983 : BaseDct%Dta%EmisL1Unit, BaseDct%Dta%EmisL2, &
2984 : BaseDct%Dta%EmisL2Unit, I, J, L, LowLL, &
2985 : UppLL, DilFact, RC )
2986 : IF ( RC /= HCO_SUCCESS ) THEN
2987 : WRITE(MSG,*) 'Error getting dilution factor at ',I,J,&
2988 : ': ', TRIM(BaseDct%cName)
2989 : ERROR = 1
2990 : EXIT
2991 : ENDIF
2992 : ENDIF
2993 :
2994 : ! Scale base emission by dilution factor
2995 : OUTARR_3D(I,J,L) = DilFact * TMPVAL
2996 : ENDIF
2997 : ENDDO !L
2998 :
2999 : ENDDO !I
3000 : ENDDO !J
3001 : !$OMP END PARALLEL DO
3002 :
3003 : ! Check for error
3004 : IF ( ERROR == 1 ) THEN
3005 : CALL HCO_ERROR( MSG, RC )
3006 : RETURN
3007 : ENDIF
3008 :
3009 : ! ----------------------------------------------------------------
3010 : ! Apply scale factors
3011 : ! The container IDs of all scale factors associated with this base
3012 : ! container are stored in vector Scal_cID.
3013 : ! ----------------------------------------------------------------
3014 :
3015 : ! Loop over scale factors
3016 : IF ( BaseDct%nScalID > 0 ) THEN
3017 :
3018 : DO N = 1, BaseDct%nScalID
3019 :
3020 : ! Get the scale factor container ID for the current slot
3021 : IDX = BaseDct%Scal_cID(N)
3022 :
3023 : ! Point to data container with the given container ID
3024 : CALL Pnt2DataCont( HcoState, IDX, ScalDct, RC )
3025 : IF ( RC /= HCO_SUCCESS ) THEN
3026 : CALL HCO_ERROR( 'ERROR 20', RC, THISLOC=LOC )
3027 : RETURN
3028 : ENDIF
3029 :
3030 : ! Sanity check: scale field cannot be a base field
3031 : IF ( (ScalDct%DctType == HCO_DCTTYPE_BASE) ) THEN
3032 : MSG = 'Wrong scale field type: ' // TRIM(ScalDct%cName)
3033 : CALL HCO_ERROR( MSG, RC )
3034 : RETURN
3035 : ENDIF
3036 :
3037 : ! Skip this scale factor if no data defined. This is possible
3038 : ! if scale factors are only defined for a given time range and
3039 : ! the simulation datetime is outside of this range.
3040 : IF ( .NOT. FileData_ArrIsDefined(ScalDct%Dta) ) THEN
3041 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
3042 : MSG = 'Skip scale factor '//TRIM(ScalDct%cName)// &
3043 : ' because it is not defined for this datetime.'
3044 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3045 : ENDIF
3046 : CYCLE
3047 : ENDIF
3048 :
3049 : ! Verbose mode
3050 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
3051 : MSG = 'Applying scale factor ' // TRIM(ScalDct%cName)
3052 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3053 : ENDIF
3054 :
3055 : ! Get vertical extension of this scale factor array.
3056 : IF( (ScalDct%Dta%SpaceDim<=2) ) THEN
3057 : ScalLL = 1
3058 : ELSE
3059 : ScalLL = SIZE(ScalDct%Dta%V3(1)%Val,3)
3060 : ENDIF
3061 :
3062 : ! Check if there is a mask field associated with this scale
3063 : ! factor. In this case, get a pointer to the corresponding
3064 : ! mask field and evaluate scale factors only inside the mask
3065 : ! region.
3066 : IF ( ASSOCIATED(ScalDct%Scal_cID) ) THEN
3067 : CALL Pnt2DataCont( HcoState, ScalDct%Scal_cID(1), MaskDct, RC )
3068 : IF ( RC /= HCO_SUCCESS ) THEN
3069 : CALL HCO_ERROR( 'ERROR 21', RC, THISLOC=LOC )
3070 : RETURN
3071 : ENDIF
3072 :
3073 : ! Must be mask field
3074 : IF ( MaskDct%DctType /= HCO_DCTTYPE_MASK ) THEN
3075 : MSG = 'Invalid mask for scale factor: '//TRIM(ScalDct%cName)
3076 : MSG = TRIM(MSG) // '; mask: '//TRIM(MaskDct%cName)
3077 : CALL HCO_ERROR( MSG, RC )
3078 : RETURN
3079 : ENDIF
3080 : ENDIF
3081 :
3082 : ! Reinitialize error flag. Will be set to 1 or 2 if error occurs,
3083 : ! and to -1 if negative scale factor is ignored.
3084 : ERROR = 0
3085 :
3086 : ! Loop over all latitudes and longitudes
3087 : !$OMP PARALLEL DO &
3088 : !$OMP DEFAULT( SHARED ) &
3089 : !$OMP PRIVATE( I, J, tIdx, TMPVAL, L, LowLL, UppLL, tmpLL, MaskScale )&
3090 : !$OMP COLLAPSE( 2 )&
3091 : !$OMP SCHEDULE( DYNAMIC, 4 )
3092 : DO J = 1, nJ
3093 : DO I = 1, nI
3094 :
3095 : ! ------------------------------------------------------------
3096 : ! If there is a mask associated with this scale factors, check
3097 : ! if this grid box is within or outside of the mask region.
3098 : ! Values that partially fall into the mask region are either
3099 : ! treated as binary (100% inside or outside), or partially
3100 : ! (using the real grid area fractions), depending on the
3101 : ! HEMCO options.
3102 : ! ------------------------------------------------------------
3103 :
3104 : ! Default mask scaling is 1.0 (no mask applied)
3105 : MaskScale = 1.0_sp
3106 :
3107 : ! If there is a mask applied to this scale factor ...
3108 : IF ( ASSOCIATED(MaskDct) ) THEN
3109 : CALL GetMaskVal ( MaskDct, I, J, &
3110 : MaskScale, MaskFractions, RC )
3111 : IF ( RC /= HCO_SUCCESS ) THEN
3112 : ERROR = 4
3113 : EXIT
3114 : ENDIF
3115 : ENDIF
3116 :
3117 : ! We can skip this grid box if mask is completely zero
3118 : IF ( MaskScale <= 0.0_sp ) CYCLE
3119 :
3120 : ! Get current time index for this container and at this location
3121 : tIDx = tIDx_GetIndx( HcoState, ScalDct%Dta, I, J )
3122 : IF ( tIDx < 1 ) THEN
3123 : WRITE(*,*) 'Cannot get time slice index at location ',I,J,&
3124 : ': ', TRIM(ScalDct%cName), tIDx
3125 : ERROR = 3
3126 : EXIT
3127 : ENDIF
3128 :
3129 : ! Check if this is a mask. If so, add mask values to the MASK
3130 : ! array. For now, we assume masks to be binary, i.e. 0 or 1.
3131 : ! We may want to change that in future to also support values
3132 : ! in between. This is especially important when regridding
3133 : ! high resolution masks onto coarser grids!
3134 : ! ------------------------------------------------------------
3135 : IF ( ScalDct%DctType == HCO_DCTTYPE_MASK ) THEN
3136 :
3137 : ! Get mask value
3138 : CALL GetMaskVal ( ScalDct, I, J, &
3139 : TMPVAL, MaskFractions, RC )
3140 : IF ( RC /= HCO_SUCCESS ) THEN
3141 : ERROR = 4
3142 : EXIT
3143 : ENDIF
3144 :
3145 : ! Pass to output mask
3146 : MASK(I,J,:) = MASK(I,J,:) * TMPVAL
3147 :
3148 : ! testing only
3149 : IF ( HCO_IsVerb(HcoState%Config%Err) .AND. I==1 .AND. J==1 ) THEN
3150 : write(MSG,*) 'Mask field ', TRIM(ScalDct%cName), &
3151 : ' found and added to temporary mask.'
3152 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3153 : ENDIF
3154 :
3155 : ! Advance to next grid box
3156 : CYCLE
3157 : ENDIF! DctType=MASK
3158 :
3159 : ! ------------------------------------------------------------
3160 : ! For non-mask fields, apply scale factors to all levels
3161 : ! of the base field individually. If the scale factor
3162 : ! field has more than one vertical level, use the
3163 : ! vertical level closest to the corresponding vertical
3164 : ! level of the base emission field
3165 : ! ------------------------------------------------------------
3166 :
3167 : ! Get lower and upper vertical index
3168 : CALL GetVertIndx( HcoState, BaseDct, isLevDct1, &
3169 : LevDct1, LevDct1_Unit, isLevDct2, &
3170 : LevDct2, LevDct2_Unit, I, &
3171 : J, LowLL, UppLL, RC )
3172 : IF ( RC /= HCO_SUCCESS ) THEN
3173 : ERROR = 1 ! Will cause error
3174 : EXIT
3175 : ENDIF
3176 :
3177 : ! Loop over all vertical levels of the base field
3178 : DO L = LowLL,UppLL
3179 : ! If the vertical level exceeds the number of available
3180 : ! scale factor levels, use the highest available level.
3181 : IF ( L > ScalLL ) THEN
3182 : TmpLL = ScalLL
3183 : ! Otherwise use the same vertical level index.
3184 : ELSE
3185 : TmpLL = L
3186 : ENDIF
3187 :
3188 : ! Get scale factor for this grid box. Use same uniform
3189 : ! value if it's a scalar field
3190 : IF ( ScalDct%Dta%SpaceDim == 1 ) THEN
3191 : TMPVAL = ScalDct%Dta%V2(tidx)%Val(1,1)
3192 : ELSEIF ( ScalDct%Dta%SpaceDim == 2 ) THEN
3193 : TMPVAL = ScalDct%Dta%V2(tidx)%Val(I,J)
3194 : ELSE
3195 : TMPVAL = ScalDct%Dta%V3(tidx)%Val(I,J,TmpLL)
3196 : ENDIF
3197 :
3198 : ! Set missing value to one
3199 : IF ( TMPVAL == HCO_MISSVAL ) TMPVAL = 1.0_sp
3200 :
3201 : ! Eventually apply mask scaling
3202 : IF ( MaskScale /= 1.0_sp ) THEN
3203 : TMPVAL = TMPVAL * MaskScale
3204 : ENDIF
3205 :
3206 : ! For negative scale factor, proceed according to the
3207 : ! negative value setting specified in the configuration
3208 : ! file (NegFlag = 2: use this value):
3209 : IF ( TMPVAL < 0.0_sp .AND. HcoState%Options%NegFlag /= 2 ) THEN
3210 :
3211 : ! NegFlag = 1: ignore and show warning
3212 : IF ( HcoState%Options%NegFlag == 1 ) THEN
3213 : ERROR = -1 ! Will prompt warning
3214 : CYCLE
3215 :
3216 : ! Return w/ error otherwise
3217 : ELSE
3218 : WRITE(*,*) 'Negative scale factor at ',I,J,TmpLL,tidx,&
3219 : ': ', TRIM(ScalDct%cName), TMPVAL
3220 : ERROR = 1 ! Will cause error
3221 : EXIT
3222 : ENDIF
3223 : ENDIF
3224 :
3225 : ! -------------------------------------------------------
3226 : ! Apply scale factor in accordance to field operator
3227 : ! -------------------------------------------------------
3228 :
3229 : ! Oper 1: multiply
3230 : IF ( ScalDct%Oper == 1 ) THEN
3231 : OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL
3232 :
3233 : ! Oper -1: divide
3234 : ELSEIF ( ScalDct%Oper == -1 ) THEN
3235 : ! Ignore zeros to avoid NaN
3236 : IF ( TMPVAL /= 0.0_sp ) THEN
3237 : OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) / TMPVAL
3238 : ENDIF
3239 :
3240 : ! Oper 2: square
3241 : ELSEIF ( ScalDct%Oper == 2 ) THEN
3242 : OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL * TMPVAL
3243 :
3244 : ! Return w/ error otherwise (Oper 3 is only allowed for masks!)
3245 : ELSE
3246 : WRITE(*,*) 'Illegal operator for ', TRIM(ScalDct%cName), ScalDct%Oper
3247 : ERROR = 2 ! Will cause error
3248 : EXIT
3249 : ENDIF
3250 :
3251 : ENDDO !LL
3252 :
3253 : ! Verbose mode
3254 : if ( HCO_IsVerb(HcoState%Config%Err) .and. i == ix .and. j == iy ) then
3255 : write(MSG,*) 'Scale field ', TRIM(ScalDct%cName)
3256 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3257 : write(MSG,*) 'Time slice: ', tIdx
3258 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3259 : write(MSG,*) 'IX, IY: ', IX, IY
3260 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3261 : write(MSG,*) 'Scale factor (IX,IY,L1): ', TMPVAL
3262 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3263 : write(MSG,*) 'Mathematical operation : ', ScalDct%Oper
3264 : CALL HCO_MSG(HcoState%Config%Err,MSG)
3265 : ! write(lun,*) 'Updt (IX,IY,L1): ', OUTARR_3D(IX,IY,1)
3266 : endif
3267 :
3268 : ENDDO !I
3269 : ENDDO !J
3270 : !$OMP END PARALLEL DO
3271 :
3272 : ! error check
3273 : IF ( ERROR > 0 ) THEN
3274 : IF ( ERROR == 1 ) THEN
3275 : MSG = 'Negative scale factor found (aborted): ' // TRIM(ScalDct%cName)
3276 : ELSEIF ( ERROR == 2 ) THEN
3277 : MSG = 'Illegal mathematical operator for scale factor: ' // TRIM(ScalDct%cName)
3278 : ELSEIF ( ERROR == 3 ) THEN
3279 : MSG = 'Encountered negative time index for scale factor: ' // TRIM(ScalDct%cName)
3280 : ELSEIF ( ERROR == 4 ) THEN
3281 : MSG = 'Mask error in ' // TRIM(ScalDct%cName)
3282 : ELSE
3283 : MSG = 'Error when applying scale factor: ' // TRIM(ScalDct%cName)
3284 : ENDIF
3285 : ScalDct => NULL()
3286 : CALL HCO_ERROR( MSG, RC )
3287 : RETURN
3288 : ENDIF
3289 :
3290 : ! eventually prompt warning for negative values
3291 : IF ( ERROR == -1 ) THEN
3292 : MSG = 'Negative scale factor found (ignored): ' // TRIM(ScalDct%cName)
3293 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
3294 : ENDIF
3295 :
3296 : ! Free pointer
3297 : MaskDct => NULL()
3298 :
3299 : ENDDO ! N
3300 : ENDIF ! N > 0
3301 :
3302 : ! Update optional variables
3303 : IF ( PRESENT(UseLL) ) THEN
3304 : UseLL = 1
3305 : IF ( nnLL > 0 ) UseLL = NINT(REAL(TotLL,kind=sp)/REAL(nnLL,kind=sp))
3306 : ENDIF
3307 :
3308 : ! Weight output emissions by mask
3309 : OUTARR_3D = OUTARR_3D * MASK
3310 :
3311 : ! Cleanup and leave w/ success
3312 : ScalDct => NULL()
3313 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
3314 :
3315 : END SUBROUTINE Get_Current_Emissions_Adj
3316 : !EOC
3317 : #endif
3318 : !EOC
3319 : END MODULE HCO_Calc_Mod
|
