forked from geoschem/geos-chem
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmixing_mod.F90
1021 lines (883 loc) · 40.7 KB
/
mixing_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
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: mixing_mod.F90
!
! !DESCRIPTION: Module mixing\_mod.F90 is a wrapper module for the PBL mixing
! in GEOS-Chem.
!\\
!\\
! !INTERFACE:
!
MODULE Mixing_Mod
!
! !USES:
!
USE Precision_Mod
IMPLICIT NONE
PRIVATE
!
! !PUBLIC MEMBER FUNCTIONS:
!
PUBLIC :: DO_MIXING
PUBLIC :: DO_TEND
!
! !REVISION HISTORY:
! 04 Mar 2015 - C. Keller - Initial version.
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
CONTAINS
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: do_mixing
!
! !DESCRIPTION: Subroutine DO\_MIXING performs the PBL mixing.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE Do_Mixing( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, RC )
!
! !USES:
!
USE ErrCode_Mod
USE Input_Opt_Mod, ONLY : OptInput
USE Pbl_Mix_Mod, ONLY : Do_Full_Pbl_Mixing
USE State_Chm_Mod, ONLY : ChmState
USE State_Diag_Mod, ONLY : DgnState
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
USE Vdiff_Mod, ONLY : Do_Vdiff
!
! !INPUT PARAMETERS:
!
TYPE(OptInput), INTENT(IN ) :: Input_Opt ! Input Options
TYPE(GrdState), INTENT(IN ) :: State_Grid ! Grid State
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(MetState), INTENT(INOUT) :: State_Met ! Meteorology State
TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State
TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diagnostics State
INTEGER, INTENT(INOUT) :: RC ! Failure or success
!
! !REMARKS
! (A) While all dry deposition rates are calculated either in
! DO_PBL_MIX2 or DO_TEND, settling of aerosols is still
! computed in the dust/seasalt modules.
!
! !REVISION HISTORY:
! 04 Mar 2015 - C. Keller - Initial version
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! LOCAL VARIABLES:
!
! Scalars
LOGICAL :: OnlyAbovePBL
CHARACTER(LEN=255) :: ErrMsg, ThisLoc
!=======================================================================
! DO_MIXING begins here!
!=======================================================================
! Initialize
RC = GC_SUCCESS
ErrMsg = ''
ThisLoc = ' -> at DO_MIXING (in module GeosCore/mixing_mod.F90)'
!-----------------------------------------------------------------------
! Do non-local PBL mixing. This will apply the species tendencies
! (emission fluxes and dry deposition rates) below the PBL.
! This is done for all species with defined emissions / dry
! deposition rates, including dust.
!
! Set OnlyAbovePBL flag (used below by DO_TEND) to indicate that
! fluxes within the PBL have already been applied.
! ----------------------------------------------------------------------
IF ( Input_Opt%LTURB .AND. Input_Opt%LNLPBL ) THEN
!--------------------------------------------------------------------
! %%%%% HISTORY (aka netCDF diagnostics) %%%%%
!
! Initialize the diagnostic array for the History Component. This will
! prevent leftover values from being carried over to this timestep.
! (For example, if on the last iteration, the PBL height was higher than
! it is now, then we will have stored drydep fluxes up to that height,
! so we need to zero these out.)
!--------------------------------------------------------------------
! Non-local mixing
CALL Do_Vdiff( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, RC )
! Trap potential errors
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Error encountred in "DO_PBL_MIX_2"!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
! Fluxes in PBL have been applied (non-local PBL mixing)
OnlyAbovePBL = .TRUE.
ELSE
! Fluxes in PBL have not been applied (full PBL mixing)
OnlyAbovePBL = .FALSE.
ENDIF
!-----------------------------------------------------------------------
! Apply tendencies. This will apply dry deposition rates and
! emission fluxes below the PBL if it has not yet been done
! via the non-local PBL mixing. It also adds the emissions above
! the PBL to the species array. Emissions of some species may be
! capped at the tropopause to avoid build-up in stratosphere.
!-----------------------------------------------------------------------
! Apply tendencies
CALL DO_TEND( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, OnlyAbovePBL, RC )
! Trap potential error
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Error encountred in "DO_TEND"!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
!-----------------------------------------------------------------------
! Do full pbl mixing. This fully mixes the updated species
! concentrations within the PBL.
!
! Now also archive concentrations and calculate turbulence
! tendencies (ckeller, 7/15/2015)
!-----------------------------------------------------------------------
IF ( Input_Opt%LTURB .AND. .NOT. Input_Opt%LNLPBL ) THEN
! Full PBL mixing
CALL Do_Full_Pbl_Mixing( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, RC )
! Trap potential error
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Error encountred in "DO_PBL_MIX"!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
ENDIF
END SUBROUTINE DO_MIXING
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: do_tend
!
! !DESCRIPTION: Subroutine DO\_TEND adds the species tendencies (dry deposition
! and emissions) to the species array.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE DO_TEND( Input_Opt, State_Chm, State_Diag, State_Grid, &
State_Met, OnlyAbovePBL, RC, DT )
!
! !USES:
!
USE Diagnostics_Mod, ONLY : Compute_Budget_Diagnostics
USE ErrCode_Mod
USE ERROR_MOD, ONLY : SAFE_DIV
USE GET_NDEP_MOD, ONLY : SOIL_DRYDEP
USE HCO_Utilities_GC_Mod, ONLY : HCO_GC_GetDiagn
USE HCO_Utilities_GC_Mod, ONLY : GetHcoValEmis, GetHcoValDep, InquireHco
USE HCO_Utilities_GC_Mod, ONLY : LoadHcoValEmis, LoadHcoValDep
USE Input_Opt_Mod, ONLY : OptInput
USE PhysConstants, ONLY : AVO
USE Species_Mod, ONLY : Species
USE State_Chm_Mod, ONLY : ChmState
USE State_Chm_Mod, ONLY : Ind_
USE State_Diag_Mod, ONLY : DgnState
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
USE TIME_MOD, ONLY : GET_TS_DYN, GET_TS_CONV, GET_TS_CHEM
USE UnitConv_Mod
#ifdef MODEL_CLASSIC
use hco_utilities_gc_mod, only: TMP_MDL ! danger
#endif
!
! !INPUT PARAMETERS:
!
TYPE(OptInput), INTENT(IN ) :: Input_Opt ! Input opts
TYPE(MetState), INTENT(IN ) :: State_Met ! Met state
TYPE(GrdState), INTENT(IN ) :: State_Grid ! Grid state
LOGICAL, INTENT(IN ) :: OnlyAbovePBL ! Only above PBL?
REAL(fp), INTENT(IN ), OPTIONAL :: DT ! Time step [s]
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry state
TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diags State
INTEGER, INTENT(INOUT) :: RC ! Success/Failure
!
! !REVISION HISTORY:
! 04 Mar 2015 - C. Keller - Initial version
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
! Scalars
INTEGER :: I, J, L, L1, L2, N, D, NN, NA, nAdvect, S
INTEGER :: DRYDEPID, origUnit
INTEGER :: PBL_TOP, DRYD_TOP, EMIS_TOP
REAL(fp) :: TS, TMP, FRQ, RKT, FRAC, FLUX, AREA_M2
REAL(fp) :: MWkg, DENOM
LOGICAL :: FND
LOGICAL :: PBL_DRYDEP, LINEAR_CHEM, ChemGridOnly
LOGICAL :: LEMIS, LDRYD
LOGICAL :: DryDepSpec, EmisSpec
REAL(f8) :: DT_Tend
! PARANOX loss fluxes (kg/m2/s). These are obtained from the
! HEMCO PARANOX extension via the diagnostics module.
REAL(fp) :: PNOXLOSS
REAL(f4), POINTER :: Ptr2D (:,:) => NULL()
REAL(f4), POINTER :: PNOXLOSS_O3 (:,:)
REAL(f4), POINTER :: PNOXLOSS_HNO3(:,:)
! SAVEd scalars (defined on first call only)
LOGICAL, SAVE :: FIRST = .TRUE.
INTEGER, SAVE :: id_MACR, id_RCHO, id_ACET, id_ALD2
INTEGER, SAVE :: id_ALK4, id_C2H6, id_C3H8, id_CH2O
INTEGER, SAVE :: id_PRPE, id_O3, id_HNO3, id_BrO
INTEGER, SAVE :: id_Br2, id_Br, id_HOBr, id_HBr
INTEGER, SAVE :: id_BrNO3, id_CH4_SAB, id_CO2
! Pointers and objects
TYPE(Species), POINTER :: SpcInfo
REAL(fp), POINTER :: DepFreq(:,:,: ) ! IM, JM, nDryDep
! Temporary save for total ch4 (Xueying Yu, 12/08/2017)
REAL(fp) :: total_ch4_pre_soil_absorp(State_Grid%NX, &
State_Grid%NY, &
State_Grid%NZ)
! Strings
CHARACTER(LEN=255) :: ErrMsg, ErrorMsg, ThisLoc
#ifdef ADJOINT
LOGICAL :: IS_ADJ
#endif
!=================================================================
! DO_TEND begins here!
!=================================================================
! Assume success
RC = GC_SUCCESS
ErrMsg = ''
ThisLoc = ' -> at DO_TEND (in module GeosCore/mixing_mod.F90)'
! Special case that there is no dry deposition and emissions
IF ( .NOT. Input_Opt%LDRYD .AND. .NOT. Input_Opt%DoEmissions ) RETURN
! Initialize
LINEAR_CHEM = Input_Opt%LINEAR_CHEM
LEMIS = Input_Opt%DoEmissions
LDRYD = Input_Opt%LDRYD
PBL_DRYDEP = Input_Opt%PBL_DRYDEP
nAdvect = State_Chm%nAdvect
! Initialize pointer
SpcInfo => NULL()
DepFreq => State_Chm%DryDepFreq
PNOxLoss_O3 => NULL()
PNOxLoss_HNO3 => NULL()
!------------------------------------------------------------------------
! Emissions/dry deposition budget diagnostics - Part 1 of 2
!------------------------------------------------------------------------
IF ( State_Diag%Archive_BudgetEmisDryDep ) THEN
! Get initial column masses (full, trop, PBL)
CALL Compute_Budget_Diagnostics( &
Input_Opt = Input_Opt, &
State_Chm = State_Chm, &
State_Grid = State_Grid, &
State_Met = State_Met, &
isFull = State_Diag%Archive_BudgetEmisDryDepFull, &
diagFull = NULL(), &
mapDataFull = State_Diag%Map_BudgetEmisDryDepFull, &
isTrop = State_Diag%Archive_BudgetEmisDryDepTrop, &
diagTrop = NULL(), &
mapDataTrop = State_Diag%Map_BudgetEmisDryDepTrop, &
isPBL = State_Diag%Archive_BudgetEmisDryDepPBL, &
diagPBL = NULL(), &
mapDataPBL = State_Diag%Map_BudgetEmisDryDepPBL, &
colMass = State_Diag%BudgetColumnMass, &
before_op = .TRUE., &
RC = RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Emissions/dry deposition budget diagnostics error 1'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
ENDIF
#if defined( ADJOINT ) && defined ( DEBUG )
IF (Input_Opt%is_adjoint .and. Input_Opt%IS_FD_SPOT_THIS_PET) THEN
WRITE(*,*) ' SpcAdj(IFD,JFD) before unit converstion: ', &
State_Chm%SpeciesAdj(Input_Opt%IFD, Input_Opt%JFD, &
Input_Opt%LFD, Input_Opt%NFD)
WRITE(*,*) ' Spc(IFD,JFD) before unit converstion: ', &
State_Chm%Species(Input_Opt%NFD)%Conc(Input_Opt%IFD, Input_Opt%JFD, Input_Opt%LFD)
ENDIF
#endif
! DO_TEND previously operated in units of kg. The species arrays are in
! v/v for mixing, hence needed to convert before and after.
! Now use units kg/m2 as State_Chm%SPECIES units in DO_TEND to
! remove area-dependency (ewl, 9/30/15)
CALL Convert_Spc_Units( &
Input_Opt = Input_Opt, &
State_Chm = State_Chm, &
State_Grid = State_Grid, &
State_Met = State_Met, &
outUnit = KG_SPECIES_PER_M2, &
origUnit = origUnit, &
RC = RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Unit conversion error!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
#if defined( ADJOINT ) && defined ( DEBUG )
IF (Input_Opt%is_adjoint .and. Input_Opt%IS_FD_SPOT_THIS_PET) THEN
WRITE(*,*) ' SpcAdj(IFD,JFD) after unit converstion: ', &
State_Chm%SpeciesAdj(Input_Opt%IFD, Input_Opt%JFD, &
Input_Opt%LFD, Input_Opt%NFD)
WRITE(*,*) ' Spc(IFD,JFD) after unit converstion: ', &
State_Chm%Species(Input_Opt%NFD)%Conc(Input_Opt%IFD, Input_Opt%JFD, Input_Opt%LFD)
ENDIF
#endif
! Trap potential error
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Unit conversion error!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
! Get time step [s]
IF ( PRESENT(DT) ) THEN
TS = DT
ELSE
TS = GET_TS_DYN()
ENDIF
#ifdef ADJOINT
if (Input_Opt%Is_Adjoint) then
TS = TS * -1
endif
#endif
! First-time setup
IF ( FIRST ) THEN
! Define species indices on the first call
id_MACR = Ind_('MACR' )
id_RCHO = Ind_('RCHO' )
id_ACET = Ind_('ACET' )
id_ALD2 = Ind_('ALD2' )
id_ALK4 = Ind_('ALK4' )
id_C2H6 = Ind_('C2H6' )
id_C3H8 = Ind_('C3H8' )
id_CH2O = Ind_('CH2O' )
id_CO2 = Ind_('CO2' )
id_PRPE = Ind_('PRPE' )
id_O3 = Ind_('O3' )
id_HNO3 = Ind_('HNO3' )
id_BrO = Ind_('BrO' )
id_Br2 = Ind_('Br2' )
id_Br = Ind_('Br' )
id_HOBr = Ind_('HOBr' )
id_HBr = Ind_('HBr' )
id_BrNO3= Ind_('BrNO3')
id_CH4_SAB = Ind_('CH4_SAB')
FIRST = .FALSE.
ENDIF
! On first call, get pointers to the PARANOX loss fluxes. These are
! stored in diagnostics 'PARANOX_O3_DEPOSITION_FLUX' and
! 'PARANOX_HNO3_DEPOSITION_FLUX'. The call below links pointers
! PNOXLOSS_O3 and PNOXLOSS_HNO3 to the data values stored in the
! respective diagnostics. The pointers will remain unassociated if
! the diagnostics do not exist.
! This is only needed if non-local PBL scheme is not being used.
! Otherwise, PARANOX fluxes are applied in vdiff_mod.F90.
! (ckeller, 4/10/2015)
!
! If using HEMCO Intermediate grid feature, then the call needs to be
! refreshed at every time step for regridding. (hplin, 6/21/20)
IF ( .NOT. Input_Opt%LNLPBL ) THEN
CALL HCO_GC_GetDiagn( Input_Opt, State_Grid, 'PARANOX_O3_DEPOSITION_FLUX', &
.FALSE., RC, Ptr2D = Ptr2D )
IF( ASSOCIATED( Ptr2D )) THEN
ALLOCATE ( PNOxLoss_O3( State_Grid%NX, State_Grid%NY ), STAT=RC )
PNOxLoss_O3(:,:) = Ptr2D(:,:)
ENDIF
Ptr2D => NULL()
CALL HCO_GC_GetDiagn( Input_Opt, State_Grid, 'PARANOX_HNO3_DEPOSITION_FLUX',&
.FALSE., RC, Ptr2D = Ptr2D )
IF( ASSOCIATED( Ptr2D )) THEN
ALLOCATE ( PNOxLoss_HNO3( State_Grid%NX, State_Grid%NY ), STAT=RC )
PNOxLoss_HNO3(:,:) = Ptr2D(:,:)
ENDIF
Ptr2D => NULL()
ENDIF
!-----------------------------------------------------------------------
! For tagged CH4 simulations
! Save the total CH4 concentration before apply soil absorption
!-----------------------------------------------------------------------
IF ( Input_Opt%ITS_A_TAGCH4_SIM ) THEN
total_ch4_pre_soil_absorp = State_Chm%Species(1)%Conc(:,:,:)
ENDIF
!=======================================================================
! Do for every advected species and grid box
!=======================================================================
! Note: For GEOS-Chem Classic HEMCO "Intermediate Grid" feature,
! where HEMCO is running on a distinct grid from the model, the
! on-demand regridding is most optimized when contiguous accesses
! to GetHcoValEmis and GetHcoValDep are performed for the given species.
! Therefore, it is most optimal to call in the following fashion (IJKN)
! Emis(1,1,1,1) -> Emis(1,1,2,1) -> ... -> Dep(1,1,1,1) -> Dep(1,1,2,1)
! By switching emis/dep or the species # LAST, as either of these changing
! WILL trigger a new regrid and thrashing of the old buffer.
!
! Therefore, the loop below has been adjusted to run serially for each
! species, and parallelizing the inner I, J loop instead (hplin, 6/27/20)
! Also, moved some non-I,J specific variables outside of the loop for optimization
DO NA = 1, nAdvect
! Initialize PRIVATE error-handling variables
ErrorMsg = ''
! Get the species ID from the advected species ID
N = State_Chm%Map_Advect(NA)
! Get info about this species from the species database
SpcInfo => State_Chm%SpcData(N)%Info
! Molecular weight in kg
MWkg = SpcInfo%MW_g * 1.e-3_fp
!--------------------------------------------------------------------
! Check if we need to do dry deposition for this species
!--------------------------------------------------------------------
! Initialize
DryDepSpec = .FALSE.
DryDepID = -1
! Only if dry deposition is turned on and we do want to consider
! processes below the PBL...
IF ( LDRYD .AND. .NOT. OnlyAbovePBL ) THEN
! Get dry deposition ID (used by drydep_mod.F90) for this species.
! This is now stored in the species database object. (bmy, 7/6/16)
DryDepID = SpcInfo%DryDepId
! Check if this is a HEMCO drydep species
DryDepSpec = ( DryDepId > 0 )
IF ( .NOT. DryDepSpec ) THEN
CALL InquireHco ( N, Dep=DryDepSpec )
ENDIF
! Special case for O3 or HNO3: include PARANOX loss
IF ( N == id_O3 .AND. ASSOCIATED(PNOXLOSS_O3 ) ) &
DryDepSpec = .TRUE.
IF ( N == id_HNO3 .AND. ASSOCIATED(PNOXLOSS_HNO3) ) &
DryDepSpec = .TRUE.
ENDIF
! Set emissions top level:
! This is the top of atmosphere unless concentration build-up
! in stratosphere wants to be avoided.
ChemGridOnly = .FALSE.
! Set emissions to zero above chemistry grid for the following VOCs
IF ( N == id_MACR .OR. N == id_RCHO .OR. &
N == id_ACET .OR. N == id_ALD2 .OR. &
N == id_ALK4 .OR. N == id_C2H6 .OR. &
N == id_C3H8 .OR. N == id_CH2O .OR. &
N == id_PRPE ) THEN
ChemGridOnly = .TRUE.
ENDIF
! Bry concentrations become prescribed in lin. strat. chemistry.
! Therefore avoid any emissions of these compounds above the
! chemistry grid (lin. strat. chem. applies above chemistry grid
! only).
IF ( LINEAR_CHEM ) THEN
IF ( N == id_BrO .OR. N == id_Br2 .OR. &
N == id_Br .OR. N == id_HOBr .OR. &
N == id_HBr .OR. N == id_BrNO3 ) THEN
ChemGridOnly = .TRUE.
ENDIF
ENDIF
!--------------------------------------------------------------------
! Check if we need to do emissions for this species
!--------------------------------------------------------------------
IF ( LEMIS ) THEN
CALL InquireHco ( N, Emis=EmisSpec )
ELSE
EmisSpec = .FALSE.
ENDIF
! If there is emissions for this species, it must be loaded into memory first.
! This is achieved by attempting to retrieve a grid box while NOT in a parallel
! loop. Failure to load this will result in severe performance issues!! (hplin, 9/27/20)
IF ( EmisSpec ) THEN
CALL LoadHcoValEmis ( Input_Opt, State_Grid, N )
ENDIF
IF ( DryDepSpec ) THEN
CALL LoadHcoValDep ( Input_Opt, State_Grid, N )
ENDIF
!--------------------------------------------------------------------
! Can go to next species if this species does not have
! dry deposition and/or emissions
!--------------------------------------------------------------------
IF ( .NOT. DryDepSpec .AND. .NOT. EmisSpec ) CYCLE
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, L1, L2 ) &
!$OMP PRIVATE( PBL_TOP, FND, TMP ) &
!$OMP PRIVATE( FRQ, RKT, FRAC, FLUX, Area_m2 ) &
!$OMP PRIVATE( DRYD_TOP, EMIS_TOP, PNOXLOSS, DENOM ) &
!$OMP PRIVATE( S, ErrorMsg )
! Loop over all grid boxes
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
!-----------------------------------------------------------------
! Define various quantities before computing tendencies
!-----------------------------------------------------------------
! Get PBL_TOP at this grid box
PBL_TOP = MAX( 1, FLOOR( State_Met%PBL_TOP_L(I,J) ) )
! Determine lower level L1 to be used:
! If specified so, apply emissions only above the PBL_TOP.
! This will also disable dry deposition.
IF ( OnlyAbovePBL ) THEN
L1 = PBL_TOP + 1
ELSE
L1 = 1
ENDIF
! Set dry deposition top level based on PBL_DRYDEP flag of
! Input_Opt.
IF ( PBL_DRYDEP ) THEN
DRYD_TOP = PBL_TOP
ELSE
DRYD_TOP = 1
ENDIF
! Restrict to chemistry grid
IF ( ChemGridOnly ) THEN
EMIS_TOP = State_Met%ChemGridLev(I,J)
EMIS_TOP = MIN(State_Grid%NZ,EMIS_TOP)
ELSE
EMIS_TOP = State_Grid%NZ
ENDIF
! L2 is the upper level index to loop over
L2 = MAX(DRYD_TOP, EMIS_TOP)
! This should not happen:
IF ( L2 < L1 ) CYCLE
! Loop over selected vertical levels
DO L = L1, L2
!--------------------------------------------------------------
! Apply dry deposition frequencies to all levels below the
! PBL top.
!--------------------------------------------------------------
IF ( DryDepSpec .AND. ( L <= DRYD_TOP ) ) THEN
! Init
FRQ = 0.0_fp
! Dry deposition frequency from drydep_mod.F90. This is
! stored in State_Chm%DryDepFreq. Units are [s-1].
IF ( DRYDEPID > 0 ) THEN
FRQ = DepFreq(I,J,DRYDEPID)
ENDIF
! Dry deposition frequency from HEMCO. HEMCO calculates
! dry deposition frequencies for air-sea exchange and
! from ship NOx plume parameterization (PARANOx). The
! units are [s-1].
CALL GetHcoValDep ( Input_Opt, State_Grid, N, I, J, 1, FND, TMP )
! Add to dry dep frequency from drydep_mod.F90
IF ( FND ) FRQ = FRQ + TMP
! Get PARANOX deposition loss. Apply to surface level only.
! PNOXLOSS is in kg/m2/s. (ckeller, 4/10/15)
PNOXLOSS = 0.0_fp
IF ( L == 1 ) THEN
IF ( N == id_O3 .AND. ASSOCIATED(PNOXLOSS_O3) ) THEN
PNOXLOSS = PNOXLOSS_O3(I,J)
ENDIF
IF ( N == id_HNO3 .AND. ASSOCIATED(PNOXLOSS_HNO3) ) THEN
PNOXLOSS = PNOXLOSS_HNO3(I,J)
ENDIF
ENDIF
! Apply dry deposition
IF ( FRQ > 0.0_fp .OR. PNOXLOSS > 0.0_fp ) THEN
! Compute exponential loss term
RKT = FRQ * TS
FRAC = EXP(-RKT)
! Loss in kg/m2
FLUX = ( 1.0_fp - FRAC ) * State_Chm%Species(N)%Conc(I,J,L)
! Apply dry deposition
State_Chm%Species(N)%Conc(I,J,L) = FRAC * &
State_Chm%Species(N)%Conc(I,J,L)
#ifdef ADJOINT
if (Input_Opt%Is_Adjoint) then
State_Chm%SpeciesAdj(I,J,L,N) = FRAC * &
State_Chm%SpeciesAdj(I,J,L,N)
endif
#endif
! Eventually add PARANOX loss. PNOXLOSS is in kg/m2/s.
! Make sure PARANOx loss is applied to tracers. (ckeller,
! 3/29/16).
IF ( PNOXLOSS > 0 ) THEN
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) - ( PNOXLOSS * TS )
FLUX = FLUX + ( PNOXLOSS * TS )
ENDIF
! Loss in [molec/cm2/s]
! Added a safe_div due to small parallelization error
! (mdy, 5/15)
!
! NOTE: The original computation was:
! FLUX = FLUX / MWkg * AVO / TS / ( AREA_M2 * 1.0e4_fp ) ]
! so the denominator as we had it was wrong.
! Now corrected (elundgren, bmy, 6/12/15)
DENOM = ( MWkg * TS * 1.0e+4_fp ) / AVO
FLUX = SAFE_DIV( FLUX, DENOM, 0.0e+0_fp ) ! molec/cm2/s
! Eventually add to SOIL_DRYDEP
IF ( Input_Opt%LSOILNOX ) THEN
CALL SOIL_DRYDEP( I, J, L, N, FLUX, State_Chm )
ENDIF
!--------------------------------------------------------
! HISTORY: Archive drydep flux loss from mixing
! Units = molec/cm2/s
!
! NOTE: we don't need to multiply by the ratio of
! TS_CONV / TS_CHEM, as the updating frequency for
! HISTORY is determined by the "frequency" setting in
! the "HISTORY.rc"input file. The old bpch diagnostics
! archived the drydep due to chemistry every chemistry
! timestep = 2X the dynamic timestep. So in order to
! avoid double-counting the drydep flux from mixing,
! you had to multiply by TS_CONV / TS_CHEM.
!
! ALSO NOTE: When comparing History output to bpch
! output, you must use an updating frequency equal to
! the dynamic timestep so that the drydep fluxes due to
! mixing will be equivalent w/ the bpch output. It is
! also recommended to turn off chemistry so as to be
! able to compare the drydep fluxes due to mixing in
! bpch vs. History as an "apples-to-apples" comparison.
!
! -- Bob Yantosca (yantosca@seas.harvard.edu)
!--------------------------------------------------------
IF ( ( State_Diag%Archive_DryDepMix .or. &
State_Diag%Archive_DryDep ) .and. &
DryDepID > 0 ) THEN
S = State_Diag%Map_DryDepMix%id2slot(DryDepID)
IF ( S > 0 ) THEN
State_Diag%DryDepMix(I,J,S) = Flux
ENDIF
ENDIF
ENDIF ! apply drydep
ENDIF ! L <= PBLTOP
!--------------------------------------------------------------
! Apply emissions.
! These are always taken from HEMCO
!--------------------------------------------------------------
IF ( EmisSpec .AND. ( L <= EMIS_TOP ) ) THEN
! Get HEMCO emissions. Units are [kg/m2/s].
! Fix hplin: for intermediate grid, pass SkipCheck in a tight loop. Note that this assumes that adjacent
! calls to GetHcoValEmis are from the same species ID, or there will be big trouble. (hplin, 10/10/20)
#ifdef MODEL_CLASSIC
IF ( Input_Opt%LIMGRID ) THEN
FND = .true.
TMP = TMP_MDL(I,J,L) ! this is a kludge for the tight loop optimization
ELSE
#endif
CALL GetHcoValEmis ( Input_Opt, State_Grid, N, I, J, L, FND, TMP, SkipCheck=.true. )
#ifdef MODEL_CLASSIC
ENDIF
#endif
! Add emissions (if any)
! Bug fix: allow negative fluxes. (ckeller, 4/12/17)
!IF ( FND .AND. (TMP > 0.0_fp) ) THEN
IF ( FND ) THEN
! Flux: [kg/m2] = [kg m-2 s-1 ] x [s]
FLUX = TMP * TS
#ifdef ADJOINT
IF ( I .eq. Input_Opt%IFD .and. J .eq. Input_Opt%JFD .and. &
L .eq. Input_Opt%LFD .and. N .eq. Input_Opt%NFD) THEN
WRITE(*,*) ' GetHcoVal(IFD,JFD) = ', TMP, ' FLUX = ', FLUX
IF ( Input_Opt%is_adjoint ) THEN
WRITE(*,*) ' SpeciesAdj(FD) = ', State_Chm%SpeciesAdj(I,J,L,N)
ENDIF
ENDIF
#endif
! Add to species array
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) + FLUX
ENDIF
ENDIF
!--------------------------------------------------------------
! Special handling for tagged CH4 simulations
!
! Tagged CH4 species are split off into a separate loop to
! ensure we remove soil absorption from NA=1 (total CH4) first
!--------------------------------------------------------------
IF ( Input_Opt%ITS_A_TAGCH4_SIM ) THEN
! If we are in the chemistry grid
IF ( L <= EMIS_TOP ) THEN
! Total CH4 species
IF ( NA == 1 ) THEN
! Get soil absorption from HEMCO. Units are [kg/m2/s].
! CH4_SAB is species #15
CALL GetHcoValEmis ( Input_Opt, State_Grid, 15, I, J, L, FND, TMP )
! Remove soil absorption from total CH4 emissions
IF ( FND ) THEN
! Flux: [kg/m2] = [kg m-2 s-1 ] x [s]
FLUX = TMP * TS
! Apply soil absorption as loss
State_Chm%Species(N)%Conc(I,J,L) = &
State_Chm%Species(N)%Conc(I,J,L) - FLUX
ENDIF
ENDIF
ENDIF
ENDIF
! Check for negative concentrations
IF ( State_Chm%Species(N)%Conc(I,J,L) < 0.0_fp ) THEN
#ifdef TOMAS
! For TOMAS simulations only, look for negative and reset
! to small positive. This prevents the run from dying,
! while we look for the root cause of the issue.
! -- Betty Croft, Bob Yantosca (21 Jan 2022)
print *, 'Found negative ', N, State_Chm%Species(N)%Conc(I,J,L)
State_Chm%Species(N)%Conc(I,J,L) = 1e-26_fp
#else
IF ( ( N /= id_CH4_SAB ) .and. ( N /= id_CO2 ) ) THEN
! KLUDGE: skip the warning message for CH4_SAB, which can be
! negative (it's a soil absorption flux). The TagCH4
! simulation is not used regularly as of Feb 2021 -- fix this
! later if need by. (bmy, 2/25/21)
Print*, 'WARNING: Negative concentration for species ', &
TRIM( SpcInfo%Name), ' at (I,J,L) = ', I, J, L
ErrorMsg = 'Negative species concentations encountered.' // &
' This may be fixed by increasing the' // &
' background concentration or by shortening' // &
' the transport time step.'
RC = GC_FAILURE
ENDIF
#endif
ENDIF
ENDDO !L
ENDDO !J
ENDDO !I
!$OMP END PARALLEL DO
! Exit with error condition
IF ( RC /= GC_SUCCESS ) THEN
CALL GC_Error( ErrorMsg, RC, ThisLoc )
RETURN
ENDIF
! Nullify pointer
SpcInfo => NULL()
ENDDO !N
!--------------------------------------------------------------
! Special handling for tagged CH4 simulations
!--------------------------------------------------------------
IF ( Input_Opt%ITS_A_TAGCH4_SIM ) THEN
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, N, NA, ErrorMsg )
DO NA = 1, nAdvect
! Get the species ID from the advected species ID
N = State_Chm%Map_Advect(NA)
! Loop over all grid boxes
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX
! Tagged CH4 tracers
IF ( NA >= 2 .and. NA <= nAdvect-1 ) THEN
! Apply soil absorption to each tagged CH4 species
! (Xueying Yu, 12/08/2017)
State_Chm%Species(N)%Conc(I,J,L) = &
SAFE_DIV(State_Chm%Species(N)%Conc(I,J,L), &
total_ch4_pre_soil_absorp(I,J,L), &
0.e+0_fp) * &
State_Chm%Species(1)%Conc(I,J,L)
ENDIF
! Check for negative concentrations
! KLUDGE: skip the warning message for CH4_SAB, which can be
! negative (it's a soil absorption flux). The TagCH4 simulation
! is not used regularly as of Feb 2021 -- fix this later if
! need by. (bmy, 2/25/21)
IF ( State_Chm%Species(N)%Conc(I,J,L) < 0.0_fp ) THEN
IF ( N /= id_CH4_SAB ) THEN
Print*, 'WARNING: Negative concentration for species ', &
TRIM( State_Chm%SpcData(N)%Info%Name), &
' at (I,J,L) = ', I, J, L
ErrorMsg = 'Negative species concentations encountered.' // &
' This may be fixed by increasing the' // &
' background concentration or by shortening' // &
' the transport time step.'
RC = GC_FAILURE
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Exit with error condition
IF ( RC /= GC_SUCCESS ) THEN
CALL GC_Error( ErrorMsg, RC, ThisLoc )
RETURN
ENDIF
ENDIF
#if defined( ADJOINT ) && defined ( DEBUG )
IF (Input_Opt%is_adjoint .and. Input_Opt%IS_FD_SPOT_THIS_PET) THEN
WRITE(*,*) ' SpcAdj(IFD,JFD) before unit converstion: ', &
State_Chm%SpeciesAdj(Input_Opt%IFD, Input_Opt%JFD, &
Input_Opt%LFD, Input_Opt%NFD)
WRITE(*,*) ' Spc(IFD,JFD) before unit converstion: ', &
State_Chm%Species(Input_Opt%NFD)%Conc(Input_Opt%IFD, Input_Opt%JFD, Input_Opt%LFD)
ENDIF
#endif
! Convert State_Chm%Species back to original units
CALL Convert_Spc_Units( &
Input_Opt = Input_Opt, &
State_Chm = State_Chm, &
State_Grid = State_Grid, &
State_Met = State_Met, &
outUnit = origUnit, &
RC = RC )
IF ( RC /= GC_SUCCESS ) THEN
ErrMsg = 'Unit conversion error!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
#if defined( ADJOINT ) && defined ( DEBUG )
IF (Input_Opt%is_adjoint .and. Input_Opt%IS_FD_SPOT_THIS_PET) THEN
WRITE(*,*) ' SpcAdj(IFD,JFD) after unit converstion: ', &
State_Chm%SpeciesAdj(Input_Opt%IFD, Input_Opt%JFD, &
Input_Opt%LFD, Input_Opt%NFD)
WRITE(*,*) ' Spc(IFD,JFD) after unit converstion: ', &
State_Chm%Species(Input_Opt%NFD)%Conc(Input_Opt%IFD, Input_Opt%JFD, Input_Opt%LFD)
ENDIF
#endif
!------------------------------------------------------------------------
! Emissions/dry deposition budget diagnostics - Part 2 of 2
!------------------------------------------------------------------------
IF ( State_Diag%Archive_BudgetEmisDryDep ) THEN
! Timestep for diagnostics [s]
DT_Tend = DBLE( TS )
! Compute change in column masses (after emis/dryd - before emis/dryd)
! and store in diagnostic arrays. Units are [kg/s].
CALL Compute_Budget_Diagnostics( &
Input_Opt = Input_Opt, &
State_Chm = State_Chm, &
State_Grid = State_Grid, &
State_Met = State_Met, &
isFull = State_Diag%Archive_BudgetEmisDryDepFull, &
diagFull = State_Diag%BudgetEmisDryDepFull, &
mapDataFull = State_Diag%Map_BudgetEmisDryDepFull, &
isTrop = State_Diag%Archive_BudgetEmisDryDepTrop, &
diagTrop = State_Diag%BudgetEmisDryDepTrop, &
mapDataTrop = State_Diag%Map_BudgetEmisDryDepTrop, &
isPBL = State_Diag%Archive_BudgetEmisDryDepPBL, &
diagPBL = State_Diag%BudgetEmisDryDepPBL, &
mapDataPBL = State_Diag%Map_BudgetEmisDryDepPBL, &
colMass = State_Diag%BudgetColumnMass, &
timeStep = DT_Tend, &
RC = RC )