-
Notifications
You must be signed in to change notification settings - Fork 37
/
Copy pathhco_calc_mod.F90
3447 lines (3049 loc) · 122 KB
/
hco_calc_mod.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!------------------------------------------------------------------------------
! Harmonized Emissions Component (HEMCO) !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: hco_calc_mod.F90
!
! !DESCRIPTION: Module HCO\_Calc\_Mod contains routines to calculate
! HEMCO core emissions based on the content of the HEMCO EmisList
! object. All emissions are in [kg/m2/s].
!\\
!\\
! Emissions for the current datetime are calculated by multiplying base
! emissions fields with the associated scale factors. Different
! inventories are merged/overlayed based upon the category and hierarchy
! attributes assigned to the individual base fields. Within the same
! category, fields of higher hierarchy overwrite lower-hierarchy fields.
! Emissions of different categories are always added.
!\\
!\\
! The assembled emission array is written into the corresponding emission
! rates array of the HEMCO state object: HcoState%Spc(HcoID)%Emis, where
! HcoID denotes the corresponding species ID. Emis covers dimension lon,
! lat, lev on the HEMCO grid, i.e. unlike the emission arrays in EmisList
! that only cover the levels defined in the source files, Emis extends
! over all vertical model levels.
!\\
!\\
! Negative emissions are not supported and are ignored. Surface
! deposition velocities are stored in HcoState%Spc(HcoID)%Depv and can
! be added therein.
!\\
!\\
! In addition to emissions and surface deposition rates, HEMCO also
! supports concentrations (kg/m3). Data is automatically written into
! the concentration array HcoState%Spc(HcoID)%Conc if the source data
! is marked as being concentration data (i.e. if Dta%IsConc is .TRUE.).
! The latter is automatically determined by HEMCO based upon the data
! units.
!\\
!\\
! All emission calculation settings are passed through the HcoState
! options object (HcoState%Options). These include:
!
! \begin{itemize}
! \item ExtNr: extension number to be considered.
! \item SpcMin: lower species ID (HEMCO ID) to be considered.
! \item SpcMax: upper species ID (HEMCO ID) to be considered. If set
! to -1, all species above or equal to SpcMin are considered.
! \item CatMin: lower emission category to be considered.
! \item CatMax: upper emission category to be considered. If set to
! -1, all categories above or equal to CatMin are considered.
! \item FillBuffer: if set to TRUE, the emissions will be written into
! buffer array HcoState%Buffer3D instead of HcoState%Spc(ID)%Emis.
! If this option is enabled, only one species can be calculated at
! a time (by setting SpcMin/SpcMax, CatMin/CatMax and/or ExtNr
! accordingly). This option is useful for extensions, e.g. if
! additional scalings need to be done on some emission fields
! assembled by HEMCO (e.g. PARANOX extension).
! \end{itemize}
!
! !INTERFACE:
!
MODULE HCO_Calc_Mod
!
! !USES:
!
USE HCO_Diagn_Mod
USE HCO_Error_Mod
USE HCO_Types_Mod
USE HCO_DataCont_Mod, ONLY : Pnt2DataCont
IMPLICIT NONE
PRIVATE
!
! !PUBLIC MEMBER FUNCTIONS:
!
PUBLIC :: HCO_CalcEmis
PUBLIC :: HCO_CheckDepv
PUBLIC :: HCO_EvalFld
PUBLIC :: HCO_MaskFld
#ifdef ADJOINT
PUBLIC :: GET_CURRENT_EMISSIONS_ADJ
#endif
!
! !PRIVATE MEMBER FUNCTIONS:
!
PRIVATE :: GET_CURRENT_EMISSIONS
PRIVATE :: GetMaskVal
PRIVATE :: GetDilFact
PRIVATE :: GetVertIndx
PRIVATE :: GetIdx
!
! !PARAMETER
!
! Mask threshold. All mask values below this value will be evaluated
! as zero (= outside of mask), and all values including and above this
! value as inside the mask. This is only of relevance if the MaskFractions
! option is false. If MaskFractions is true, the fractional mask values are
! considered, e.g. a grid box can contribute 40% to a mask region, etc.
! The MaskFractions toggle can be set in the settings section of the HEMCO
! configuration file (Use mask fractions: true/false). It defaults to false.
REAL(sp), PARAMETER :: MASK_THRESHOLD = 0.5_sp
!
! ============================================================================
!
! !REVISION HISTORY:
! 25 Aug 2012 - C. Keller - Initial version.
! See https://github.com/geoschem/hemco for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
INTERFACE HCO_EvalFld
MODULE PROCEDURE HCO_EvalFld_2D
MODULE PROCEDURE HCO_EvalFld_3D
END INTERFACE HCO_EvalFld
CONTAINS
!EOC
!------------------------------------------------------------------------------
! Harmonized Emissions Component (HEMCO) !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: HCO_CalcEmis
!
! !DESCRIPTION: Subroutine HCO\_CalcEmis calculates the 3D emission
! fields at current datetime for the specified species, categories, and
! extension number.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC )
!
! !USES:
!
USE HCO_STATE_MOD, ONLY : HCO_State
USE HCO_ARR_MOD, ONLY : HCO_ArrAssert
USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont
USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined
USE HCO_Scale_Mod, ONLY : HCO_ScaleArr
!
! !INPUT PARAMETERS:
!
LOGICAL, INTENT(IN ) :: UseConc ! Use concentration fields?
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
INTEGER, INTENT(INOUT) :: RC ! Return code
!
! !REVISION HISTORY:
! 25 Aug 2012 - C. Keller - Initial Version
! See https://github.com/geoschem/hemco for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
! Working pointers: list and data container
TYPE(ListCont), POINTER :: Lct
TYPE(DataCont), POINTER :: Dct
! Temporary emission arrays
REAL(hp), POINTER :: OutArr(:,:,:) => NULL()
REAL(hp), TARGET :: SpcFlx( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
REAL(hp), TARGET :: CatFlx( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
REAL(hp), TARGET :: TmpFlx( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
REAL(hp) :: Mask ( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
REAL(hp) :: HirFlx( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
REAL(hp) :: HirMsk( HcoState%NX, &
HcoState%NY, &
HcoState%NZ )
! Integers
INTEGER :: ThisSpc, PrevSpc ! current and previous species ID
INTEGER :: ThisCat, PrevCat ! current and previous category
INTEGER :: ThisHir, PrevHir ! current and previous hierarchy
INTEGER :: SpcMin, SpcMax ! range of species to be considered
INTEGER :: CatMin, CatMax ! range of categories to be considered
INTEGER :: ExtNr ! Extension Nr to be used
INTEGER :: nI, nJ, nL
INTEGER :: nnSpec, FLAG
LOGICAL :: Found, DoDiagn, EOL, UpdateCat
! For error handling & verbose mode
CHARACTER(LEN=255) :: MSG, LOC
! testing / debugging
integer :: ix,iy
!=================================================================
! HCO_CalcEmis begins here!
!=================================================================
! testing only
ix = 30
iy = 34
! Initialize
LOC = 'HCO_CalcEmis (HCO_CALC_MOD.F90)'
Lct => NULL()
Dct => NULL()
! Enter routine
CALL HCO_ENTER (HcoState%Config%Err, LOC, RC )
IF(RC /= HCO_SUCCESS) RETURN
!-----------------------------------------------------------------
! Initialize variables
!-----------------------------------------------------------------
! Initialize
SpcFlx(:,:,:) = 0.0_hp
CatFlx(:,:,:) = 0.0_hp
HirFlx(:,:,:) = 0.0_hp
HirMsk(:,:,:) = 0.0_hp
PrevSpc = -1
PrevHir = -1
PrevCat = -1
nnSpec = 0
! Pass emission grid dimensions
nI = HcoState%NX
nJ = HcoState%NY
nL = HcoState%NZ
! Pass calculation options
SpcMin = HcoState%Options%SpcMin !Lower species ID
SpcMax = HcoState%Options%SpcMax !Upper species ID
CatMin = HcoState%Options%CatMin !Lower emission category
CatMax = HcoState%Options%CatMax !Upper emission category
ExtNr = HcoState%Options%ExtNr !Extension number
DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics?
! Verbose mode
IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN
WRITE (MSG, *) 'Run HEMCO calculation w/ following options:'
CALL HCO_MSG ( HcoState%Config%Err, MSG )
WRITE (MSG, "(A20,I5)") 'Extension number:', ExtNr
CALL HCO_MSG ( HcoState%Config%Err, MSG )
WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax
CALL HCO_MSG ( HcoState%Config%Err, MSG )
WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax
CALL HCO_MSG ( HcoState%Config%Err, MSG )
WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn
CALL HCO_MSG ( HcoState%Config%Err, MSG )
ENDIF
!=================================================================
! Walk through all containers of EmisList and determine the
! emissions for all containers that qualify for calculation.
! The containers in EmisList are sorted by species, category and
! hierarchy. This enables a straightforward, piece-by-piece
! assembly of the final emission array (start with lowest
! hierarchy emissions, then overwrite piece-by-piece with higher
! hierarchy values).
!=================================================================
! Point to the head of the emissions linked list
EOL = .FALSE. ! End of list
Lct => NULL()
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
! Do until end of EmisList (==> loop over all emission containers)
DO
! Have we reached the end of the list?
IF ( FLAG /= HCO_SUCCESS ) THEN
EOL = .TRUE.
ELSE
EOL = .FALSE.
ENDIF
! ------------------------------------------------------------
! Select container and update all working variables & arrays.
! ------------------------------------------------------------
IF ( .NOT. EOL ) THEN
! Dct is the current data container
Dct => Lct%Dct
! Check if this is a base field
IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Sanity check: Make sure this container holds data.
! 'Empty' containers are possible if the simulation time
! is outside of the specified data time range and time
! slice cycling is deactivated (CycleFlag > 1).
IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Check if this is the specified extension number
IF ( Dct%ExtNr /= ExtNr ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Advance to next container if the species ID is outside
! the specified species range (SpcMin - SpcMax). Consider
! all species above SpcMin if SpcMax is negative!
IF( ( Dct%HcoID < SpcMin ) .OR. &
( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Advance to next emission field if the emission category of
! the current container is outside of the specified species
! range (CatMin - CatMax). Consider all categories above CatMin
! if CatMax is negative!
IF( ( Dct%Cat < CatMin ) .OR. &
( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Check if this container holds data in the desired unit format,
! i.e. concentration data if UseConc is enabled, emission data
! otherwise.
IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN
CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
CYCLE
ENDIF
! Update working variables
ThisSpc = Dct%HcoID
ThisCat = Dct%Cat
ThisHir = Dct%Hier
! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir
! to make sure that emissions are added to HEMCO in the section
! below!
ELSE
ThisSpc = -1
ThisCat = -1
ThisHir = -1
ENDIF
!--------------------------------------------------------------------
! Before computing emissions of current data container make sure that
! emissions of previous container are properly archived.
!--------------------------------------------------------------------
! Add emissions on hierarchy level to the category flux array. Do
! this only if this is a new species, a new category or a new
! hierarchy level.
! Note: no need to add to diagnostics because hierarchy level
! diagnostics are filled right after computing the emissions of
! a given data container (towards the end of the DO loop).
IF ( (ThisHir /= PrevHir) .OR. &
(ThisSpc /= PrevSpc) .OR. &
(ThisCat /= PrevCat) ) THEN
! Add hierarchy level emissions to category array over the
! covered regions.
CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx
! Reset
HirFlx = 0.0_hp
HirMsk = 0.0_hp
ENDIF
!--------------------------------------------------------------------
! If this is a new species or category, pass the previously collected
! emissions to the species array. Update diagnostics at category level.
! Skip this step for first species, i.e. if PrevSpc is still -1.
!--------------------------------------------------------------------
UpdateCat = .FALSE.
IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE.
IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE.
IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE.
IF ( UpdateCat ) THEN
! CatFlx holds the emissions for this category. Pass this to
! the species array SpcFlx.
SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:)
! verbose
IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN
WRITE(MSG,*) 'Added category emissions to species array: '
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'Species : ', PrevSpc
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'Category : ', PrevCat
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx)
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx)
CALL HCO_MSG(HcoState%Config%Err,MSG)
ENDIF
! Add category emissions to diagnostics at category level
! (only if defined in the diagnostics list).
IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN
! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
! are updated, including those manually defined in other models
! (mps, 11/30/21)
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=PrevCat, Hier=-1, HcoID=PrevSpc, &
AutoFill=1, Array3D=CatFlx, COL=-1, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
RETURN
ENDIF
#ifdef ADJOINT
IF (HcoState%IsAdjoint) THEN
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=PrevCat, Hier=-1, HcoID=PrevSpc, &
AutoFill=1, Array3D=CatFlx, &
COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
RETURN
ENDIF
ENDIF
#endif
ENDIF
! Reset CatFlx array and the previously used hierarchy
! ==> Emission hierarchies are only important within the
! same category, hence always start over at lowest hierarchy
! when entering a new category.
CatFlx(:,:,:) = 0.0_hp
PrevHir = -1
ENDIF
!--------------------------------------------------------------------
! If this is a new species, pass previously calculated emissions
! to the final emissions array in HcoState.
! Update diagnostics at extension number level.
! Don't do before first emission calculation, i.e. if PrevSpc
! is still the initialized value of -1!
!--------------------------------------------------------------------
IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN
! Add to OutArr
OutArr(:,:,:) = OutArr(:,:,:) + SpcFlx(:,:,:)
! testing only
IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN
WRITE(MSG,*) 'Added total emissions to output array: '
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'Species: ', PrevSpc
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx)
CALL HCO_MSG(HcoState%Config%Err,MSG)
WRITE(MSG,*) 'OutArr : ', SUM(OutArr)
CALL HCO_MSG(HcoState%Config%Err,MSG)
ENDIF
! Add to diagnostics at extension number level.
! The same diagnostics may be updated multiple times during
! the same time step, continuously adding emissions to it.
IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN
! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
! are updated, including those manually defined in other models
! (mps, 11/30/21)
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=-1, Hier=-1, HcoID=PrevSpc, &
AutoFill=1,Array3D=SpcFlx, COL=-1, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
RETURN
ENDIF
#ifdef ADJOINT
IF (HcoState%IsAdjoint) THEN
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=-1, Hier=-1, HcoID=PrevSpc, &
AutoFill=1,Array3D=SpcFlx, &
COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
RETURN
ENDIF
ENDIF
#endif
ENDIF
! Reset arrays and previous hierarchy.
SpcFlx(:,:,:) = 0.0_hp
PrevCat = -1
PrevHir = -1
OutArr => NULL()
ENDIF
!--------------------------------------------------------------------
! Exit DO loop here if end of list
!--------------------------------------------------------------------
IF ( EOL ) EXIT
!--------------------------------------------------------------------
! Update/archive information on species level if needed
!--------------------------------------------------------------------
IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN
! Update number of species for which emissions have been
! calculated.
nnSpec = nnSpec + 1
! To write emissions into temporary array, make OutArr point
! to the buffer array HcoState%Buffer3D.
IF ( HcoState%Options%FillBuffer ) THEN
! Cannot use temporary array for more than one species!
IF ( nnSpec > 1 ) THEN
MSG = 'Cannot fill buffer for more than one species!'
CALL HCO_ERROR( MSG, RC )
RETURN
ENDIF
! Point to array and check allocation status as well as
! array size.
OutArr => HcoState%Buffer3D%Val
IF ( .NOT. ASSOCIATED( OutArr ) ) THEN
MSG = 'Buffer array is not associated'
CALL HCO_ERROR( MSG, RC )
RETURN
ENDIF
IF ( (SIZE(OutArr,1) /= nI) .OR. &
(SIZE(OutArr,2) /= nJ) .OR. &
(SIZE(OutArr,3) /= nL) ) THEN
MSG = 'Buffer array has wrong dimension!'
CALL HCO_ERROR( MSG, RC )
RETURN
ENDIF
! To write emissions directly into HcoState, make OutArr
! point to current species' array in HcoState. Use emission
! array for emissions, and concentration array for concentrations.
ELSE
! For concentrations:
IF ( UseConc ) THEN
CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Conc, &
nI, nJ, nL, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
RETURN
ENDIF
OutArr => HcoState%Spc(ThisSpc)%Conc%Val
! For emissions:
ELSE
CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Emis, &
nI, nJ, nL, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
RETURN
ENDIF
OutArr => HcoState%Spc(ThisSpc)%Emis%Val
ENDIF
ENDIF
! verbose mode
IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN
WRITE(MSG,*) 'Calculating emissions for species ', &
TRIM(HcoState%Spc(ThisSpc)%SpcName)
CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' )
ENDIF
ENDIF
!--------------------------------------------------------------------
! Get current emissions and write into TmpFlx array. The array Mask
! denotes all valid grid boxes for this inventory.
!--------------------------------------------------------------------
TmpFlx(:,:,:) = 0.0_hp
CALL GET_CURRENT_EMISSIONS( HcoState, Dct, nI, nJ, nL, TmpFlx, Mask, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
RETURN
ENDIF
! Eventually add universal scale factor
CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
RETURN
ENDIF
! Check for negative values according to the corresponding setting
! in the configuration file: 2 means allow negative values, 1 means
! set to zero and prompt a warning, else return with error.
IF ( HcoState%Options%NegFlag /= 2 ) THEN
IF ( ANY(TmpFlx < 0.0_hp) ) THEN
! Set to zero and prompt warning
IF ( HcoState%Options%NegFlag == 1 ) THEN
WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp
MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName)
CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
! Return with error
ELSE
MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // &
'To allow negatives, edit settings in the configuration file.'
CALL HCO_ERROR( MSG, RC )
RETURN
ENDIF
ENDIF
ENDIF
! ------------------------------------------------------------
! Collect all emissions of the same category (and species) on
! the hierarchy level into array HirFlx. HirMsk contains the
! combined covered region. That is, if there are two regional
! inventories with the same hierarchy HirMsk will cover both
! of these regions.
! The specified field hierarchies determine whether the
! temporary emissions are added (if hierarchy is the same
! as the previously used hierarchy), or if they overwrite the
! previous values in HirFlx (if hierarchy is higher than the
! previous hierarchy).
! ------------------------------------------------------------
! Add emissions to the hierarchy array HirFlx if this hierarchy
! is the same as previous hierarchy
IF ( ThisHir == PrevHir ) THEN
HirFlx = HirFlx + TmpFlx
HirMsk = HirMsk + Mask
! Make sure mask values do not exceed 1.0
WHERE(HirMsk > 1.0 ) HirMsk = 1.0
! If hierarchy is larger than those of the previously used
! fields, overwrite HirFlx with new values.
ELSE
HirFlx = TmpFlx
HirMsk = Mask
ENDIF
! Update diagnostics at hierarchy level. Make sure that only
! positive values are used.
! The same diagnostics may be updated multiple times
! during the same time step, continuously adding
! emissions to it.
! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's
! ok to keep negative values (ckeller, 7/12/15).
IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN
! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
! are updated, including those manually defined in other models
! (mps, 11/30/21)
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, &
AutoFill=1, Array3D=TmpFlx, &
COL=-1, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
RETURN
ENDIF
#ifdef ADJOINT
IF (HcoState%IsAdjoint) THEN
! I don't know why I chose collection=-1 instead of
! collection=HcoState%Diagn%HcoDiagnIDAdjoint like in the other
! parts of the adjoint code here, but it's what worked in the
! old repo so I'm keeping it for now. May need to change
CALL Diagn_Update( HcoState, ExtNr=ExtNr, &
Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, &
AutoFill=1, Array3D=TmpFlx, &
COL=-1, RC=RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
RETURN
ENDIF
ENDIF
#endif
ENDIF
! Update previously used species, category and hierarchy
PrevSpc = ThisSpc
PrevCat = ThisCat
PrevHir = ThisHir
! Advance to next emission container
CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG )
ENDDO ! Loop over EmisList
! Make sure internal pointers are nullified
Lct => NULL()
Dct => NULL()
OutArr => NULL()
! verbose
IF ( HCO_IsVerb(HcoState%Config%Err,1) ) THEN
WRITE (MSG, *) 'HEMCO emissions successfully calculated!'
CALL HCO_MSG ( HcoState%Config%Err, MSG )
ENDIF
! Leave w/ success
CALL HCO_LEAVE ( HcoState%Config%Err, RC )
END SUBROUTINE HCO_CalcEmis
!EOC
!------------------------------------------------------------------------------
! Harmonized Emissions Component (HEMCO) !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: HCO_CheckDepv
!
! !DESCRIPTION: Subroutine HCO\_CheckDepv is a simple routine to check the
! dry deposition frequency value. This is to avoid unrealistically high
! deposition frequencies that may occur if grid box concentrations are very
! low. The deposition frequency is limited to a value that will make sure
! that the drydep exponent ( exp( -depfreq * dt ) ) is still small enough to
! remove all species mass. The maximum limit of depfreq * dt can be defined
! as a HEMCO option (MaxDepExp). Its default value is 20.0.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE HCO_CheckDepv( HcoState, Depv, RC )
!
! !USES:
!
USE HCO_STATE_MOD, ONLY : HCO_State
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
REAL(hp), INTENT(INOUT) :: Depv ! Deposition velocity
INTEGER, INTENT(INOUT) :: RC ! Return code
!
! !REVISION HISTORY:
! 11 Mar 2015 - C. Keller - Initial Version
! See https://github.com/geoschem/hemco for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
REAL(hp) :: ExpVal
!=================================================================
! HCO_CheckDepv begins here!
!=================================================================
ExpVal = Depv * HcoState%TS_EMIS
IF ( ExpVal > HcoState%Options%MaxDepExp ) THEN
Depv = HcoState%Options%MaxDepExp / HcoState%TS_EMIS
ENDIF
END SUBROUTINE HCO_CheckDepv
!EOC
!------------------------------------------------------------------------------
! Harmonized Emissions Component (HEMCO) !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: Get_Current_Emissions
!
! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current
! emissions for the specified emission container.
! This subroutine is only called by HCO\_CalcEmis and for base emission
! containers, i.e. containers of type 1.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE Get_Current_Emissions( HcoState, BaseDct, nI, nJ, &
nL, OUTARR_3D, MASK, RC, UseLL )
!
! !USES:
!
USE HCO_State_Mod, ONLY : HCO_State
USE HCO_tIdx_MOD, ONLY : tIDx_GetIndx
USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: nI ! # of lons
INTEGER, INTENT(IN) :: nJ ! # of lats
INTEGER, INTENT(IN) :: nL ! # of levs
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
TYPE(DataCont), POINTER :: BaseDct ! base emission
! container
REAL(hp), INTENT(INOUT) :: OUTARR_3D(nI,nJ,nL) ! output array
REAL(hp), INTENT(INOUT) :: MASK (nI,nJ,nL) ! mask array
INTEGER, INTENT(INOUT) :: RC
!
! !OUTPUT PARAMETERS:
!
INTEGER, INTENT( OUT), OPTIONAL :: UseLL
!
! !REMARKS:
! This routine uses multiple loops over all grid boxes (base emissions
! and scale factors use separate loops). In an OMP environment, this approach
! seems to be faster than using only one single loop (but repeated calls to
! point to containers, etc.). The alternative approach is used in routine
! Get\_Current\_Emissions\_B at the end of this module and may be employed
! on request.
!
! !REVISION HISTORY:
! 25 Aug 2012 - C. Keller - Initial Version
! See https://github.com/geoschem/hemco for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
! Pointers
TYPE(DataCont), POINTER :: ScalDct
TYPE(DataCont), POINTER :: MaskDct
TYPE(DataCont), POINTER :: LevDct1
TYPE(DataCont), POINTER :: LevDct2
! Scalars
REAL(sp) :: TMPVAL, MaskScale
REAL(hp) :: DilFact
REAL(hp) :: ScalFact
INTEGER :: tIDx, IDX
INTEGER :: totLL, nnLL
INTEGER :: I, J, L, N
INTEGER :: LowLL, UppLL, ScalLL, TmpLL
INTEGER :: ERROR
CHARACTER(LEN=255) :: MSG, LOC
LOGICAL :: NegScalExist
LOGICAL :: MaskFractions
LOGICAL :: isLevDct1
LOGICAL :: isLevDct2
LOGICAL :: isMaskDct
LOGICAL :: isPblHt
LOGICAL :: isBoxHt
INTEGER :: LevDct1_Unit
INTEGER :: LevDct2_Unit
! testing only
INTEGER, PARAMETER :: IX=25, IY=25
!=================================================================
! GET_CURRENT_EMISSIONS begins here
!=================================================================
! Initialize
ScalDct => NULL()
MaskDct => NULL()
LOC = 'GET_CURRENT_EMISSIONS (hco_calc_mod.F90)'
! Enter
CALL HCO_ENTER(HcoState%Config%Err, LOC, RC )
IF(RC /= HCO_SUCCESS) RETURN
! Check if container contains data
IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN
MSG = 'Array not defined: ' // TRIM(BaseDct%cName)
CALL HCO_ERROR( MSG, RC )
RETURN
ENDIF
! Initialize mask. By default, assume that we use all grid boxes.
MASK(:,:,:) = 1.0_hp
MaskFractions = HcoState%Options%MaskFractions
! Verbose
IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN
WRITE(MSG,*) 'Evaluate field ', TRIM(BaseDct%cName)
CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
ENDIF
! Put check for PBLHEIGHT here (bmy, 3/4/21)
#if !defined ( ESMF_ )
IF ( .NOT. ASSOCIATED(HcoState%Grid%PBLHEIGHT%Val) ) THEN
MSG = 'PBLHEIGHT (in meters) is missing in HEMCO state'
CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
RETURN
ENDIF
#endif
! ----------------------------------------------------------------
! Set base emissions
! ----------------------------------------------------------------
! Initialize ERROR. Will be set to 1 if error occurs below
ERROR = 0
! Initialize variables to compute average vertical level index
totLL = 0.0
nnLL = 0.0
!-----------------------------------------------------------------
! Check for level index containers
! Move error checks here, outside of the parallel DO loop
!-----------------------------------------------------------------
IF ( BaseDct%levScalID1 > 0 ) THEN
CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
RETURN
ENDIF
ELSE
LevDct1 => NULL()
ENDIF
IF ( BaseDct%levScalID2 > 0 ) THEN
CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC )
IF ( RC /= HCO_SUCCESS ) THEN
CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
RETURN
ENDIF
ELSE
LevDct2 => NULL()
ENDIF
! Test whether LevDct1 and LevDct2 are associated
isLevDct1 = ASSOCIATED( LevDct1 )
isLevDct2 = ASSOCIATED( LevDct2 )
! Get the units of LevDct1 (if it exists)
IF ( isLevDct1 ) THEN
LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
IF ( LevDct1_Unit < 0 ) THEN
MSG = 'LevDct1 units are not defined!'
CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
RC = HCO_FAIL
RETURN
ENDIF
ELSE
LevDct1_Unit = -1
ENDIF
! Get the units of LevDct2 (if it exists)
IF ( isLevDct2 ) THEN
LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
IF ( LevDct2_Unit < 0 ) THEN
MSG = 'LevDct2_Units are not defined!'
CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
RETURN
ENDIF
ELSE
LevDct2_Unit = -1
ENDIF
! Throw an error if boxheight is missing and the units are in meters
IF ( LevDct1_Unit == HCO_EMISL_M .or. &
LevDct2_Unit == HCO_EMISL_M ) THEN
IF ( .NOT. ASSOCIATED(HcoState%Grid%BXHEIGHT_M%Val) ) THEN
MSG = 'Boxheight (in meters) is missing in HEMCO state'
CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
RETURN
ENDIF
ENDIF
! Throw an error if boxheight is missing and the units are in PBL frac
IF ( LevDct1_Unit == HCO_EMISL_PBL .or. &
LevDct2_Unit == HCO_EMISL_PBL ) THEN
IF ( .NOT. ASSOCIATED(HcoState%Grid%PBLHEIGHT%Val) ) THEN
MSG = 'Boundary layer height is missing in HEMCO state'
CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
RETURN
ENDIF
ENDIF
!------------------------------------------------------------------------
! Loop over all latitudes and longitudes
!
! NOTE: It is OK to exit from the "I" loop, because only
! the "J" loop is being parallelized (bmy, 3/8/21)
!------------------------------------------------------------------------
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, tIdx, TMPVAL, DilFact, LowLL, UppLL )&
!$OMP REDUCTION( +:totLL )&
!$OMP REDUCTION( +:nnLL )
DO J = 1, nJ
DO I = 1, nI
! Zero private variables for safety's sake
tmpVal = 0.0_hp
dilFact = 0.0_hp
lowLL = 0
uppLL = 0