source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/exp_slv.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: 11.3 KB
Line 
1!$Id: exp_slv.F90 123 2009-03-27 10:38:52Z 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!! Stacy Walters, NCAR, stacy@ucar.edu
12!!
13!! Anne Cozic, LSCE, anne.cozic@cea.fr
14!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
15!!
16!! This software is a computer program whose purpose is to simulate the
17!! atmospheric gas phase and aerosol composition. The model is designed to be
18!! used within a transport model or a general circulation model. This version
19!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
20!! for emissions, transport (resolved and sub-grid scale), photochemical
21!! transformations, and scavenging (dry deposition and washout) of chemical
22!! species and aerosols interactively in the GCM. Several versions of the INCA
23!! model are currently used depending on the envisaged applications with the
24!! chemistry-climate model.
25!!
26!! This software is governed by the CeCILL  license under French law and
27!! abiding by the rules of distribution of free software.  You can  use,
28!! modify and/ or redistribute the software under the terms of the CeCILL
29!! license as circulated by CEA, CNRS and INRIA at the following URL
30!! "http://www.cecill.info".
31!!
32!! As a counterpart to the access to the source code and  rights to copy,
33!! modify and redistribute granted by the license, users are provided only
34!! with a limited warranty  and the software's author,  the holder of the
35!! economic rights,  and the successive licensors  have only  limited
36!! liability.
37!!
38!! In this respect, the user's attention is drawn to the risks associated
39!! with loading,  using,  modifying and/or developing or reproducing the
40!! software by the user in light of its specific status of free software,
41!! that may mean  that it is complicated to manipulate,  and  that  also
42!! therefore means  that it is reserved for developers  and  experienced
43!! professionals having in-depth computer knowledge. Users are therefore
44!! encouraged to load and test the software's suitability as regards their
45!! requirements in conditions enabling the security of their systems and/or
46!! data to be ensured and,  more generally, to use and operate it in the
47!! same conditions as regards security.
48!!
49!! The fact that you are presently reading this means that you have had
50!! knowledge of the CeCILL license and that you accept its terms.
51!! =========================================================================
52
53#include <inca_define.h>
54
55SUBROUTINE EXP_SOL( &
56   base_sol        ,&
57   reaction_rates  ,&
58   het_rates       ,&
59   extfrc          ,& 
60#ifdef NMHC
61   co_prod         ,&
62   co_loss         ,&
63   ch4_loss        ,&
64   n2o_loss        ,&
65   hnm             ,&
66#endif
67   nstep           ,&
68   delt )
69
70  !-----------------------------------------------------------------------
71  !             ... Exp_sol advances the volumetric mixing ratio
72  !           forward one time step via the fully explicit
73  !           Euler scheme
74  ! Stacy Walters, NCAR, 1998.
75  ! Modified by Didier Hauglustaine, IPSL, for LMDz/INCA, 2000.
76  !-----------------------------------------------------------------------
77
78  USE INCA_DIM
79  USE CHEM_MODS, ONLY: clsmap
80  USE SPECIES_NAMES
81
82  USE RATE_INDEX_MOD
83
84  IMPLICIT NONE
85
86  !-----------------------------------------------------------------------
87  !             ... Dummy arguments
88  !-----------------------------------------------------------------------
89  INTEGER, INTENT(in) ::  nstep            ! time step index
90  REAL, INTENT(in)    ::  delt             ! time step in seconds
91  REAL, INTENT(in)    ::  reaction_rates(PLNPLV,RXNCNT)
92#ifdef NMHC
93  REAL, INTENT(in)    ::  hnm(PLNPLV)
94  REAL, INTENT(out)   ::  co_prod(PLNPLV), co_loss(PLNPLV)
95  REAL, INTENT(out)   ::  ch4_loss(PLNPLV)
96  REAL, INTENT(out)   ::  n2o_loss(PLNPLV)
97#endif
98# if HETCNT != 0
99  REAL, INTENT(in)   ::  het_rates(PLNPLV,HETCNT)
100# else
101  REAL, INTENT(in)   ::  het_rates(1)
102# endif
103# if EXTCNT != 0
104  REAL, INTENT(out)   ::  extfrc(PLNPLV,EXTCNT)
105# else
106  REAL, INTENT(out)   ::  extfrc(1)
107# endif
108  REAL, INTENT(inout) ::  base_sol(PLNPLV,PCNST)
109
110  !-----------------------------------------------------------------------
111  !             ... Local variables
112  !-----------------------------------------------------------------------
113  INTEGER  ::  k, l, m
114  REAL     ::  prod(PLNPLV,CLSCNT1)
115  REAL     ::  loss(PLNPLV,CLSCNT1)
116# if CLSINDPRD1 != 0 || EXTCNT != 0
117  REAL     ::  ind_prd(PLNPLV,CLSCNT1)
118# endif
119
120
121# if CLSINDPRD1 != 0 || EXTCNT != 0
122  !-----------------------------------------------------------------------     
123  !        ... Put "independent" production in the forcing
124  !-----------------------------------------------------------------------     
125  CALL INDPRD( &
126     1          ,&
127     ind_prd    ,&
128     base_sol   ,&
129     extfrc     ,&
130     reaction_rates )
131
132# endif
133  !-----------------------------------------------------------------------     
134  !             ... Form F(y)
135  !-----------------------------------------------------------------------     
136  CALL EXP_PROD_LOSS( &
137     prod            ,&
138     loss            ,&
139     base_sol        ,&
140     reaction_rates  ,&
141     het_rates )
142
143# if CLSINDPRD1 != 0 || EXTCNT != 0
144  DO k = 1,CLSCNT1
145    prod(:,k) = prod(:,k) + ind_prd(:,k)
146  END DO
147# endif
148
149#if defined(AERONLY) || defined(GES)
150  !-----------------------------------------------------------------------     
151  !     ... Solve for the mixing ratio at t(n+1)
152  !-----------------------------------------------------------------------     
153  DO k = 1,CLSCNT1
154    l = clsmap(k,1)
155    base_sol(:,l) = base_sol(:,l)+delt*(prod(:,k)-loss(:,k))
156  END DO
157 
158#else
159  !-----------------------------------------------------------------------     
160  !     ... Solve for the mixing ratio at t(n+1)
161  !-----------------------------------------------------------------------     
162  ! Warning: CLSCNT1 is O3S and CLSCNT1-1 is O3I
163  DO k = 1,CLSCNT1 - 2
164    l = clsmap(k,1)
165    base_sol(:,l) = base_sol(:,l)+delt*(prod(:,k)-loss(:,k))
166  END DO
167 
168  !-----------------------------------------------------------------------
169  !       ... Special code for O3S;
170  !           note O3S is assumed to be the last "explicit" species
171  !-----------------------------------------------------------------------
172  l = clsmap(CLSCNT1,1)
173
174
175#ifdef NMHC
176  !     reactions included in the loss processes of O3S
177  !     defined by
178  !     O3S-loss = J1 * ( k3 / (k1 + k2 + k3 + k4 + k5 + k6)) +
179  !                     (k7 + k8 + k9 + k10 + k11 + k12 + k13 + k14 + k15 + k16)
180  !     where
181  !     J1:  O3 + hv --> O1D + O2        ( reaction_rates(:,2) )
182  !     k1:  O1D + N2                    ( reaction_rates(:,110) )
183  !     k2:  O1D + O2                    ( reaction_rates(:,111) )
184  !     k3:  O1D + H2O                   ( reaction_rates(:,112) )
185  !     k4:  O1D + H2                    ( reaction_rates(:,113) )
186  !     k5:  O1D + N2O --> 2*NO          ( reaction_rates(:,114) )
187  !     k6:  O1D + N2O --> N2 + O2       ( reaction_rates(:,115) )
188  !     k7:  O3 + OH                     ( reaction_rates(:,126) )
189  !     k8:  O3 + HO2                    ( reaction_rates(:,127) )
190  !     k9:  O3 + H                      ( reaction_rates(:,134) )
191  !    k10:  O3 + C3H6                   ( reaction_rates(:,177) - reaction_rates(:,179))
192  !    k11:  O3 + C2H4                   ( reaction_rates(:,180) - reaction_rates(:,181))
193  !    k12:  O3 + ISOP                   ( reaction_rates(:,198) - reaction_rates(:,200))
194  !    k13:  O3 + APIN                   ( reaction_rates(:,221))
195  !    k14:  O3 + MACR                   ( reaction_rates(:,240) - reaction_rates(:,242))
196  !    k15:  O3 + MVK                    ( reaction_rates(:,244) - reaction_rates(:,246))
197  !    k16:  O3 + ALKEN                  ( reaction_rates(:,297) - reaction_rates(:,299))
198#ifdef AER 
199  loss(:,1) = reaction_rates(:,jO3O1D)*reaction_rates(:,kO1DH2O)   &
200     /(reaction_rates(:,kO1DN2) + reaction_rates(:,kO1DO2)   &
201     + reaction_rates(:,kO1DH2O) + reaction_rates(:,kO1DH2)   &
202     + reaction_rates(:,kN2OO1D2NO) + reaction_rates(:,kN2OO1DN2))  &
203     + reaction_rates(:,kOHO3)*base_sol(:,id_oh)               &
204     + reaction_rates(:,kHO2O3)*base_sol(:,id_ho2)              &
205     + reaction_rates(:,kHO3)*base_sol(:,id_h)                &
206     + reaction_rates(:,kC3H6O31)*base_sol(:,id_c3h6)             &
207     + reaction_rates(:,kC3H6O32)*base_sol(:,id_c3h6)             &
208     + reaction_rates(:,kC3H6O33)*base_sol(:,id_c3h6)             &
209     + reaction_rates(:,kC2H4O31)*base_sol(:,id_c2h4)             &
210     + reaction_rates(:,kC2H4O32)*base_sol(:,id_c2h4)             &
211     + reaction_rates(:,kISOPO31)*base_sol(:,id_isop)             &
212     + reaction_rates(:,kISOPO32)*base_sol(:,id_isop)             &
213     + reaction_rates(:,kISOPO33)*base_sol(:,id_isop)             &
214     + reaction_rates(:,kAPINO33)*base_sol(:,id_apin)             &
215     + reaction_rates(:,kMACRO31)*base_sol(:,id_macr)             &
216     + reaction_rates(:,kMACRO32)*base_sol(:,id_macr)             &
217     + reaction_rates(:,kMVKO31)*base_sol(:,id_mvk)              &
218     + reaction_rates(:,kMVKO32)*base_sol(:,id_mvk)              &
219     + reaction_rates(:,kALKENO31)*base_sol(:,id_alken)            &
220     + reaction_rates(:,kALKENO32)*base_sol(:,id_alken)            &
221     + reaction_rates(:,kALKENO33)*base_sol(:,id_alken)
222#else
223  loss(:,1) = reaction_rates(:,jO3O1D)*reaction_rates(:,kO1DH2O)   &
224     /(reaction_rates(:,kO1DN2) + reaction_rates(:,kO1DO2)   &
225     + reaction_rates(:,kO1DH2O) + reaction_rates(:,kO1DH2)   &
226     + reaction_rates(:,kN2OO1D2NO) + reaction_rates(:,kN2OO1DN2))  &
227     + reaction_rates(:,kOHO3)*base_sol(:,id_oh)               &
228     + reaction_rates(:,kHO2O3)*base_sol(:,id_ho2)              &
229     + reaction_rates(:,kHO3)*base_sol(:,id_h)                &
230     + reaction_rates(:,kC3H6O31)*base_sol(:,id_c3h6)             &
231     + reaction_rates(:,kC3H6O32)*base_sol(:,id_c3h6)             &
232     + reaction_rates(:,kC3H6O33)*base_sol(:,id_c3h6)             &
233     + reaction_rates(:,kC2H4O31)*base_sol(:,id_c2h4)             &
234     + reaction_rates(:,kC2H4O32)*base_sol(:,id_c2h4)             &
235     + reaction_rates(:,kISOPO31)*base_sol(:,id_isop)             &
236     + reaction_rates(:,kISOPO32)*base_sol(:,id_isop)             &
237     + reaction_rates(:,kISOPO33)*base_sol(:,id_isop)             &
238     + reaction_rates(:,kAPINO3)*base_sol(:,id_apin)             &
239     + reaction_rates(:,kMACRO31)*base_sol(:,id_macr)             &
240     + reaction_rates(:,kMACRO32)*base_sol(:,id_macr)             &
241     + reaction_rates(:,kMVKO31)*base_sol(:,id_mvk)              &
242     + reaction_rates(:,kMVKO32)*base_sol(:,id_mvk)              &
243     + reaction_rates(:,kALKENO31)*base_sol(:,id_alken)            &
244     + reaction_rates(:,kALKENO32)*base_sol(:,id_alken)            &
245     + reaction_rates(:,kALKENO33)*base_sol(:,id_alken)
246
247#endif
248
249#endif
250   ! ... Special code for O3S;
251  base_sol(:,l) = base_sol(:,l)*EXP( -delt*loss(:,1) )
252
253#ifdef NMHC
254  !     diagnostics for production and destruction terms     
255  co_prod(:)  = prod(:,8) * hnm(:) 
256  co_loss(:)  = loss(:,8) * hnm(:) 
257  ch4_loss(:) = loss(:,6) * hnm(:) 
258  n2o_loss(:) = loss(:,7) * hnm(:) 
259#endif
260
261#endif
262
263END SUBROUTINE EXP_SOL
Note: See TracBrowser for help on using the repository browser.