-
Notifications
You must be signed in to change notification settings - Fork 0
/
MARS_mod.f90
executable file
·3913 lines (3294 loc) · 154 KB
/
MARS_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
! <MARS_mod.f90 - A component of the EMEP MSC-W Unified Eulerian
! Chemical transport Model>
!*****************************************************************************!
!*
!* Copyright (C) 2007-2011 met.no
!*
!* Contact information:
!* Norwegian Meteorological Institute
!* Box 43 Blindern
!* 0313 OSLO
!* NORWAY
!* email: [email protected]
!* http://www.emep.int
!*
!* This program is free software: you can redistribute it and/or modify
!* it under the terms of the GNU General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
!* This program is distributed in the hope that it will be useful,
!* but WITHOUT ANY WARRANTY; without even the implied warranty of
!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!* GNU General Public License for more details.
!*
!* You should have received a copy of the GNU General Public License
!* along with this program. If not, see <http://www.gnu.org/licenses/>.
!*****************************************************************************!
module MARS_mod
! -----------------------------------------------------------------------
! Calculates gas-aerosol equilibrium for SO4, HNO3-NO3 and NH3-NH4 system
! Made available by Frank Binkowski (originally from EPA's RPM model)
! Presently not in use, EQSAM is used for inorganic equilibrium stuff
!------------------------------------------------------------------------
use CheckStop_mod, only: CheckStop
use Io_mod, only: ios, datewrite
use MARS_Aero_water_mod, only: Awater!,awater_2900
use Config_module, only: NPROC
use Debug_module, only: DEBUG ! -> DEBUG%EQUIB
use Par_mod, only: me ! processor number
implicit none
private
!to test without the smoothing between different TNH4/TSO4 regimes
!before rv4.2beta (Jan 2013), the code was (or should be) equivalent to MARS_RATIO_SMOOTH=.false.
logical, parameter :: MARS_RATIO_SMOOTH=.true.
real, parameter :: FLOOR = 1.0E-30 ! minimum concentration
! -30 from RPM
!/- subroutines:
public :: rpmares, &
rpmares_2900, &!svn version 2908, for testing if there are significant differences. will be deleted
rpmares_new, &!GEOSCHEM adapted version
DO_RPMARES_new, &!driver for rpmares_new
awater_new, &!almost identical to awater. Only difference Y2 = ( 1.0d0 - MFSSO4 ) / MFSSO4?
cubic, &
actcof, &
IS_SAFE_DIV
integer, private, save :: MAXNNN1 = 0
integer, private, save :: MAXNNN2 = 0
!ds real MWNO3 ! molecular weight for NO3
real, private, parameter :: MWNO3 = 62.0049
!ds real MWHNO3 ! molecular weight for HNO3
real, private, parameter :: MWHNO3 = 63.01287
!ds real MWSO4 ! molecular weight for SO4
real, private, parameter :: MWSO4 = 96.0576
!ds real MWHSO4 ! molecular weight for HSO4
real, private, parameter :: MWHSO4 = MWSO4 + 1.0080
!ds real MH2SO4 ! molecular weight for H2SO4
real, private, parameter :: MH2SO4 = 98.07354
!ds real MWNH3 ! molecular weight for NH3
real, private, parameter :: MWNH3 = 17.03061
!ds real MWNH4 ! molecular weight for NH4
real, private, parameter :: MWNH4 = 18.03858
contains
real function poly4 (a,x)
! Calculates the polynomial based on 4 coefficients a(1:4):
!-- arguments
real, dimension(4), intent(in) :: a
real, intent(in) :: x
poly4 = a(1) + x * ( a(2) + x * ( a(3) + x * ( a(4) ) ) )
end function poly4
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
real function poly6(a,x)
! Calculates the polynomial based on 6 coefficients a(1:6):
!-- arguments
real, dimension(6), intent(in) :: a
real, intent(in) :: x
poly6 = a(1) + x * ( a(2) + x * ( a(3) + x * ( a(4) + &
x * ( a(5) + x * (a(6) ) ) ) ) )
end function poly6
!from GEOS-Chem
!Could/should be moved to another module
FUNCTION IS_SAFE_DIV( N, D, R4 ) RESULT( F )
!
! !INPUT PARAMETERS:
!
REAL*8, INTENT(IN) :: N ! Numerator
REAL*8, INTENT(IN) :: D ! Denominator
LOGICAL, INTENT(IN), OPTIONAL :: R4 ! Logical flag to use the limits
! of REAL*4 to define underflow
! or overflow. Extra defensive.
!
! !OUTPUT PARAMETERS:
!
LOGICAL :: F ! =F if division isn't allowed
! =T otherwise
!
! !REMARKS:
! UnderFlow, OverFlow and NaN are tested for. If you need to
! differentiate between the three, use the SAFE_DIV (phs, 4/14/09)
!
! !REVISION HISTORY:
! 11 Jun 2008 - P. Le Sager - Initial version
! 20 Nov 2009 - R. Yantosca - Added ProTeX header
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
INTEGER MaxExp, MinExp
REAL*4 RR
!==================================================================
! IS_SAFE_DIV begins here!
!==================================================================
MaxExp = MAXEXPONENT( N )
MinExp = MINEXPONENT( N )
IF ( PRESENT( R4 ) ) THEN
IF ( R4 ) THEN
MaxExp = MAXEXPONENT( RR )
MinExp = MINEXPONENT( RR )
ENDIF
ENDIF
IF ( EXPONENT(N) - EXPONENT(D) >= MaxExp .or. D==0 .or.&
EXPONENT(N) - EXPONENT(D) <= MinExp ) THEN
F = .FALSE.
ELSE
F = .TRUE.
ENDIF
! Return to calling program
END FUNCTION IS_SAFE_DIV
!>-------------------------------------------------------------------------------<
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
subroutine rpmares ( SO4, HNO3, NO3, NH3, NH4, RH, TEMP, &
ASO4, ANO3, AH2O, ANH4, GNH3, GNO3, &
ERRMARK,debug_flag)
!-----------------------------------------------------------------------
!C
!C Description:
!C
!C ARES calculates the chemical composition of a sulfate/nitrate/
!C ammonium/water aerosol based on equilibrium thermodynamics.
!C
!C This code considers two regimes depending upon the molar ratio
!C of ammonium to sulfate.
!C
!C For values of this ratio less than 2,the code solves a cubic for
!C hydrogen ion molality, HPLUS, and if enough ammonium and liquid
!C water are present calculates the dissolved nitric acid. For molal
!C ionic strengths greater than 50, nitrate is assumed not to be present.
!C
!C For values of the molar ratio of 2 or greater, all sulfate is assumed
!C to be ammonium sulfate and a calculation is made for the presence of
!C ammonium nitrate.
!C
!C The Pitzer multicomponent approach is used in subroutine ACTCOF to
!C obtain the activity coefficients. Abandoned -7/30/97 FSB
!c The Bromley method of calculating the activity coefficients is s used
!c in this version
!c The calculation of liquid water
!C is done in subroutine water. Details for both calculations are given
!C in the respective subroutines.
!C
!C Based upon MARS due to
!C P. Saxena, A.B. Hudischewskyj, C. Seigneur, and J.H. Seinfeld,
!C Atmos. Environ., vol. 20, Number 7, Pages 1471-1483, 1986.
!C
!C and SCAPE due to
!C Kim, Seinfeld, and Saxeena, Aerosol Ceience and Technology,
!C Vol 19, number 2, pages 157-181 and pages 182-198, 1993.
!C
!C NOTE: All concentrations supplied to this subroutine are TOTAL
!C over gas and aerosol phases
!C
!C Parameters:
!C
!C SO4 : Total sulfate in MICROGRAMS/M**3 as sulfate (IN)
!C HNO3 : Nitric Acid in MICROGRAMS/M**3 as nitric acid (IN)
!C NO3 : Total nitrate in MICROGRAMS/M**3 as nitric acid (IN)
!C NH3 : Total ammonia in MICROGRAMS/M**3 as ammonia (IN)
!C NH4 : Ammonium in MICROGRAMS/M**3 as ammonium (IN)
!C RH : Fractional relative humidity (IN)
!C TEMP : Temperature in Kelvin (IN)
!C GNO3 : Gas phase nitric acid in MICROGRAMS/M**3 (OUT)
!C GNH3 : Gas phase ammonia in MICROGRAMS/M**3 (OUT)
!C ASO4 : Aerosol phase sulfate in MICROGRAMS/M**3 (OUT)
!C ANO3 : Aerosol phase nitrate in MICROGRAMS/M**3 (OUT)
!C ANH4 : Aerosol phase ammonium in MICROGRAMS/M**3 (OUT)
!C AH2O : Aerosol phase water in MICROGRAMS/M**3 (OUT)
!C NITR : Number of iterations for obtaining activity coefficients (OUT)
!C NR : Number of real roots to the cubic in the low ammonia case (OUT)
!C
!C Revision History:
!C Who When Detailed description of changes
!C --------- -------- -------------------------------------------
!C S.Roselle 11/10/87 Received the first version of the MARS code
!C S.Roselle 12/30/87 Restructured code
!C S.Roselle 2/12/88 Made correction to compute liquid-phase
!C concentration of H2O2.
!C S.Roselle 5/26/88 Made correction as advised by SAI, for
!C computing H+ concentration.
!C S.Roselle 3/1/89 Modified to operate with EM2
!C S.Roselle 5/19/89 Changed the maximum ionic strength from
!C 100 to 20, for numerical stability.
!C F.Binkowski 3/3/91 Incorporate new method for ammonia rich case
!C using equations for nitrate budget.
!C F.Binkowski 6/18/91 New ammonia poor case which
!C omits letovicite.
!C F.Binkowski 7/25/91 Rearranged entire code, restructured
!C ammonia poor case.
!C F.Binkowski 9/9/91 Reconciled all cases of ASO4 to be output
!C as SO4--
!C F.Binkowski 12/6/91 Changed the ammonia defficient case so that
!C there is only neutralized sulfate (ammonium
!C sulfate) and sulfuric acid.
!C F.Binkowski 3/5/92 Set RH bound on AWAS to 37 % to be in agreement
!C with the Cohen et al. (1987) maximum molality
!C of 36.2 in Table III.( J. Phys Chem (91) page
!C 4569, and Table IV p 4587.)
!C F.Binkowski 3/9/92 Redid logic for ammonia defficient case to remove
!C possibility for denomenator becoming zero;
!C this involved solving for HPLUS first.
!C Note that for a relative humidity
!C less than 50%, the model assumes that there is no
!C aerosol nitrate.
!C F.Binkowski 4/17/95 Code renamed ARES (AeRosol Equilibrium System)
!C Redid logic as follows
!C 1. Water algorithm now follows Spann & Richardson
!C 2. Pitzer Multicomponent method used
!C 3. Multicomponent practical osmotic coefficient
!C use to close iterations.
!C 4. The model now assumes that for a water
!C mass fraction WFRAC less than 50% there is
!C no aerosol nitrate.
!C F.Binkowski 7/20/95 Changed how nitrate is calculated in ammonia poor
!C case, and changed the WFRAC criterion to 40%.
!C For ammonium to sulfate ratio less than 1.0
!C all ammonium is aerosol and no nitrate aerosol
!C exists.
!C F.Binkowski 7/21/95 Changed ammonia-ammonium in ammonia poor case to
!C allow gas-phase ammonia to exist.
!C F.Binkowski 7/26/95 Changed equilibrium constants to values from
!C Kim et al. (1993)
!C F.Binkowski 6/27/96 Changed to new water format
!c F.Binkowski 7/30/97 Changed to Bromley method for multicomponent
!c activity coefficients. The binary activity coefficients
!c are the same as the previous version
!c F.Binkowski 8/1/97 Chenged minimum sulfate from 0.0 to 1.0e-6 i.e.
!c 1 picogram per cubic meter
!C I.Ackermann 2/23/98 modification for solving precision problems
!C on workstations
!C I.Ackermann 2/25/99 changed logic as follows:
!c If iterations fail, initial values of nitrate
!c are retained.
!c Total mass budgets are changed to account for gas
!c phase returned. (incorporated from FZB's models3
!c framework)
!C eliminated ratio=5 !! for low to zero sulfate
!C I.Ackermann 3/17/99 modified ratio = 5 treatment see RB3,p.125
!C
!C-----------------------------------------------------------------------
!...........ARGUMENTS and their descriptions
real, intent(in) :: SO4 & ! Total sulfate in micrograms / m**3
,HNO3 & ! Total nitric acid in micrograms / m**3
,NO3 & ! Total nitrate in micrograms / m**3
,NH3 & ! Total ammonia in micrograms / m**3
,NH4 & ! Total ammonium in micrograms / m**3
,RH & ! Fractional relative humidity
,TEMP ! Temperature in Kelvin
real, intent(out):: ASO4 & ! Aerosol sulfate in micrograms / m**3
,ANO3 & ! Aerosol nitrate in micrograms / m**3
,AH2O & ! Aerosol liquid water content water in micrograms / m**3
,ANH4 & ! Aerosol ammonium in micrograms / m**3
,GNO3 & ! Gas-phase nitric acid in micrograms / m**3
,GNH3 ! Gas-phase ammonia in micrograms / m**3
logical, intent(in) :: debug_flag
!C...........INCLUDES and their descriptions
!! INCLUDE SUBST_CONST ! constants
!...........PARAMETERS and their descriptions:
real MWNACL ! molecular weight for NaCl
parameter ( MWNACL = 58.44277 )
!ds moved a bunch upstairs.
real MWORG ! molecular weight for Organic Species
parameter ( MWORG = 16.0 )
real MWCL ! molecular weight for Chloride
parameter ( MWCL = 35.453 )
real MWAIR ! molecular weight for AIR
parameter ( MWAIR = 28.964 )
real MWLCT ! molecular weight for Letovicite
parameter ( MWLCT = 3.0 * MWNH4 + 2.0 * MWSO4 + 1.0080 )
real MWAS ! molecular weight for Ammonium Sulfate
parameter ( MWAS = 2.0 * MWNH4 + MWSO4 )
real MWABS ! molecular weight for Ammonium Bisulfate
parameter ( MWABS = MWNH4 + MWSO4 + 1.0080 )
!...........SCRATCH LOCAL VARIABLES and their descriptions:
REAL fRH ! Index set to percent relative humidity
INTEGER NITR ! Number of iterations for activity coefficients
INTEGER NNN ! Loop index for iterations
INTEGER NR ! Number of roots to cubic equation for HPLUS
INTEGER ERRMARK
real A0 ! Coefficients and roots of
real A1 ! Coefficients and roots of
real A2 ! Coefficients and roots of
real AA ! Coefficients and discriminant for quadratic equation for ammonium nitrate
real BAL ! internal variables ( high ammonia case)
real BB ! Coefficients and discriminant for quadratic equation for ammonium nitrate
real BHAT ! Variables used for ammonia solubility
real CC ! Coefficients and discriminant for quadratic equation for ammonium nitrate
real CONVT ! Factor for conversion of units
real DD ! Coefficients and discriminant for quadratic equation for ammonium nitrate
real DISC ! Coefficients and discriminant for quadratic equation for ammonium nitrate
real EROR ! Relative error used for convergence test
real FNH3 ! "Free ammonia concentration", that which exceeds TWOSO4
real GAMAAB ! Activity Coefficient for (NH4+, HSO4-)GAMS( 2,3 )
real GAMAAN ! Activity coefficient for (NH4+, NO3-) GAMS( 2,2 )
real GAMAHAT ! Variables used for ammonia solubility
real GAMANA ! Activity coefficient for (H+ ,NO3-) GAMS( 1,2 )
real GAMAS1 ! Activity coefficient for (2H+, SO4--) GAMS( 1,1 )
real GAMAS2 ! Activity coefficient for (H+, HSO4-) GAMS( 1,3 )
real GAMOLD ! used for convergence of iteration
real GASQD ! internal variables ( high ammonia case)
real HPLUS ! Hydrogen ion (low ammonia case) (moles / kg water)
real K1A ! Equilibrium constant for ammoniua to ammonium
real K2SA ! Equilibrium constant for sulfate-bisulfate (aqueous)
real K3 ! Dissociation constant for ammonium nitrate
real KAN ! Equilibrium constant for ammonium nitrate (aqueous)
real KHAT ! Variables used for ammonia solubility
real KNA ! Equilibrium constant for nitric acid (aqueous)
real KPH ! Henry's Law Constant for ammonia
real KW ! Equilibrium constant for water dissociation
real KW2 ! Internal variable using KAN
real MAN ! Nitrate (high ammonia case) (moles / kg water)
real MAS ! Sulfate (high ammonia case) (moles / kg water)
real MHSO4 ! Bisulfate (low ammonia case) (moles / kg water)
real MNA ! Nitrate (low ammonia case) (moles / kg water)
real MNH4 ! Ammonium (moles / kg water)
real MOLNU ! Total number of moles of all ions
real MSO4 ! Sulfate (low ammonia case) (moles / kg water)
real PHIBAR ! Practical osmotic coefficient
real PHIOLD ! Previous value of practical osmotic coefficient used for convergence of iteration
real RATIO ! Molar ratio of ammonium to sulfate
real RK2SA ! Internal variable using K2SA
real RKNA ! Internal variables using KNA
real RKNWET ! Internal variables using KNA
real RR1
real RR2
real STION ! Ionic strength
real T1 ! Internal variables for temperature corrections
real T2 ! Internal variables for temperature corrections
real T21 ! Internal variables of convenience (low ammonia case)
real T221 ! Internal variables of convenience (low ammonia case)
real T3 ! Internal variables for temperature corrections
real T4 ! Internal variables for temperature corrections
real T6 ! Internal variables for temperature corrections
real TNH4 ! Total ammonia and ammonium in micromoles / meter ** 3
real TNO3 ! Total nitrate in micromoles / meter ** 3
real TOLER1 ! Tolerances for convergence test
real TOLER2 ! Tolerances for convergence test
real TSO4 ! Total sulfate in micromoles / meter ** 3
real TWOSO4 ! 2.0 * TSO4 (high ammonia case) (moles / kg water)
real WFRAC ! Water mass fraction
real WH2O ! Aerosol liquid water content (internally)
! micrograms / meter **3 on output
! internally it is 10 ** (-6) kg (water) / meter ** 3
! the conversion factor (1000 g = 1 kg) is applied
! for AH2O output
real WSQD ! internal variables ( high ammonia case)
real XNO3 ! Nitrate aerosol concentration in micromoles / meter ** 3
real XXQ ! Variable used in quadratic solution
real YNH4 ! Ammonium aerosol concentration in micromoles / meter** 3
real ZH2O ! Water variable saved in case ionic strength too high.
real ZSO4 ! Total sulfate molality - mso4 + mhso4 (low ammonia case) (moles / kg water)
real CAT( 2 ) ! Array for cations (1, H+); (2, NH4+) (moles / kg water)
real AN ( 3 ) ! Array for anions (1, SO4--); (2, NO3-); (3, HSO4-) (moles / kg water)
real CRUTES( 3 ) ! Coefficients and roots of
real GAMS( 2, 3 ) ! Array of activity coefficients
real MINSO4 ! Minimum value of sulfate laerosol concentration
parameter( MINSO4 = 1.0E-6 / MWSO4 )
real MINNO3
parameter( MINNO3 = 1.0E-6 / MWNO3 ) !2/25/99 IJA
!st real FLOOR
!st parameter( FLOOR = 1.0E-30) ! minimum concentration
!2/25/99 IJA
! FSB New variables Total ammonia and nitrate mass concentrations
real TMASSNH3 ! Total ammonia (gas and particle)
! as ammonia gas [ug /m**3]
real TMASSHNO3 ! Total nitrate (vapor and particle) as
! nitric acid vapor [ug /m**3]
!for RATIO between RATIO_Low and RATIO_High, consider the species in two boxes one with
! RATIO_Low and the other with RATIO_High. Then take a linear combination of the results in each box
logical ::DO_RATIO_High_2,DO_RATIO_Low_2
real, parameter ::RATIO_Low = 1.95,RATIO_High= 2.05 !somewhat arbitrarily
real ::ASO4_High, ANO3_High, AH2O_High, ANH4_High, GNH3_High, GNO3_High
real ::ASO4_High1, ANO3_High1, AH2O_High1, ANH4_High1, GNH3_High1, GNO3_High1,diff1
real ::ASO4_Low, ANO3_Low, AH2O_Low, ANH4_Low, GNH3_Low, GNO3_Low
real ::ASO4_Low1, ANO3_Low1, AH2O_Low1, ANH4_Low1, GNH3_Low1, GNO3_Low1,diff2
real::TSO4_HighA,TSO4_LowA,High_Factor,X
!-----------------------------------------------------------------------
integer, save :: nMarsErrors = 0 ! tracks solution failures DSMARS
!-----------------------------------------------------------------------
! begin body of subroutine RPMARES
!ASO4=FLOOR;ANO3=FLOOR;AH2O=FLOOR;ANH4=FLOOR;GNO3=FLOOR;GNH3=FLOOR
!Initialise the output variables
ASO4=0.0;ANO3=0.0;AH2O=0.0;ANH4=0.0;GNO3=0.0;GNH3=0.0
ASO4 = SO4 !ds from RPM
!...convert into micromoles/m**3
!..iamodels3 merge NH3/NH4 , HNO3,NO3 here
TSO4 = MAX( 0.0, SO4 / MWSO4 )
TSO4 = MAX( FLOOR, TSO4 ) !GEOS added
TNO3 = MAX( 0.0, (NO3 / MWNO3 + HNO3 / MWHNO3) )
TNH4 = MAX( 0.0, (NH3 / MWNH3 + NH4 / MWNH4) )
!2/25/99 IJA
! TMASSNH3 = MAX(0.0, NH3 + (MWNH3 / MWNH4) * NH4 )
! TMASSHNO3 = MAX(0.0, NO3 + (MWHNO3 / MWNO3) * NO3 )
TMASSNH3 = MAX(0.0, NH3 + NH4 )
TMASSHNO3 = MAX(0.0, HNO3 + NO3 )
!...now set humidity index fRH as a percent
!st IRH = NINT( 100.0 * RH )
fRH = RH
!GEOS added
! For extremely low relative humidity ( less than 1% ) set the
! water content to a minimum and skip the calculation.
IF ( fRH .LT. 0.01 ) THEN
AH2O = FLOOR
!not from GEOS:
ASO4 = SO4
ANO3 = NO3
WH2O = 0.0
GNH3 = NH3
ANH4 = NH4
GNO3 = HNO3
RETURN
ENDIF
!...Check for valid fRH
fRH = MAX( 0.01, fRH )
fRH = MIN( 0.99, fRH )
!...Specify the equilibrium constants at correct
!... temperature. Also change units from ATM to MICROMOLE/M**3 (for KAN,
!... KPH, and K3 )
!... Values from Kim et al. (1993) except as noted.
CONVT = 1.0 / ( 0.082 * TEMP )
T6 = 0.082E-9 * TEMP
T1 = 298.0 / TEMP
T2 = ALOG( T1 )
T3 = T1 - 1.0
T4 = 1.0 + T2 - T1
KNA = 2.511E+06 * EXP( 29.17 * T3 + 16.83 * T4 ) * T6
K1A = 1.805E-05 * EXP( -1.50 * T3 + 26.92 * T4 )
K2SA = 1.015E-02 * EXP( 8.85 * T3 + 25.14 * T4 )
KW = 1.010E-14 * EXP( -22.52 * T3 + 26.92 * T4 )
KPH = 57.639 * EXP( 13.79 * T3 - 5.39 * T4 ) * T6
!!! K3 = 5.746E-17 * EXP( -74.38 * T3 + 6.12 * T4 ) * T6 * T6
KHAT = KPH * K1A / KW
KAN = KNA * KHAT
!...Compute temperature dependent equilibrium constant for NH4NO3
!... ( from Mozurkewich, 1993)
K3 = EXP( 118.87 - 24084.0 / TEMP - 6.025 * ALOG( TEMP ) )
!...Convert to (micromoles/m**3) **2
K3 = K3 * CONVT * CONVT
WH2O = 0.0
STION = 0.0
AH2O = 0.0
MAS = 0.0
MAN = 0.0
HPLUS = 0.0
TOLER1 = 0.00001
TOLER2 = 0.001
NITR = 0
NR = 0
RATIO = 0.0
GAMAAN = 1.0
GAMOLD = 1.0
!ds from RPM, but removed FLOOR : (slighty different logic)
if( (TSO4 <= MINSO4 ) .and. (TNO3 < MINNO3) ) then
!print *, "DSX HERE"
ASO4 = SO4 ! MAX..
ANO3 = NO3 ! MAX..
WH2O = 0.0
AH2O = 0.0
GNH3 = NH3 ! MAX(FLOOR,NH3)
GNO3 = HNO3 ! MAX(FLOOR,NO3)
ANH4 = NH4
RETURN
END IF
!ds end rpm
!...set the ratio according to the amount of sulfate and nitrate
IF ( TSO4 > MINSO4 ) THEN
RATIO = TNH4 / TSO4
!...If there is no sulfate and no nitrate, there can be no ammonium
!... under the current paradigm. Organics are ignored in this version.
ELSE
IF ( TNO3 <= MINNO3 ) THEN
! *** If there is very little sulfate and nitrate set concentrations
! to a very small value and return.
! Jun 2012, Note these values are set in the initialisation
ASO4 = SO4 ! MAX(FLOOR, ASO4)
ANO3 = NO3 ! MAX(FLOOR, ANO3 )
WH2O = 0.0
AH2O = 0.0
GNH3 = NH3 ! MAX(FLOOR,GNH3)
GNO3 = HNO3 ! MAX(FLOOR,GNO3)
ANH4 = NH4
RETURN
END IF
!...For the case of no sulfate and nonzero nitrate, set ratio to 5
!... to send the code to the high ammonia case if there is more
!... ammonia than sulfate, otherwise send to low ammonia case.
IF (TNH4 > TSO4) THEN
RATIO = 5.0 !this is a high ammonia case with low sulfur
ELSE
RATIO = 1. !this is a low ammonia case with low sulfur
END IF
END IF
DO_RATIO_High_2=.false.
DO_RATIO_Low_2=.false.
ASO4_High = 0.0
ANO3_High = 0.0
AH2O_High = 0.0
ANH4_High = 0.0
GNH3_High = 0.0
GNO3_High = 0.0
ASO4_Low = 0.0
ANO3_Low = 0.0
AH2O_Low = 0.0
ANH4_Low = 0.0
GNH3_Low = 0.0
GNO3_Low = 0.0
TSO4_HighA = TSO4
TSO4_LowA = TSO4
IF ( RATIO > 2.0 )then
DO_RATIO_High_2=.true.
High_Factor=1.0
else
DO_RATIO_Low_2=.true.
High_Factor=0.0
end if
IF ( RATIO >RATIO_Low .and. RATIO < RATIO_High .and. MARS_RATIO_SMOOTH)then
DO_RATIO_High_2=.true.
DO_RATIO_Low_2=.true.
TSO4_HighA=TSO4*Ratio/RATIO_High
TSO4_LowA=TSO4*Ratio/RATIO_Low
! High_Factor=(RATIO-RATIO_Low)/(RATIO_High-RATIO_Low)
High_Factor=(RATIO*RATIO_High-RATIO_Low*RATIO_High)/(RATIO*RATIO_High-RATIO*RATIO_Low)
end if
!....................................
!......... High Ammonia Case ........
!....................................
! IF ( RATIO > 2.0 ) THEN ! NH4/SO4 > 2
IF ( DO_RATIO_High_2 ) THEN ! NH4/SO4 > 2
GAMAAN = 0.1
!...Set up twice the sulfate for future use.
TWOSO4 = 2.0 * TSO4_HighA
XNO3 = 0.0
YNH4 = TWOSO4
!...Treat different regimes of relative humidity
!...ZSR relationship is used to set water levels. Units are
!... 10**(-6) kg water/ (cubic meter of air)
!... start with ammomium sulfate solution without nitrate
CALL awater(fRH,TSO4_HighA,YNH4,TNO3,AH2O ) !**** note TNO3
WH2O = 1.0E-3 * AH2O
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = 0.0
ANH4_High = YNH4 * MWNH4
if(debug_flag) call datewrite("MARS debug ", -1,(/ ASO4_High, ANH4_High, AH2O /) )
if ( DEBUG%EQUIB ) then
if( ASO4_High + ANH4_High + AH2O < 1.0-10 ) then
call datewrite("MARS failing? ", -1,(/ ASO4_High, ANH4_High, AH2O /) )
print *, "MARS PROB ", ASO4_High, ANH4_High, AH2O, TSO4_HighA, YNH4
call CheckStop("MARS")
end if
end if
WFRAC = (AH2O + FLOOR) / ( ASO4_High + ANH4_High + AH2O + FLOOR )
!ds WFRAC = AH2O / ( ASO4 + ANH4 + AH2O + FLOOR )
!CRUDE FIX? WFRAC = AH2O / ( ASO4 + ANH4 + AH2O )
!!!! IF ( WFRAC == 0.0 ) RETURN ! No water
IF ( WFRAC < 0.2 ) THEN
!..."dry" ammonium sulfate and ammonium nitrate
!... compute free ammonia
FNH3 = TNH4 - TWOSO4
CC = TNO3 * FNH3 - K3
!...check for not enough to support aerosol
!dsjIF ( CC <= 0.0 ) THEN
IF ( CC < FLOOR ) THEN
XNO3 = 0.0
ELSE
AA = 1.0
BB = -( TNO3 + FNH3 )
DISC = BB * BB - 4.0 * CC
!...Check for complex roots of the quadratic
!... set nitrate to zero and RETURN if complex roots are found
!2/25/99 IJA
!DS IF ( DISC < 0.0 ) THEN
!dsDSdIF ( DISC < FLOOR ) THEN
IF ( DISC < 0.0 ) THEN
if( DEBUG%EQUIB .and. debug_flag ) print *, "MARS DISC NEG ", XNO3, WH2O, DISC
XNO3 = 0.0
AH2O_High = 1000.0 * WH2O
YNH4 = TWOSO4
GNO3_High = HNO3
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = NO3
ANH4_High = YNH4 * MWNH4
GNH3_High = TMASSNH3 - ANH4
if( GNH3 < 0.0 ) then
print *, " NEG GNH3", TWOSO4, ANH4, TMASSNH3
call CheckStop("NEG GNH3")
end if
! RETURN
goto 333
END IF
!...to get here, BB .lt. 0.0, CC .gt. 0.0 always
DD = SQRT( DISC )
XXQ = -0.5 * ( BB + SIGN ( 1.0, BB ) * DD )
!...Since both roots are positive, select smaller root.
XNO3 = MIN( XXQ / AA, CC / XXQ )
END IF
!2/25/99 IJA
AH2O_High = 1000.0 * WH2O
YNH4 = TWOSO4 + XNO3
ASO4_High = TSO4_HighA * MWSO4
!dsSAFE ANO3 = XNO3 * MWNO3
ANO3_High = min(XNO3 * MWNO3, TMASSHNO3 )
!dsSAFE ANH4 = YNH4 * MWNH4 ! ds should be safe as NH4/SO4 >2, but anyway:
ANH4_High = min(YNH4 * MWNH4, TMASSNH3 ) ! ds should be safe as NH4/SO4 >2, but anyway:
GNH3_High = TMASSNH3 - ANH4_High
GNO3_High = TMASSHNO3 - ANO3_High
! if( GNH3 < 0.0 .or. GNO3 < 0.0 ) then
! print *, " NEG GNH3 GNO3", TWOSO4, ANH4, TMASSNH3, ANO3, TMASSHNO3
! call CheckStop("NEG GNH3 GNO3")
! end if
! RETURN
goto 333
END IF
!...liquid phase containing completely neutralized sulfate and
!... some nitrate. Solve for composition and quantity.
MAS = TSO4_HighA / max(WH2O,FLOOR)
MAN = 0.0
XNO3 = 0.0
YNH4 = TWOSO4
PHIOLD = 1.0
!...Start loop for iteration
!...The assumption here is that all sulfate is ammonium sulfate,
!... and is supersaturated at lower relative humidities.
diff1=1.0/TOLER1
DO 1501 NNN = 1, 150
NITR = NNN
GASQD = GAMAAN * GAMAAN
WSQD = WH2O * WH2O
KW2 = KAN * WSQD / GASQD
AA = 1.0 - KW2
BB = TWOSO4 + KW2 * ( TNO3 + TNH4 - TWOSO4 )
CC = -KW2 * TNO3 * ( TNH4 - TWOSO4 )
!...This is a quadratic for XNO3 [MICROMOLES / M**3] of nitrate in solution.
DISC = BB * BB - 4.0 * AA * CC
if( DEBUG%EQUIB ) then
MAXNNN1 = NNN
!if( MAXNNN1 > 140) print "(a,i4,9es12.3)", "NNN1 ", NNN, DISC, TNO3, TNH4, TWOSO4
end if
!...Check for complex roots, retain inital values and RETURN
!2/25/99 IJA
!DS IF ( DISC < 0.0 ) THEN
!dsDS IF ( DISC < FLOOR ) THEN
IF ( DISC < 0.0 ) THEN
if( DEBUG%EQUIB .and. debug_flag ) print *, "MARS DISC NEG2 ", XNO3, WH2O, DISC
XNO3 = 0.0
AH2O_High = 1000.0 * WH2O
YNH4 = TWOSO4
GNO3_High = HNO3
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = NO3
!ANH4 = YNH4 * MWNH4
ANH4_High = min( YNH4 * MWNH4, TMASSNH3) ! ds added "min"
GNH3_High = TMASSNH3 - ANH4_High
!WRITE( 10, * ) ' COMPLEX ROOTS '
! RETURN
goto 333
END IF
! 2/25/99 IJA
! Deal with degenerate case (yoj)
!DS IF ( AA /= 0.0 ) THEN
IF ( abs(AA) > FLOOR ) THEN
if( DEBUG%EQUIB .and. debug_flag ) print "(a,9es11.3)", "MARS DEGEN ", XNO3, WH2O, DISC, AA, BB, CC
DD = SQRT( DISC )
XXQ = -0.5 * ( BB + SIGN ( 1.0, BB ) * DD )
RR1 = XXQ / AA
RR2 = CC / XXQ! Normally: BB>0, xxQ<0
!...choose minimum positve root
IF ( ( RR1 * RR2 ) < 0.0 ) THEN
if( DEBUG%EQUIB .and. debug_flag ) print "(a,10es10.3)", "MARS RR1*RR2 ", XNO3, WH2O, DISC, RR1, RR2
XNO3 = MAX( RR1, RR2 )
ELSE if(MIN( RR1, RR2 )>0.0)then
XNO3 = MIN( RR1, RR2 )
ELSE!two negative roots !DS PW added 4th July 2012
!--------------------- return copied from above
XNO3 = 0.0
AH2O_High = 1000.0 * WH2O
YNH4 = TWOSO4
GNO3_High = HNO3
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = NO3
!ds ANH4 = YNH4 * MWNH4
ANH4_High = min( YNH4 * MWNH4, TMASSNH3) ! ds added "min"
GNH3_High = TMASSNH3 - ANH4_High
if( DEBUG%EQUIB .and. debug_flag ) WRITE( *, * ) ' TWO NEG ROOTS '
! RETURN
goto 333
END IF
ELSE
XNO3 = - CC / BB ! AA equals zero here
if( DEBUG%EQUIB .and. debug_flag ) print "(a,4es10.3)", "MARS NONDEGEN ", AA, BB, CC, XNO3
END IF
XNO3 = MIN( XNO3, TNO3 )
!...This version assumes no solid sulfate forms (supersaturated )
!... Now update water
CALL AWATER ( fRH, TSO4_HighA, YNH4, XNO3, AH2O)
!...ZSR relationship is used to set water levels. Units are
!... 10**(-6) kg water/ (cubic meter of air)
!... The conversion from micromoles to moles is done by the units of WH2O.
WH2O = 1.0E-3 * AH2O
!...Ionic balance determines the ammonium in solution.
MAN = XNO3 / max(WH2O,FLOOR)
MAS = TSO4_HighA / max(WH2O,FLOOR)
MNH4 = 2.0 * MAS + MAN
YNH4 = MNH4 * WH2O
!st ... FIXING
if(MNH4<0. .or. MAS<0. .or. MAN<0.) then
MNH4 = 1.e-30
MAS = 1.e-30
MAN = 1.e-30
end if
!...MAS, MAN and MNH4 are the aqueous concentrations of sulfate, nitrate,
!... and ammonium in molal units (moles/(kg water) ).
STION = 3.0 * MAS + MAN
CAT( 1 ) = 0.0
CAT( 2 ) = MNH4
AN ( 1 ) = MAS
AN ( 2 ) = MAN
AN ( 3 ) = 0.0
CALL ACTCOF ( CAT, AN, GAMS, MOLNU, PHIBAR , ERRMARK,1,debug_flag)
GAMAAN = GAMS( 2, 2 )
!...Use GAMAAN for convergence control
EROR = ABS( GAMOLD - GAMAAN ) / GAMOLD
GAMOLD = GAMAAN
!...Check to see if we have a solution
IF ( EROR <= TOLER1 ) THEN
!!! WRITE( 11, * ) RH, STION, GAMS( 1, 1 ),GAMS( 1, 2 ), GAMS( 1, 3 ),
!!! & GAMS( 2, 1 ), GAMS( 2, 2 ), GAMS( 2, 3 ), PHIBAR
! 2/25/99 IJA
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = min(XNO3 * MWNO3,TMASSHNO3)
ANH4_High = min( YNH4 * MWNH4, TMASSNH3 )
GNO3_High = TMASSHNO3 - ANO3_High
GNH3_High = TMASSNH3 - ANH4_High
AH2O_High = 1000.0 * WH2O
! RETURN
goto 333
END IF
1501 CONTINUE
!...If after NITR iterations no solution is found, then:
! FSB retain the initial values of nitrate particle and vapor
! 2/25/99 IJA
ASO4_High = TSO4_HighA * MWSO4
ANO3_High = NO3
XNO3 = NO3 / MWNO3
YNH4 = TWOSO4
! ANH4 = YNH4 * MWNH4
ANH4_High = min( YNH4 * MWNH4, TMASSNH3 ) ! ds pw added "min"
CALL AWATER ( fRH, TSO4_HighA, YNH4, XNO3, AH2O_High)
GNO3_High = HNO3
GNH3_High = TMASSNH3 - ANH4_High
! RETURN
goto 333
333 CONTINUE !finished DO_RATIO_High_2
! ELSE
ENDIF
IF ( DO_RATIO_Low_2 ) THEN ! NH4/SO4 <= 2
!......................................
!......... Low Ammonia Case ...........
!......................................
!...coded by Dr. Francis S. Binkowski 12/8/91.(4/26/95)
! modified 8/28/98
!...All cases covered by this logic
WH2O = 0.0
CALL AWATER ( fRH, TSO4_LowA, TNH4, TNO3, AH2O )
WH2O = 1.0E-3 * AH2O
ZH2O = AH2O
!...convert 10**(-6) kg water/(cubic meter of air) to micrograms of water
!... per cubic meter of air (1000 g = 1 kg)
! 2/25/99 IJA
ASO4_Low = TSO4_LowA * MWSO4
ANH4_Low = TNH4 * MWNH4
!dsSAFE ANO3 = NO3
ANO3_Low = min( NO3, TMASSHNO3 )
GNO3_Low = TMASSHNO3 - ANO3_Low
GNH3_Low = FLOOR
AH2O_Low = 1.0E3 *WH2O
!...Check for zero water.
!ds IF ( WH2O == 0.0 ) RETURN
IF ( abs(WH2O) < FLOOR ) goto 111!RETURN
ZSO4 = TSO4_LowA / max(WH2O,FLOOR)
!...ZSO4 is the molality of total sulfate i.e. MSO4 + MHSO4
!!! IF ( ZSO4 > 11.0 ) THEN
!...do not solve for aerosol nitrate for total sulfate molality
!... greater than 11.0 because the model parameters break down
!... greater than 9.0 because the model parameters break down
IF ( ZSO4 > 9.0 ) THEN ! 18 June 97
goto 111
!RETURN
END IF
!...First solve with activity coeffs of 1.0, then iterate.
PHIOLD = 1.0
GAMANA = 1.0
GAMAS1 = 1.0
GAMAS2 = 1.0
GAMAAB = 1.0
GAMOLD = 1.0
!...All ammonia is considered to be aerosol ammonium.
MNH4 = TNH4 / max(WH2O,FLOOR)
!...MNH4 is the molality of ammonium ion.
YNH4 = TNH4