source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/nitrates.F90 @ 6610

Last change on this file since 6610 was 6610, checked in by acosce, 10 months ago

INCA used for ICOLMDZORINCA_CO2_Transport_GMD_2023

File size: 8.8 KB
Line 
1!$Id: nitrates.F90 104 2008-12-23 10:28:51Z acosce $
2!! =========================================================================
3!! INCA - INteraction with Chemistry and Aerosols
4!!
5!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
6!!           Unite mixte CEA-CNRS-UVSQ
7!!
8!! Contributors to this INCA subroutine:
9!!
10!! Didier Hauglustaine, LSCE, hauglustaine@cea.fr
11!!
12!! This software is a computer program whose purpose is to simulate the
13!! atmospheric gas phase and aerosol composition. The model is designed to be
14!! used within a transport model or a general circulation model. This version
15!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
16!! for emissions, transport (resolved and sub-grid scale), photochemical
17!! transformations, and scavenging (dry deposition and washout) of chemical
18!! species and aerosols interactively in the GCM. Several versions of the INCA
19!! model are currently used depending on the envisaged applications with the
20!! chemistry-climate model.
21!!
22!! This software is governed by the CeCILL  license under French law and
23!! abiding by the rules of distribution of free software.  You can  use,
24!! modify and/ or redistribute the software under the terms of the CeCILL
25!! license as circulated by CEA, CNRS and INRIA at the following URL
26!! "http://www.cecill.info".
27!!
28!! As a counterpart to the access to the source code and  rights to copy,
29!! modify and redistribute granted by the license, users are provided only
30!! with a limited warranty  and the software's author,  the holder of the
31!! economic rights,  and the successive licensors  have only  limited
32!! liability.
33!!
34!! In this respect, the user's attention is drawn to the risks associated
35!! with loading,  using,  modifying and/or developing or reproducing the
36!! software by the user in light of its specific status of free software,
37!! that may mean  that it is complicated to manipulate,  and  that  also
38!! therefore means  that it is reserved for developers  and  experienced
39!! professionals having in-depth computer knowledge. Users are therefore
40!! encouraged to load and test the software's suitability as regards their
41!! requirements in conditions enabling the security of their systems and/or
42!! data to be ensured and,  more generally, to use and operate it in the
43!! same conditions as regards security.
44!!
45!! The fact that you are presently reading this means that you have had
46!! knowledge of the CeCILL license and that you accept its terms.
47!! =========================================================================
48
49#include <inca_define.h>
50
51#ifdef AER
52#ifndef DUSS
53SUBROUTINE AERTHERM (&
54   delt             ,&
55   temp             ,&
56   relhum           ,&
57   pmid             ,&
58   hnm              ,&
59   asno3m_p_nh3hno3 ,&
60   asnh4m_p_nh3hno3 ,&
61   hno3_p_nh3hno3   ,&
62   nh3_p_nh3hno3    ,&
63   mmr              ,&
64   vmr ) 
65
66  USE INCA_DIM
67  USE CONST_MOD
68  USE AEROSOL_DIAG
69  USE AEROSOL_MOD
70  USE SPECIES_NAMES
71  USE CHEM_MODS, ONLY : invariants
72  USE PRINT_INCA
73  USE RATE_INDEX_MOD
74
75  IMPLICIT NONE
76 
77  !-----------------------------------------------------------------
78  !        ... Dummy arguments
79  !-----------------------------------------------------------------
80  REAL, INTENT(in)    ::  delt                        ! timestep in seconds
81  REAL, INTENT(in)    ::  temp(PLON,PLEV)             ! temperature
82  REAL, INTENT(in)    ::  pmid(PLON,PLEV)             ! midpoint pressure in Pa
83  REAL, INTENT(in)    ::  hnm(PLON,PLEV)              ! total concentration
84  REAL, INTENT(inout) ::  vmr(PLON,PLEV,PCNST)        ! xported species ( vmr )
85  REAL, INTENT(inout) ::  mmr(PLON,PLEV,PCNST)        ! xported species ( mmr )
86  REAL, INTENT(inout) ::  relhum(PLON,PLEV)           ! relative humidity
87  REAL, INTENT(inout) ::  asno3m_p_nh3hno3(PLON,PLEV) ! for diagnostics
88  REAL, INTENT(inout) ::  asnh4m_p_nh3hno3(PLON,PLEV) ! for diagnostics
89  REAL, INTENT(inout) ::  hno3_p_nh3hno3(PLON,PLEV)   ! for diagnostics
90  REAL, INTENT(inout) ::  nh3_p_nh3hno3(PLON,PLEV)    ! for diagnostics
91     
92!-----------------------------------------------------------------
93!        ... Local variables
94!-----------------------------------------------------------------
95  INTEGER  ::  i, j, k
96
97  REAL, DIMENSION(PLON,PLEV) :: tinv
98  REAL, DIMENSION(PLON,PLEV) :: relhumloc, relhum1
99  REAL, DIMENSION(PLON,PLEV) :: hno3, nh3, nh4p, no3m, so42m, nh4pini
100  REAL, DIMENSION(PLON,PLEV) :: tn, ta, tadisp, taini, tnta, ts
101  REAL, DIMENSION(PLON,PLEV) :: tam, tsm
102  REAL, DIMENSION(PLON,PLEV) :: drh !0 to 1 as relhum
103  REAL, DIMENSION(PLON,PLEV) :: kps, kpl, kpl1, kpl2, kpl3
104  REAL, DIMENSION(PLON,PLEV) :: zrho
105  REAL, DIMENSION(PLON,PLEV) :: vmr0_no3m,vmr0_nh4p,vmr0_nh3,vmr0_hno3 
106
107  REAL, PARAMETER :: mwnh3  = 17.e-3 !Kg/mol
108  REAL, PARAMETER :: mwhno3 = 63.e-3
109  REAL, PARAMETER :: mwnh4  = 18.e-3
110  REAL, PARAMETER :: mwno3  = 62.e-3
111  REAL, PARAMETER :: mwso4  = 96.e-3
112  REAL, PARAMETER :: mwa    = 29.e-3
113
114  REAL :: zso4
115  REAL :: wrk1, wrk2
116  REAL :: tfac, tautot, hno3eq, nh3eq
117
118  zrho(:,:) = pmid(:,:)/(temp(:,:)*287.04)
119  tinv(:,:)    = 1. / temp(:,:) 
120
121  relhumloc = relhum
122  WHERE( relhumloc < 0.e0 )
123      relhumloc = 0.e0
124  END WHERE
125  WHERE( relhumloc > 0.98)
126      relhumloc = 0.98
127  END WHERE
128  relhum1(:,:) = 1. - relhumloc(:,:) 
129
130!Equilibrium constant based on Mozurkewich, 1993
131  drh(:,:)  = EXP(723.7*tinv(:,:)+1.6954) * 1.e-2
132  kps(:,:)  = EXP(118.87-24084.*tinv(:,:)-6.025*LOG(temp(:,:)))
133  kpl1(:,:) = EXP(-135.94+8763.*tinv(:,:)+19.12*LOG(temp(:,:)))
134  kpl2(:,:) = EXP(-122.65+9969.*tinv(:,:)+16.22*LOG(temp(:,:)))
135  kpl3(:,:) = EXP(-182.61+13875.*tinv(:,:)+24.46*LOG(temp(:,:)))
136
137  DO i  = 1, PLON
138    DO j  = 1, PLEV
139
140      kpl(i,j) = kps(i,j)
141
142      IF ( relhumloc(i,j) >= drh(i,j) ) THEN
143        kpl(i,j) = (kpl1(i,j)-kpl2(i,j)*relhum1(i,j)+kpl3(i,j)*relhum1(i,j)**2.)&
144                 * relhum1(i,j)**1.75*kpl(i,j)
145      ENDIF
146
147    ENDDO
148  ENDDO
149 
150!Initial mixing ratio for diagnostics purpose
151  vmr0_no3m(:,:) = vmr(:,:,id_ASNO3M)
152  vmr0_nh4p(:,:) = vmr(:,:,id_ASNH4M)
153  vmr0_hno3(:,:) = vmr(:,:,id_HNO3)
154  vmr0_nh3(:,:)  = vmr(:,:,id_NH3)
155
156!Volume mixing ratios
157  hno3(:,:)  = vmr(:,:,id_HNO3)*1.e9
158
159  no3m(:,:)  = vmr(:,:,id_ASNO3M)*1.e9
160  nh4p(:,:)  = vmr(:,:,id_ASNH4M)*1.e9
161  so42m(:,:) = vmr(:,:,id_ASSO4M)*1.e9
162  nh3(:,:)   = vmr(:,:,id_NH3)*1.e9
163
164!Total in moles
165  tn(:,:)    = hno3(:,:) + no3m(:,:)
166  ta(:,:)    = nh3(:,:)  + nh4p(:,:)
167  ts(:,:)    = so42m(:,:)
168
169  tam(:,:)   = mmr(:,:,id_ASNH4M)*1.e9 + mmr(:,:,id_NH3)*1.e9
170  tsm(:,:)   = mmr(:,:,id_ASSO4M)*1.e9
171
172  DO i  = 1, PLON
173    DO j  = 1, PLEV
174
175!Sulfate state. Metzger et al. 2002
176
177    zso4 = 2.0
178    IF (tsm(i,j)> 0.5*tam(i,j)) zso4 = 1.5
179    IF (tsm(i,j)>     tam(i,j)) zso4 = 1.0
180    so4state(i,j) = zso4
181
182    tadisp(i,j)  = MAX((ta(i,j)-zso4*so42m(i,j)),0.)
183    tnta(i,j)    = tn(i,j)*tadisp(i,j)
184    nh4pini(i,j) = MIN(ta(i,j),zso4*so42m(i,j))
185
186!Step 1: Equilibrium concentrations
187      IF   ( tnta(i,j) > kpl(i,j) ) THEN
188
189        wrk1    = (tadisp(i,j)+tn(i,j))**2. - 4.*(tnta(i,j)-kpl(i,j))
190        wrk1    = MAX(wrk1,0.)
191        wrk2    = 0.5*(tadisp(i,j)+tn(i,j)-SQRT(wrk1))
192        wrk2    = MAX(wrk2,0.)
193        wrk2    = MIN(wrk2,tn(i,j))
194        wrk2    = MIN(wrk2,tadisp(i,j))
195
196        hno3eq  = MAX(tn(i,j)-wrk2,0.)
197        nh3eq   = MAX(tadisp(i,j)-wrk2,0.)
198
199!Step 2a: Time dependence -time constants (tau and delt in sec)
200!This needs to be calculated explicitely based on Wexler and Seinfeld (1990) and Ackermann et al(1995)
201!for now use the value provided by Ackermann et al. for alpha=0.5 and Radius=0.1um. This should be calculated explicitely but little effect.
202
203        tautot = 2.05*60.
204        tfac = 1.-EXP(-delt/tautot)
205
206!Step 2b: Time dependence -gas phase concentrations
207       
208        hno3(i,j) = MAX((hno3(i,j) - tfac * (hno3(i,j)-hno3eq)),0.)
209        nh3(i,j)  = MAX((nh3(i,j)  - tfac * (nh3(i,j)-nh3eq)),  0.)
210
211!Step 3: Update aerosol phase
212
213        no3m(i,j) = MAX((tn(i,j)-hno3(i,j)),0.)
214        nh4p(i,j) = MAX((tadisp(i,j)-nh3(i,j)),0.)
215
216        vmr(i,j,id_HNO3)   = hno3(i,j) * 1.e-9
217        vmr(i,j,id_NH3)    = nh3(i,j)  * 1.e-9
218        vmr(i,j,id_ASNO3M) = no3m(i,j) * 1.e-9
219        vmr(i,j,id_ASNH4M) = nh4p(i,j) * 1.e-9
220
221      ELSE
222
223        vmr(i,j,id_HNO3)   = tn(i,j)     * 1.e-9
224        vmr(i,j,id_NH3)    = tadisp(i,j) * 1.e-9
225        vmr(i,j,id_ASNH4M) = 1.e-19
226        vmr(i,j,id_ASNO3M) = 1.e-19
227
228      ENDIF
229
230!Add back the ammonium sulfate
231 
232      vmr(i,j,id_ASNH4M) = vmr(i,j,id_ASNH4M) + nh4pini(i,j) * 1.e-9
233
234    END DO
235  END DO
236
237!Store changes for diagnostics (molec/cm3/s)
238  asno3m_p_nh3hno3(:,:) = (vmr(:,:,id_ASNO3M)-vmr0_no3m(:,:)) * hnm(:,:)/delt
239  asnh4m_p_nh3hno3(:,:) = (vmr(:,:,id_ASNH4M)-vmr0_nh4p(:,:)) * hnm(:,:)/delt
240#ifdef NMHC
241  hno3_p_nh3hno3(:,:)   = (vmr(:,:,id_HNO3)-vmr0_hno3(:,:))   * hnm(:,:)/delt
242#endif
243  nh3_p_nh3hno3(:,:)    = (vmr(:,:,id_NH3)-vmr0_nh3(:,:))     * hnm(:,:)/delt
244 
245END SUBROUTINE AERTHERM
246#endif
247#endif
Note: See TracBrowser for help on using the repository browser.