source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/build/ppsrc/INCA_SRC/mknoprod.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: 10.5 KB
Line 
1
2
3
4
5
6
7
8
9
10
11
12!$Id: mknoprod.F90 157 2010-01-18 14:21:59Z acosce $
13!! =========================================================================
14!! INCA - INteraction with Chemistry and Aerosols
15!!
16!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
17!!           Unite mixte CEA-CNRS-UVSQ
18!!
19!! Contributors to this INCA subroutine:
20!!
21!! Didier Hauglustaine, LSCE, hauglustaine@cea.fr
22!! Line Jourdain, SA
23!!
24!! Anne Cozic, LSCE, anne.cozic@cea.fr
25!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
26!!
27!! This software is a computer program whose purpose is to simulate the
28!! atmospheric gas phase and aerosol composition. The model is designed to be
29!! used within a transport model or a general circulation model. This version
30!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
31!! for emissions, transport (resolved and sub-grid scale), photochemical
32!! transformations, and scavenging (dry deposition and washout) of chemical
33!! species and aerosols interactively in the GCM. Several versions of the INCA
34!! model are currently used depending on the envisaged applications with the
35!! chemistry-climate model.
36!!
37!! This software is governed by the CeCILL  license under French law and
38!! abiding by the rules of distribution of free software.  You can  use,
39!! modify and/ or redistribute the software under the terms of the CeCILL
40!! license as circulated by CEA, CNRS and INRIA at the following URL
41!! "http://www.cecill.info".
42!!
43!! As a counterpart to the access to the source code and  rights to copy,
44!! modify and redistribute granted by the license, users are provided only
45!! with a limited warranty  and the software's author,  the holder of the
46!! economic rights,  and the successive licensors  have only  limited
47!! liability.
48!!
49!! In this respect, the user's attention is drawn to the risks associated
50!! with loading,  using,  modifying and/or developing or reproducing the
51!! software by the user in light of its specific status of free software,
52!! that may mean  that it is complicated to manipulate,  and  that  also
53!! therefore means  that it is reserved for developers  and  experienced
54!! professionals having in-depth computer knowledge. Users are therefore
55!! encouraged to load and test the software's suitability as regards their
56!! requirements in conditions enabling the security of their systems and/or
57!! data to be ensured and,  more generally, to use and operate it in the
58!! same conditions as regards security.
59!!
60!! The fact that you are presently reading this means that you have had
61!! knowledge of the CeCILL license and that you accept its terms.
62!! =========================================================================
63
64SUBROUTINE MKNOPROD(oro,  & 
65   lat,  &
66   lon,  &
67   area, &
68   pmid, &
69   zmid, &
70   temp, &
71   ctop, &
72   cbot, &
73   nx,   &
74   ny)
75  !----------------------------------------------------------------------
76  !     ... NOx from lightning
77  ! Line Jourdain, SA, 2001.
78  ! Modified, Didier Hauglustaine, IPSL,  08-2001.
79  !--------------------------------------------------------
80  USE LIGHTNING
81  USE CHEM_CONS
82  USE CHEM_MODS
83  USE SPECIES_NAMES
84  USE INCA_DIM
85
86  IMPLICIT NONE
87 
88  !--------------------------------------------------------
89  !       ... Dummy arguments
90  !--------------------------------------------------------
91  INTEGER, INTENT(in) :: ctop(PLON) !convective clouds top       
92  INTEGER, INTENT(in) :: cbot(PLON) !convective clouds bot
93  INTEGER, INTENT(in) :: nx, ny
94  REAL, INTENT(in)    :: oro(PLON)         !land-sea mask
95  REAL, INTENT(in)    :: area(PLON)        !Grid area (km2)
96  REAL, INTENT(in)    :: lat(PLON)         !Not used
97  REAL, INTENT(in)    :: pmid(PLON,PLEV)   !Pressure (Pa)
98  REAL, INTENT(in)    :: temp(PLON,PLEV)   !Temperature
99  REAL, INTENT(in)    :: lon(PLON)         !Not used
100  REAL, INTENT(in)    :: zmid(PLON,PLEV)   !Altitude (m)
101 
102  !--------------------------------------------------------
103  !       ... Local variables
104  !--------------------------------------------------------
105  INTEGER  :: i, k, l, lp
106  INTEGER  :: lmax, lmin
107 
108  REAL, PARAMETER :: secpyr = dayspy * 8.64e4
109  REAL, PARAMETER :: ap=0.021
110  REAL, PARAMETER :: bp=-0.648
111  REAL, PARAMETER :: cp=7.493
112  REAL, PARAMETER :: dp=-36.54
113  REAL, PARAMETER :: ep=63.09
114  REAL, PARAMETER :: zero = 1.e-20
115 
116  REAL     :: dx, dy
117  REAL     :: calibration
118  REAL     :: source
119  REAL     :: ziso0(PLON)
120  REAL     :: fic(PLON), fcg(PLON)
121  REAL     :: z(PLON)
122  REAL     :: a(PLON), b(PLON)
123  REAL     :: globoprod_no(PLON)
124  REAL     :: xprod_no(PLON,PLEV)
125  REAL     :: no_col(PLON)
126  REAL     :: factor(PLON)
127  REAL     :: scal(PLON)
128  REAL     :: altp(PLON,zdim)
129  REAL     :: atc(PLON,zdim-1), atm(PLON,zdim-1), amc(PLON,zdim-1)
130  REAL     :: btc(PLON,zdim-1), btm(PLON,zdim-1), bmc(PLON,zdim-1)
131  REAL     :: coeff_tcp(PLON,PLEV), coeff_mcp(PLON,PLEV), coeff_tmp(PLON,PLEV)
132  REAL     :: somcoeff1(PLON),somcoeff2(PLON),somcoeff3(PLON)
133  REAL     :: zinterf(PLON,PLEVP)
134  REAL     :: deltaz(PLON,PLEV)
135  REAL     :: zmidkm(PLON,PLEV)
136 
137  zmidkm = zmid*1.e-3
138
139  !     ... Zero all variables for this time-step
140  ztop=0.
141  flash=0.
142  ziso0=0.
143  dcold=0.
144  pcg=0. 
145  fcg=0.
146  fic=0.
147  globoprod_no=0.
148  xprod_no = 0.
149  factor=0. 
150  scal = 1.0
151
152  !     ... Ensure cloud presence
153  DO i=1,PLON
154    IF (pmid(i,ctop(i)) < pmid(i,cbot(i))) THEN
155        ztop(i)=zmidkm(i,ctop(i))
156    ENDIF
157  ENDDO
158
159  !     ... Flashes per sec (Price and Rind, 92)
160  DO i=1,PLON
161    IF (ztop(i) > 4.) THEN
162        flash(i)=3.44e-5*ztop(i)**4.90 *oro(i) + 6.40e-4*ztop(i)**1.73 *(1.-oro(i))
163    ENDIF
164  ENDDO
165 
166  !     ... Two calibrations applied:
167  !     1. To account for cloud height averaging (see Price and Rind, 1994)
168  !     2. To ensure 5 TgN in global and annual mean.
169 
170  dx = 2.*pi/FLOAT(nx)*r2d
171  dy = pi/FLOAT(ny)*r2d
172 
173  calibration = 9.7241e-1*EXP(4.8203e-2*dx*dy)
174  flash(:) = 2. * calibration * flash(:)   !KE
175
176  !     ... Find 0 celcius isotherm
177  DO i=1,PLON
178    DO l=1,PLEV-1
179      IF ( pmid(i,l) > 2.e4 ) THEN
180          IF ( (temp(i,l) > 273.15) .AND. (temp(i,l+1) < 273.15) ) THEN
181              a(i)=(temp(i,l+1)-temp(i,l))/(zmidkm(i,l+1)-zmidkm(i,l))
182              b(i)=temp(i,l)-a(i)*zmidkm(i,l)
183              ziso0(i)=(273.15-b(i))/a(i)
184          ENDIF
185      ENDIF
186    ENDDO
187  ENDDO
188 
189  !     ... Intracloud and cloud-to-ground lightning
190
191  DO i=1,PLON
192
193    IF (ziso0(i) < ztop(i)) THEN
194
195        dcold(i)=ztop(i)-ziso0(i)
196       
197        IF ( (dcold(i) > 5.5) .AND. (dcold(i) < 14.) ) THEN
198            z(i)=(ap*dcold(i)**4+bp*dcold(i)**3+cp*dcold(i)**2+dp*dcold(i)+ep)
199            pcg(i)=1./(z(i)+1.)
200        ELSE IF (dcold(i) <= 5.5) THEN
201            pcg(i)=0.
202        ELSE IF (dcold(i) >= 14.) THEN
203            pcg(i)=2.e-2
204        ENDIF
205       
206        fcg(i) = pcg(i)*flash(i)
207        fic(i) = (1.-pcg(i))*flash(i)
208
209    ENDIF
210   
211  ENDDO
212
213!     ... Change units to fl/s/m2 DH.
214
215      flash(:) = flash(:)/60./area(:)
216      flpcg(:) = fcg(:)/60./area(:)
217 
218  !     ... Column integrated NO production
219  globoprod_no(:) = 6.7e25 * (10.*fcg(:)+fic(:)) /60./area(:)
220 
221  !     ... Pickering et al. (1998) vertical distribution
222  DO i=1,PLON
223    xprod_no(i,1:ctop(i))=1.
224  ENDDO
225 
226  zinterf(:,1)=0.
227  zinterf(:,2)=2.*zmidkm(:,1)
228  DO l=2,PLEV
229    zinterf(:,l+1)=zinterf(:,l)+2.*(zmidkm(:,l)-zinterf(:,l))
230  ENDDO
231 
232  DO l=1,PLEV
233    deltaz(:,l)=(zinterf(:,l+1)-zinterf(:,l))*1.e3
234  ENDDO
235 
236     
237  DO i=1,PLON
238    IF (flash(i) > zero) THEN
239       
240        IF ( ABS(lat(i)) >= 30.) THEN
241            scal(i)=ztop(i)/14.5
242        ELSE
243            scal(i)=ztop(i)/16.
244        ENDIF
245       
246    ENDIF
247  ENDDO
248 
249  DO l=1,zdim
250    altp(:,l)=alt(l)*scal(:)
251  ENDDO
252 
253  DO i=1,PLON
254    DO l=1,zdim-1
255
256      IF (scal(i) > zero) THEN
257          atc(i,l)=(coeff_tc(l+1)-coeff_tc(l))/(altp(i,l+1)-altp(i,l))
258          btc(i,l)=coeff_tc(l)-atc(i,l)*altp(i,l)
259          atm(i,l)=(coeff_tm(l+1)-coeff_tm(l))/(altp(i,l+1)-altp(i,l))
260          btm(i,l)=coeff_tm(l)-atm(i,l)*altp(i,l)
261          amc(i,l)=(coeff_mc(l+1)-coeff_mc(l))/(altp(i,l+1)-altp(i,l))
262          bmc(i,l)=coeff_mc(l)-amc(i,l)*altp(i,l)
263      ENDIF
264     
265    ENDDO
266  ENDDO
267 
268  coeff_mcp=0.
269  coeff_tmp=0.
270  coeff_tcp=0.
271 
272  DO i=1,PLON
273    DO lp=1,zdim-1
274      DO l=1,ctop(i)
275       
276        IF (flash(i) > zero) THEN
277           
278            IF (zmidkm(i,l) <= altp(i,1)) THEN
279                coeff_tcp(i,l)=coeff_tc(1)
280                coeff_tmp(i,l)=coeff_tm(1)
281                coeff_mcp(i,l)=coeff_mc(1)
282            ENDIF
283           
284            IF ((zmidkm(i,l) <= altp(i,lp+1)) .AND. (zmidkm(i,l) >= altp(i,lp))) THEN
285                coeff_tcp(i,l)=atc(i,lp)*zmidkm(i,l)+btc(i,lp)
286                coeff_tmp(i,l)=atm(i,lp)*zmidkm(i,l)+btm(i,lp)
287                coeff_mcp(i,l)=amc(i,lp)*zmidkm(i,l)+bmc(i,lp)
288            ENDIF
289        ENDIF
290       
291      ENDDO
292    ENDDO
293  ENDDO
294 
295  somcoeff1=0.
296  somcoeff2=0.
297  somcoeff3=0.
298 
299  DO i=1,PLON
300    DO l=1,ctop(i)
301      somcoeff1(i)=somcoeff1(i)+coeff_tcp(i,l)*deltaz(i,l)
302      somcoeff2(i)=somcoeff2(i)+coeff_tmp(i,l)*deltaz(i,l)
303      somcoeff3(i)=somcoeff3(i)+coeff_mcp(i,l)*deltaz(i,l)
304    ENDDO
305  ENDDO
306 
307  DO i=1,PLON
308    DO l=1,ctop(i)
309     
310      IF (somcoeff1(i) > zero) THEN
311          coeff_tcp(i,l)=coeff_tcp(i,l)/somcoeff1(i)*1.e2
312      ENDIF
313     
314      IF (somcoeff2(i) > zero) THEN
315          coeff_tmp(i,l)=coeff_tmp(i,l)/somcoeff2(i)*1.e2
316      ENDIF
317     
318      IF (somcoeff3(i) > zero) THEN
319          coeff_mcp(i,l)=coeff_mcp(i,l)/somcoeff3(i)*1.e2
320      ENDIF
321     
322    ENDDO
323  ENDDO
324 
325  DO i=1,PLON
326    DO l=1,ctop(i)
327      IF (flash(i) > zero) THEN
328          IF ( ABS(lat(i)) <= 30. ) THEN
329             
330              xprod_no(i,l)=coeff_tcp(i,l)*0.01 * oro(i) + coeff_tmp(i,l)*0.01 * (1.-oro(i))
331          ELSE
332             
333              xprod_no(i,l)=coeff_mcp(i,l)*0.01
334          ENDIF
335      ENDIF
336    ENDDO
337  ENDDO
338 
339  no_col=0.
340  DO i=1,PLON
341    lmax=ctop(i)
342    DO l=1,lmax
343      no_col(i) = no_col(i) + xprod_no(i,l) * deltaz(i,l) * 1.e6
344    ENDDO
345  ENDDO
346 
347  prod_light=0.
348  factor=0.
349  DO i=1,PLON
350    lmax=ctop(i)
351    IF (no_col(i) /= 0.) THEN
352        factor(i) = globoprod_no(i) / no_col(i)
353    ENDIF
354    DO l=1,lmax
355      prod_light(i,l) = xprod_no(i,l) * factor(i)
356    ENDDO
357  ENDDO
358 
359  source = 0.
360  DO i=1,PLON
361    DO l=1,PLEV
362      source = source + prod_light(i,l)*area(i)*deltaz(i,l)*1.e6
363    ENDDO
364  ENDDO
365  source = source /6.02e23 *14. *secpyr *1.e-12
366
367!Integrated lightning NOx production (should be in kg/m2/s since prod_light is in molec/cm3/s)
368  DO i=1,PLON
369    DO l=1,PLEV
370      prod_light_col(i) = prod_light_col(i) + prod_light(i,l)*1.e6*deltaz(i,l)/6.02e23*14.*1.e-3
371    ENDDO
372  ENDDO
373 
374END SUBROUTINE MKNOPROD
375
Note: See TracBrowser for help on using the repository browser.