-
Notifications
You must be signed in to change notification settings - Fork 174
/
Copy pathtomas_mod.F90
10699 lines (9746 loc) · 388 KB
/
tomas_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
#ifdef TOMAS
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: tomas_mod.F90
!
! !DESCRIPTION: Module TOMAS\_MOD contains variable specific to the TOMAS
! aerosol microphysics simulation, e.g. number of species, number of size-bins,
! mass per particle bin boundaries and arrays used inside the microphysics
! subroutine. (win, 7/9/07)
!\\
!\\
! !INTERFACE:
!
MODULE TOMAS_MOD
!
! !USES:
!
USE PRECISION_MOD ! For GEOS-Chem Precision (fp)
IMPLICIT NONE
!
! !REMARKS:
! This module also contains what used to be in sizecode.COM header file
! containing common blocks for TOMAS aerosol microphysics algorithm
! originally implemented in GISS GCM-II' by Peter Adams. Below is the original
! comment from sizecode.COM
!
! This header file includes all the variables used by the
! size-resolved aerosol microphysics code incorporated into
! the GISS GCM II' by Peter Adams. The microphysics algorithm
! conserves aerosol number and mass using the schemes developed
! by Graham Feingold and others.
!
! References:
! ============================================================================
! Tzivion, S., Feingold, G., and Levin, Z., An Efficient
! Numerical Solution to the Stochastic Collection Equation,
! J. Atmos. Sci., 44, 3139-3149, 1987.
! Feingold, G., Tzivion, S., and Levin, Z., Evolution of
! Raindrop Spectra. Part I: Solution to the Stochastic
! Collection/Breakup Equation Using the Method of Moments,
! J. Atmos. Sci., 45, 3387-3399, 1988.
! Tzivion, S., Feingold, G., and Levin, Z., The Evolution of
! Raindrop Spectra. Part II: Collisional Collection/Breakup
! and Evaporation in a Rainshaft, J. Atmos. Sci., 46, 3312-
! 3327, 1989.
! Feingold, G., Levin, Z., and Tzivion, S., The Evolution of
! Raindrop Spectra. Part III: Downdraft Generation in an
! Axisymmetrical Rainshaft Model, J. Atmos. Sci., 48, 315-
! 330, 1991.
!
! The algorithms described in these papers have been extended
! to include multicomponent aerosols and modified for a moving
! sectional approach. Using this approach, the boundaries
! between size bins are defined in terms of dry aerosol mass
! such that the actual sizes of the sections move as water
! is added to or lost from the aerosol.
!
! All of the subroutines needed for this aerosol microphysics
! algorithm use only their own internal variables or the ones
! listed here. GISS GCM II' variables are not used (a driver
! subroutine performs the necessary swapping between the GCM
! and the microphysics code). The microphysics code is,
! therefore, completely modular.
!
! !REVISION HISTORY:
! 09 Jul 2006 - W. Trivitayanurak - Initial version
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!========================================================================
! Module Variables:
! ICOMP : Number of aerosol mass species + 1 for number
! Nk : Aerosol number internal array
! Mk : Aerosol mass internal array
! Gc : Condensing gas
! Nkd : Aerosol number diagnostic array
! Mkd : Aerosol mass diagnostic array
! Gcd : Condensing gas diagnostic array
! Xk : Size bin boundary in dry mass per particle
! MOLWT : Aerosol molecular weight
! SRTSO4 : ID of sulfate
! SRTNACL : ID of sea-salt
! SRTH2O : ID of aerosol water
! SRTECOB : ID of hydrophobic EC
! SRTECIL : ID of hydrophilic EC
! SRTOCOB : ID of hydrophobic OC
! SRTOCIL : ID of hydrophilic OC
! dust?? : ID of internally mixed dust (future work)
! dust?? : ID of externally mixed dust (future work)
! BOXVOL : Grid box volume (cm3)
! BOXMASS : Grid box air mass (kg)
! TEMPTMS : Temperature (K) of each grid box
! PRES : Pressure (Pa) in grid box
! RHTOMAS : Relative humidity (0-1)
! BINACT1 : Activated bin as a function of composition
! FRACTION: Activated fraction as a fcn of composition
! IDIAG : Number of diagnostic tracer (NH4 and H2O)
!========================================================================
!
! !DEFINED PARAMETERS:
!
INTEGER, PARAMETER :: SRTSO4 = 1
INTEGER, PARAMETER :: SRTNACL = 2
INTEGER, PARAMETER :: SRTECIL = 3
INTEGER, PARAMETER :: SRTECOB = 4
INTEGER, PARAMETER :: SRTOCIL = 5
INTEGER, PARAMETER :: SRTOCOB = 6
INTEGER, PARAMETER :: SRTDUST = 7
INTEGER, PARAMETER :: SRTNH4 = 8
INTEGER, PARAMETER :: SRTH2O = 9
INTEGER, PARAMETER :: ICOMPHARD = 9
!
! !PUBLIC DATA MEMBERS:
!
! Scalars
INTEGER :: ICOMP, IDIAG
! Arrays
REAL(fp), SAVE, ALLOCATABLE, TARGET :: Xk(:)
REAL*4, SAVE, ALLOCATABLE :: MOLWT(:)
INTEGER, SAVE :: BINACT1(101,101,101)
REAL(fp), SAVE :: FRACTION1(101,101,101)
INTEGER, SAVE :: BINACT2(101,101,101)
REAL(fp), SAVE :: FRACTION2(101,101,101)
REAL(fp), ALLOCATABLE :: AVGMASS(:) ! Average mass per particle
! mid-range of size bin
! [kg/no.]
REAL(fp) :: cosmic_ions(72,46,9) !careful, this is not GCHP safe!
! [ion pairs kg^-1 s^-1]
REAL(fp), SAVE, ALLOCATABLE :: OCSCALE30(:)
REAL(fp), SAVE, ALLOCATABLE :: OCSCALE100(:)
REAL(fp), SAVE, ALLOCATABLE :: ECSCALE30(:)
REAL(fp), SAVE, ALLOCATABLE :: ECSCALE100(:)
INTEGER :: bin_nuc = 1, tern_nuc = 1 ! Switches for nucleation type.
INTEGER :: act_nuc = 0 ! in BL
INTEGER :: ion_nuc = 0 ! 1 for modgil, 2 for Yu
! (Yu currently broken, JRP 202101)
INTEGER :: absall = 1 ! 1 for soa absorb to all specnapari
! nucleation tuned by factor of 1.0D-5
REAL(fp) :: soaareafrac=1.0e+0_fp ! fraction of SOA that goes
! to SA rather than into mass
INTEGER :: xSOA = 1 !Switch for xSOA. If set to one, emit 100
! Tg/yr SOA correlated with anrtho CO
! (JKodros 6/3/15)
! switch to 0 for anthro-free
! runs (Pengfei Liu 4/18/2018)
INTEGER :: lowRH = 1 !This is to match AW more with AERONET (JKODROS 6/15)
REAL(fp), ALLOCATABLE :: H2SO4_RATE(:,:,:) ! H2SO4 prod rate [kg s-1]
REAL(fp), ALLOCATABLE :: PSO4AQ_RATE(:,:,:) ! Cld chem sulfate prod rate [kg s-1]
! Subgrid coagulation timescale (win, 10/28/08)
REAL*4 :: SGCTSCALE
!
! !PRIVATE TYPES:
!
! Number of bins (copied from State_Chm%nTomasBins)
INTEGER, PRIVATE, SAVE :: IBINS
! Species ID flags
INTEGER, PRIVATE :: id_AW01
INTEGER, PRIVATE :: id_DUST01
INTEGER, PRIVATE :: id_ECIL01
INTEGER, PRIVATE :: id_ECOB01
INTEGER, PRIVATE :: id_H2SO4
INTEGER, PRIVATE :: id_NH3
INTEGER, PRIVATE :: id_NH4
INTEGER, PRIVATE :: id_NK01
INTEGER, PRIVATE :: id_OCIL01
INTEGER, PRIVATE :: id_OCOB01
INTEGER, PRIVATE :: id_SF01
INTEGER, PRIVATE :: id_SO4
INTEGER, PRIVATE :: id_SS01
CONTAINS
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: do_tomas
!
! !DESCRIPTION: Subroutine DO\_TOMAS is the driver subroutine that calls the
! appropriate aerosol microphysics and dry deposition subroutines. (win,
! 7/23/07)
!\\
!\\
! !INTERFACE:
!
SUBROUTINE DO_TOMAS( Input_Opt, State_Chm, State_Diag, State_Grid, &
State_Met, RC )
!
! !USES:
!
USE ErrCode_Mod
USE ERROR_MOD
USE Input_Opt_Mod, ONLY : OptInput
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 UnitConv_Mod
!
! !INPUT PARAMETERS:
!
TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object
TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object
TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object
TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diagnostics State object
!
! !OUTPUT PARAMETERS:
!
INTEGER, INTENT(OUT) :: RC ! Success or failure?
!
! !REVISION HISTORY:
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
! Scalars
INTEGER :: N
! Strings
CHARACTER(LEN=255) :: MSG, LOC
!=================================================================
! DO_TOMAS begins here
!=================================================================
! Assume success
RC = GC_SUCCESS
! Check that species units are in [kg] (ewl, 8/13/15)
IF ( State_Chm%Spc_Units /= KG_SPECIES ) THEN
MSG = 'Incorrect species units: ' // TRIM(UNIT_STR(State_Chm%Spc_Units))
LOC = 'Routine DO_TOMAS in tomas_mod.F90'
CALL GC_Error( MSG, RC, LOC )
ENDIF
! Do TOMAS aerosol microphysics
CALL AEROPHYS( Input_Opt, State_Chm, State_Grid, State_Met, &
State_Diag, RC )
!print *, 'call checkmn in tomas_mod:222'
CALL CHECKMN( 0, 0, 0, Input_Opt, State_Chm, State_Grid, &
State_Met, State_Diag, 'Before Aerodrydep', RC)
! in kg
! Do dry deposition
IF ( Input_Opt%LDRYD ) THEN
CALL AERO_DRYDEP( Input_Opt, State_Chm, State_Diag, &
State_Grid, State_Met, RC )
ENDIF
! not in kg
!print *, 'call checkmn in tomas_mod:229'
CALL CHECKMN( 0, 0, 0, Input_Opt, State_Chm, State_Grid, &
State_Met, State_Diag, 'Before exiting DO_TOMAS', RC )
END SUBROUTINE DO_TOMAS
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: aerophys
!
! !DESCRIPTION: Subroutine AEROPHYS does aerosol microphysics, including
! nucleation, coagulation, and condensation.
!\\
!\\
! !INTERFACE:
!
SUBROUTINE AEROPHYS( Input_Opt, State_Chm, State_Grid, State_Met, &
State_Diag, RC )
!
! !USES:
!
#ifdef BPCH_DIAG
USE CMN_DIAG_MOD
USE DIAG_MOD, ONLY : AD61, AD61_INST
#endif
USE ErrCode_Mod
USE ERROR_MOD
USE Input_Opt_Mod, ONLY : OptInput
USE Species_Mod, ONLY : SpcConc
USE State_Chm_Mod, ONLY : ChmState
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
USE State_Diag_Mod, ONLY : DgnState
USE TIME_MOD, ONLY : GET_TS_CHEM
USE UnitConv_Mod
!
! !INPUT PARAMETERS:
!
TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object
TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object
TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object
!
! !INPUT/OUTPUT PARAMETERS:
!
TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object
TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diagnostics State object
!
! !OUTPUT PARAMETERS:
!
INTEGER, INTENT(OUT) :: RC ! Success or failure?
!
! !REVISION HISTORY:
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES
!
INTEGER :: I, J, L, N, JC, K !counters
INTEGER :: MPNUM !microphysical process id #
REAL*4 :: ADT !aerosol microphysics time step (seconds)
REAL(fp) :: QSAT !used in RH calculation
INTEGER :: TRACNUM
REAL(fp) :: FRAC
CHARACTER(LEN=255) :: MSG, LOC ! species unit check
! Arguments for CHECK_VALUE; avoids array temporaries (bmy, 1/28/14)
CHARACTER(LEN=255) :: ERR_VAR
CHARACTER(LEN=255) :: ERR_MSG
INTEGER :: ERR_IND(4)
!---------
!sfarina - move definitions of these from module common
!to within each sub and pass them around for openmp
REAL(fp) :: Nk(IBINS), Nkd(IBINS)
REAL(fp) :: Mk(IBINS, ICOMP)
REAL(fp) :: Mkd(IBINS,ICOMP)
REAL(fp) :: Gc(ICOMP - 1)
REAL(fp) :: Gcd(ICOMP - 1)
REAL*4 :: BOXVOL, BOXMASS, TEMPTMS
REAL*4 :: PRES, RHTOMAS
REAL(fp) :: surf_area ! aerosol surface area [micon^2 cm^-3]
REAL(fp) :: ionrate ! ion pair formation
! rate [ion pairs cm^-3 s^-1]
!---------
REAL(fp) :: Nkout(ibins), Mkout(ibins,icomp)
REAL(fp) :: Gcout(icomp-1)
REAL(fp) :: Nknuc(ibins), Mknuc(ibins,icomp)
REAL(fp) :: Nkcond(ibins), Mkcond(ibins,icomp)
REAL(fp) :: fn ! nucleation rate of clusters cm-3 s-1
REAL(fp) :: fn1 ! formation rate of particles to first size bin cm-3 s-1
REAL(fp) :: nucrate(State_Grid%NY,State_Grid%NZ)
REAL(fp) :: nucrate1(State_Grid%NY,State_Grid%NZ)
REAL(fp) :: tot_n_1, tot_n_1a, tot_n_2, tot_n_i ! used for nitrogen mass checks
REAL(fp) :: tot_s_1, tot_s_1a, tot_s_2 ! used for sulfur mass checks
REAL(fp) :: h2so4rate_o ! H2SO4rate for the specific grid cell
REAL(fp) :: TOT_Mk, TOT_nk ! for checking mass and number
REAL(fp) :: transfer(ibins)
LOGICAL :: PRINTNEG !<step4.0-temp> (win, 3/24/05)
logical :: ERRORSWITCH !<step4.2> To see where mnfix found negative value (win, 9/12/05)
logical :: errspot !<step4.4> To see where so4cond found errors (win, 9/21/05)
logical :: printdebug !<step4.3> Print out for debugging (win, 9/16/05)
logical :: COND, COAG, NUCL !<step5.1> switch for each process (win 4/8/06)
integer :: iob, job,lob !Just declare in case I want to debug (4/8/06)
data iob, job, lob / 1 , 1 , 1 /
real(fp) :: NH3_to_NH4, CEPS
parameter ( CEPS=1.e-17_fp )
real(fp) igR
parameter (igR=8.314) !Ideal gas constant J/mol.K
!The following are constants used in calculating rel. humidity
real(fp) axcons, bxcons, bytf, tf !for RH calculation
parameter(axcons=1.8094520287589733, &
bxcons=0.0021672473136556273, &
bytf=0.0036608580560606877, tf=273.16)
!lhe and lhs are the latent heats of evaporation and sublimation
logical, save :: firsttime = .true.
integer :: num_iter
real(fp) cplevels(9) ! cosmic ray pressure levels (for interp)
data cplevels / 959., 894., 786., 634., 468., &
321., 201., 103., 27. /
integer lev
real(fp) weight
double precision soil_ions(9) ! ion pairs cm-3 s-1 from radioactive elements in soil
data soil_ions / 5.,3.3,1.8,0.7,0.,0.,0.,0.,0./
! Pointers
TYPE(SpcConc), POINTER :: Spc(:)
!=================================================================
! AEROPHYS begins here
!=================================================================
! Assume success
RC = GC_SUCCESS
! Check that species units are in [kg]. While species units
! are now generally [kg/kg] in GEOS-Chem, they are converted to
! kg for TOMAS elsewhere in tomas_mod prior to calling this subroutine
! (ewl, 8/13/15)
IF ( State_Chm%Spc_Units /= KG_SPECIES ) THEN
MSG = 'Incorrect species units: ' // TRIM(UNIT_STR(State_Chm%Spc_Units))
LOC = 'Routine AEROPHYS in tomas_mod.F90'
CALL GC_Error( MSG, RC, LOC )
ENDIF
! Point to chemical species array [kg]
Spc => State_Chm%Species
! Initialize debugging and error-signal switches
printneg = .FALSE.
errorswitch = .FALSE.
PRINTDEBUG = .FALSE.
ERRSPOT = .FALSE.
! Initialize switches for each microphysical process
COND = .TRUE.
COAG = .TRUE.
NUCL = .TRUE.
! Initialize nucleation rate arrays
nucrate(:,:) = 0.e+0_fp
nucrate1(:,:) = 0.e+0_fp
! First-time setup
if (firsttime) then
!====================================================================
! Make sure there is access to the H2SO4 production rate array
! H2SO4RATE, which saves the H2SO4 production rate for EACH chemistry
! timestep. The prod/loss family has to be set with at least one
! with the family name PSO4 and SO4 as its one member. (win, 9/30/08)
!====================================================================
write(*,*) 'AEROPHYS: This run uses coupled condensation-', &
'nucleation scheme with pseudo-steady state H2SO4'
if(tern_nuc == 1) then
write(*,*)' Nucleation: Ternary (Napari ', &
'et al. 2002) and Binary (Vehkamaki et al. 2002)'
else
write(*,*)' Nucleation: Binary (Vehkamaki et al. 2002)'
endif
firsttime = .false.
endif
! Get chemistry timestep for use below [s]
! NOTE: This doesn't have to be !$OMP+PRIVATE (bmy, 2/7/20)
ADT = GET_TS_CHEM()
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%% NOTE: THIS PARALLEL LOOP MAY BE ABLE TO BE REVERSED TO L-J-I
!%%% WHICH IS MUCH MORE EFFICIENT (bmy, 1/28/14)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L ) &
!$OMP PRIVATE( PRES, TEMPTMS, BOXMASS, RHTOMAS, BOXVOL ) &
!$OMP PRIVATE( printneg, ionrate, lev, weight, GC, N, NK, JC ) &
!$OMP PRIVATE( MK, H2SO4rate_o, tot_n_1, k, tot_s_1, MPNUM ) &
!$OMP PRIVATE( Nkd, Mkd, TOT_NK, TOT_MK, TRANSFER ) &
!$OMP PRIVATE( Nkout,Mkout,Gcout,fn,fn1 ) &
!$OMP PRIVATE( num_iter,Nknuc,Mknuc,Nkcond ) &
!$OMP PRIVATE( Mkcond, ERRORSWITCH, tot_s_1a, tot_n_1a ) &
!$OMP PRIVATE( ERR_VAR, ERR_MSG, ERR_IND ) &
!$OMP PRIVATE( TRACNUM, NH3_TO_NH4, SURF_AREA ) &
!$OMP SCHEDULE( DYNAMIC )
DO I = 1, State_Grid%NX
DO J = 1, State_Grid%NY
DO L = 1, State_Grid%NZ
#ifdef BPCH_DIAG
! Reset the AD61_INST array used for tracking instantaneous
! certain rates. As of now, tracking nucleation (win, 10/7/08)
IF ( Input_Opt%ND61 > 0 ) THEN
AD61_INST(I,J,L,:) = 0e0
ENDIF
#endif
! Skip non-chemgrid boxes
IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE
!vbn write(890,89)I,J,L,Spc(id_H2SO4)%Conc(i,j,l)
!if(printdebug) print *,'+++++',I,J,L,'inside Aerophys'
! Get info on this grid box
PRES = State_Met%PMID(i,j,l)*100.0 ! in Pa
TEMPTMS = State_Met%T(I,J,L)
BOXMASS = State_Met%AD(I,J,L)
RHTOMAS = State_Met%RH(I,J,L)/ 1.e2
IF ( RHTOMAS > 0.99 ) RHTOMAS = 0.99
BOXVOL = State_Met%AIRVOL(I,J,L) * 1.e6 !convert from m3 -> cm3
printneg = .FALSE.
! determine ion rate
!ionrate = 10. ! set as constant now !!jrp
IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
if( (pres / 100.) .lt. cplevels(9) ) then
ionrate = cosmic_ions(i,j,9) * boxmass / boxvol
if(State_Met%FRCLND(I,J) .gt. 0.2) then
ionrate = ionrate + soil_ions(9)
endif
elseif((pres/100.) .gt. cplevels(1)) then
ionrate = cosmic_ions(i,j,1) * boxmass / boxvol
if( State_Met%FRCLND(I,J) .gt. 0.2 ) then
ionrate = ionrate + soil_ions(1)
endif
else
lev=2
do while (pres / 100. .lt. cplevels(lev))
lev=lev+1
enddo
weight=( cplevels( lev - 1 ) - pres / 100. ) / &
( cplevels( lev - 1 ) - cplevels(lev) )
ionrate=( cosmic_ions(i,j,lev ) * weight + &
cosmic_ions(i,j,lev-1) * (1.e+0_fp - weight) ) &
* boxmass / boxvol
if( State_Met%FRCLND(I,J) .gt. 0.2) then
ionrate=ionrate + ( soil_ions( lev ) * weight + &
soil_ions( lev-1 ) * (1.e+0_fp-weight) )
endif
endif
ELSE
ionrate = 0.e+0_fp
ENDIF
if(ionrate .le. 1.501) ionrate = 1.501
!print*,'i',i,'j',j,'l',l,'ionrate',ionrate
! Initialize all condensible gas values to zero
! Gc(srtso4) will remain zero until within cond_nuc where the
! pseudo steady state H2SO4 concentration will be put in this place.
DO JC=1, ICOMP-1
Gc(JC) = 0.e+0_fp
ENDDO
! Swap Spc into Nk, Mk, Gc arrays
DO N = 1, IBINS
NK(N) = Spc(id_NK01-1+N)%Conc(I,J,L)
DO JC = 1, ICOMP-IDIAG
MK(N,JC) = Spc(id_NK01-1+N+JC*IBINS)%Conc(I,J,L)
IF( IT_IS_NAN( MK(N,JC) ) ) THEN
PRINT *,'+++++++ Found NaN in AEROPHYS ++++++++'
PRINT *,'Location (I,J,L):',I,J,L,'Bin',N,'comp',JC
ENDIF
ENDDO
MK(N,SRTH2O) = Spc(id_AW01-1+N)%Conc(I,J,L)
ENDDO
! Get NH4 mass from the bulk mass and scale to bin with sulfate
IF ( SRTNH4 > 0 ) THEN
CALL NH4BULKTOBIN( MK(:,SRTSO4), Spc(id_NH4)%Conc(I,J,L), TRANSFER )
MK(1:IBINS,SRTNH4) = TRANSFER(1:IBINS)
Gc(SRTNH4) = Spc(id_NH3)%Conc(I,J,L)
ENDIF
! Give it the pseudo-steady state value instead later (win,9/30/08)
!GC(SRTSO4) = Spc(id_H2SO4)%Conc(I,J,L)
H2SO4rate_o = H2SO4_RATE(I,J,L) ! [kg s-1]
! Pengfei Liu add 2018/04/18, debug
IF ( H2SO4rate_o .lt. 0.e0 ) THEN
Print*, 'Debug TOMAS: H2SO4RATE = ', H2SO4rate_o, 'I = ', I, &
'J = ', J, 'L = ', L
H2SO4rate_o = 0.e+0_fp
ENDIF
!
IF ( I == 10 .and. J == 10 .and. L == 10 ) THEN
Print*, 'Debug TOMAS: H2SO4RATE =', H2SO4rate_o
ENDIF
! nitrogen and sulfur mass checks
! get the total mass of N
tot_n_1 = Gc(srtnh4)*14.e+0_fp/17.e+0_fp
do k=1,ibins
tot_n_1 = tot_n_1 + Mk(k,srtnh4)*14.e+0_fp/18.e+0_fp
enddo
! get the total mass of S
tot_s_1 = H2SO4rate_o*adt*32.e+0_fp/98.e+0_fp
do k=1,ibins
tot_s_1 = tot_s_1 + Mk(k,srtso4)*32.e+0_fp/96.e+0_fp
enddo
!if (printdebug.and.i==iob .and. j==job .and. l==lob ) then
! CALL DEBUGPRINT( Nk, Mk, I, J, L, 'Begin aerophys' )
! print *,'H2SO4RATE ',H2SO4rate_o
!endif
!*********************
! Aerosol dynamics
!*********************
!Do water eqm at appropriate times
CALL EZWATEREQM( MK, RHTOMAS )
!Fix any inconsistencies in M/N distribution (because of advection)
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
!if(printdebug .and. i==iob.and.j==job.and.l==lob) ERRORSWITCH =.TRUE.
!print *, 'mnfix in tomas_mod:533'
CALL MNFIX( NK, MK, ERRORSWITCH )
IF ( ERRORSWITCH ) THEN
PRINT *,'Aerophys: MNFIX found error at',I,J,L
CALL ERROR_STOP('AEROPHYS-MNFIX (1)','Enter microphys')
ENDIF
MPNUM = 11
IF ( State_Diag%Archive_TomasMNFIXezwat1mass .or. &
State_Diag%Archive_TomasMNFIXezwat1number ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
!IF ( printdebug.and.i==iob .and. j==job .and. l==lob ) THEN
! CALL DEBUGPRINT( Nk, Mk, I, J, L, 'After mnfix before cond/nucl' )
!ENDIF
! Before doing any cond/nucl/coag, check if there's any aerosol in
! the current box
TOT_NK = 0.e+0_fp
TOT_MK = 0.e+0_fp
do k = 1, ibins
TOT_NK = TOT_NK + Nk(K)
do jc=1, icomp-idiag
TOT_MK = TOT_MK + Mk(k,jc)
enddo
enddo
if(TOT_NK .lt. 1.e-10_fp) then
if( .NOT. SPINUP(5.0)) then
print *,'No aerosol in box ',I,J,L,'-->SKIP'
endif
CYCLE
endif
!-------------------------------------
! Condensation and nucleation (coupled)
!-------------------------------------
IF ( COND .AND. NUCL .AND. H2SO4rate_o > 0.e0_fp) THEN
!if(printdebug .and. i==iob.and.j==job.and.l==lob) ERRORSWITCH =.TRUE.
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
!sfdebug if(printdebug) then
!sfdebug !print*,'Before COND_NUC Gc(srtso4)=',Gc(srtso4)
!sfdebug do N = 1,IBINS
!sfdebug IF( IT_IS_NAN( NK(N) ) ) THEN
!sfdebug print*, "found NAN in nk", n, nk
!sfdebug endif
!sfdebug DO JC=1, ICOMP-1
!sfdebug IF( IT_IS_NAN( Gc(JC) ) ) THEN
!sfdebug print*, "found NAN in gc", jc, gc
!sfdebug endif
!sfdebug ENDDO
!sfdebug enddo
!sfdebug endif
CALL COND_NUC(Nk,Mk,Gc,Nkout,Mkout,Gcout,fn,fn1, &
H2SO4rate_o,adt,num_iter,Nknuc,Mknuc,Nkcond,Mkcond, &
ionrate, surf_area, BOXVOL, BOXMASS, TEMPTMS, PRES, &
RHTOMAS, ERRORSWITCH, l)
!sfdebug if(printdebug) then
!sfdebug !print*,'Before COND_NUC Gc(srtso4)=',Gc(srtso4)
!sfdebug do N = 1,IBINS
!sfdebug IF( IT_IS_NAN( NKout(N) ) ) THEN
!sfdebug print*, "found NAN in nkout", n, nkout
!sfdebug endif
!sfdebug DO JC=1, ICOMP-1
!sfdebug IF( IT_IS_NAN( Gcout(JC) ) ) THEN
!sfdebug print*, "found NAN in gcout", jc, gcout
!sfdebug endif
!sfdebug ENDDO
!sfdebug enddo
!sfdebug endif
IF ( ERRORSWITCH ) THEN
PRINT *,'Aerophys: found error at',I,J,L
CALL ERROR_STOP('AEROPHYS','After cond_nuc')
ENDIF
ERR_VAR = 'Gcout'
ERR_MSG = 'After COND_NUC'
! check for NaN and Inf (win, 10/4/08)
do jc = 1, icomp-1
ERR_IND(1) = I
ERR_IND(2) = J
ERR_IND(3) = L
ERR_IND(4) = 0
call check_value( Gcout(jc), ERR_IND, ERR_VAR, ERR_MSG )
!if( IT_IS_FINITE(Gcout(jc))) then
! print *,'xxxxxxxxx Found Inf in Gcout xxxxxxxxxxxxxx'
! print *,'Location ',I,J,L, 'comp',jc
! call debugprint( Nkout, Mkout, i,j,l,'After COND_NUC')
! stop
!endif
enddo
!get nucleation diagnostic
DO N = 1, IBINS
NK(N) = NKnuc(N)
DO JC = 1, ICOMP
MK(N,JC) = MKnuc(N,JC)
ENDDO
ENDDO
MPNUM = 3
IF ( State_Diag%Archive_TomasNUCLmass .or. &
State_Diag%Archive_TomasNUCLnumber ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
MPNUM = 7
IF ( State_Diag%Archive_TomasNUCRATEnumber) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
IF ( printdebug.and.i==iob .and. j==job .and. l==lob ) THEN
CALL DEBUGPRINT( Nk, Mk, I, J, L,'After nucleation' )
ENDIF
!get condensation diagnostic
DO N = 1, IBINS
NK(N) = NKcond(N)
DO JC = 1, ICOMP
MK(N,JC) = MKcond(N,JC)
ENDDO
ENDDO
Gc(srtnh4)=Gcout(srtnh4)
Gc(srtso4)=Gcout(srtso4)
MPNUM = 1
IF ( State_Diag%Archive_TomasNUCLmass .or. &
State_Diag%Archive_TomasNUCLnumber ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
IF ( printdebug.and.i==iob .and. j==job .and. l==lob ) THEN
CALL DEBUGPRINT( Nk, Mk, I, J, L,'After condensation' )
ENDIF
nucrate(j,l)=nucrate(j,l)+fn
nucrate1(j,l)=nucrate1(j,l)+fn1
! replaces old ND61 diagnostic!
IF ( State_Diag%Archive_TomasNUCRATEFN ) THEN
State_Diag%TomasNUCRATEFN(I,J,L) = fn
ENDIF
#ifdef BPCH_DIAG
IF ( ND61 > 0 ) THEN
IF ( L <= LD61 ) AD61(I,J,L,2) = AD61(I,J,L,2) + fn
! Tracks nucleation rates instantaneously for planeflight
AD61_INST(I,J,L,2) = fn
ENDIF
#endif
DO N = 1, IBINS
NK(N) = NKout(N)
DO JC = 1, ICOMP
MK(N,JC) = MKout(N,JC)
ENDDO
ENDDO
ENDIF ! end of cond and nuc !
! nitrogen and sulfur mass checks
! get the total mass of N
tot_n_1a = Gc(srtnh4)*14.e+0_fp/17.e+0_fp
do k=1,ibins
tot_n_1a = tot_n_1a + Mk(k,srtnh4)*14.e+0_fp/18.e+0_fp
enddo
! get the total mass of S
tot_s_1a = 0.e+0_fp
do k=1,ibins
tot_s_1a = tot_s_1a + Mk(k,srtso4)*32.e+0_fp/96.e+0_fp
enddo
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
!print *, 'mnfix in tomas_mod:677'
CALL MNFIX( Nk, Mk, ERRORSWITCH )
IF ( ERRORSWITCH ) THEN
PRINT *,'Aerophys: MNFIX found error at',I,J,L
IF( .not. SPINUP(14.0) ) THEN
CALL ERROR_STOP('AEROPHYS-MNFIX (2)','After cond/nucl')
ELSE
PRINT *,'Let error go during spin up'
ENDIF
ENDIF
MPNUM = 14
IF ( State_Diag%Archive_TomasMNFIXh2so4mass .or. &
State_Diag%Archive_TomasMNFIXh2so4number ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
!-----------------------------
! Coagulation
!-----------------------------
if(printdebug .and. i==iob.and.j==job.and.l==lob) ERRORSWITCH =.TRUE.
!if (i==iob .and. j==job .and. l==lob ) &
! CALL DEBUGPRINT( Nk, Mk, I, J, L,'Before coagulation' )
IF( COAG ) THEN
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
CALL MULTICOAG( ADT, Nk, Mk, BOXVOL, PRES, TEMPTMS, errorswitch )
if ( errorswitch ) &
CALL DEBUGPRINT( Nk, Mk, I, J, L,'After coagulation' )
!if (i==iob .and. j==job .and. l==lob ) &
! CALL DEBUGPRINT( Nk, Mk, I, J, L,'After coagulation' )
MPNUM = 2
IF ( State_Diag%Archive_TomasCOAGmass .or. &
State_Diag%Archive_TomasCOAGnumber ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
!Fix any inconsistency after coagulation (win, 4/18/06)
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
if(printdebug .and. i==iob.and.j==job.and.l==lob) &
ERRORSWITCH=.true. !4/18/06 win
!print *, 'mnfix in tomas_mod:719'
CALL MNFIX( NK, MK, ERRORSWITCH )
IF ( ERRORSWITCH ) THEN
PRINT *,'MNFIX found error at',I,J,L
IF( .not. SPINUP(14.0) ) THEN
CALL ERROR_STOP('AEROPHYS-MNFIX (3)', 'After COAGULATION' )
ELSE
PRINT *,'Let error go during spin up'
ENDIF
ENDIF
MPNUM = 15
IF ( State_Diag%Archive_TomasMNFIXcoagmass .or. &
State_Diag%Archive_TomasMNFIXcoagnumber ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
ENDIF ! Coagulation
! Do water eqm at appropriate times
CALL EZNH3EQM( Gc, Mk )
CALL EZWATEREQM ( MK, RHTOMAS )
!****************************
! End of aerosol dynamics
!****************************
!Fix any inconsistencies in M/N distribution (because of advection)
CALL STORENM(Nk, Nkd, Mk, Mkd, Gc, Gcd)
! Make sure anything that leaves AEROPHYS is free of any error
! This MNFIX call could be temporary (?) or just leave it here and
! monitor if the error fixed is significantly large meaning some
! serious problem needs to be investigated
if(printdebug .and. i==iob.and.j==job.and.l==lob) ERRORSWITCH =.true.
!print *, 'mnfix in tomas_mod:758'
CALL MNFIX(NK,MK,ERRORSWITCH)
IF ( ERRORSWITCH ) THEN
PRINT *,'End of Aerophys: MNFIX found error at',I,J,L
IF( .not. SPINUP(14.0) ) THEN
CALL ERROR_STOP('AEROPHYS-MNFIX (4)', 'End of microphysics')
ELSE
PRINT *,'Let error go during spin up'
ENDIF
ENDIF
! Accumulate changes by mnfix to diagnostic (win, 9/8/05)
MPNUM = 12
IF ( State_Diag%Archive_TomasMNFIXezwat2mass .or. &
State_Diag%Archive_TomasMNFIXezwat2number ) THEN
CALL AERODIAG( MPNUM, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
State_Grid, State_Diag )
ENDIF
! Swap Nk, Mk, and Gc arrays back to Spc
DO N = 1, IBINS
TRACNUM = id_NK01 - 1 + N
Spc(TRACNUM)%Conc(I,J,L) = NK(N)
DO JC = 1, ICOMP-IDIAG
TRACNUM = id_NK01 - 1 + N + IBINS * JC
Spc(TRACNUM)%Conc(I,J,L) = MK(N,JC)
ENDDO
Spc(id_AW01-1+N)%Conc(I,J,L) = MK(N,SRTH2O)
ENDDO
Spc(id_H2SO4)%Conc(I,J,L) = GC(SRTSO4)
! print to file to check mass conserv
!write(*,77) I,J,L, Spc(id_NH3)%Conc(I,J,L), Spc(id_NH3)%Conc(I,J,L)-GC(SRTNH4)
! Calculate NH3 gas lost to aerosol phase as NH4
NH3_to_NH4 = Spc(id_NH3)%Conc(I,J,L)-GC(SRTNH4)
! Update the bulk NH4 aerosol species
if ( NH3_to_NH4 > 0e+0_fp ) &
Spc(id_NH4)%Conc(I,J,L) = Spc(id_NH4)%Conc(I,J,L) + &
NH3_to_NH4/17.e+0_fp*18.e+0_fp
! Update NH3 gas species (win, 10/6/08)
! plus tiny amount CEPS in case zero causes some problem
Spc(id_NH3)%Conc(I,J,L) = GC(SRTNH4) + CEPS !MUST CHECK THIS!! (win,9/26/08)
!vbn write(889,89)I,J,L,Spc(id_H2SO4)%Conc(I,J,L)
89 format(3I3,'Spc(id_H2SO4)%Conc(I,J,L) kg', E13.5)
ENDDO !L loop
ENDDO !J loop
ENDDO !I loop
!$OMP END PARALLEL DO
!WRITE(777,*) '---------------------------'
77 FORMAT(3I4, ' Spc(id_NH3)%Conc(I,J,L),'E13.5,' Used', E13.5 )
IF ( COND .and. Input_Opt%Verbose ) THEN
PRINT *,'### AEROPHYS: SO4 CONDENSATION'
ENDIF
IF ( COAG .and. Input_Opt%Verbose ) THEN
PRINT *,'### AEROPHYS: COAGULATION'
ENDIF
IF ( NUCL .and. Input_Opt%Verbose ) THEN
PRINT *,'### AEROPHYS: NUCLEATION'
ENDIF
! Free pointer memory
Spc => NULL()
END SUBROUTINE AEROPHYS
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: cond_nuc
!
! !DESCRIPTION: This subroutine calculates the change in the aerosol size
! distribution due to so4 condensation and binary/ternary nucleation during
! the overal microphysics timestep.
! WRITTEN BY Jeff Pierce, May 2007 for GISS GCM-II'