Context Navigation

source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/build/ppsrc/INCA_SRC/chemmain.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: 18.9 KB

Line
1
2
3
4
5
6
7
8
9
10
11
12	!$Id: chemmain.F90 125 2009-04-21 10:22:57Z 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	!! Stacy Walters, NCAR, stacy@ucar.edu
23	!! C. Textor, LSCE
24	!!
25	!! Anne Cozic, LSCE, anne.cozic@cea.fr
26	!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
27	!!
28	!! This software is a computer program whose purpose is to simulate the
29	!! atmospheric gas phase and aerosol composition. The model is designed to be
30	!! used within a transport model or a general circulation model. This version
31	!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
32	!! for emissions, transport (resolved and sub-grid scale), photochemical
33	!! transformations, and scavenging (dry deposition and washout) of chemical
34	!! species and aerosols interactively in the GCM. Several versions of the INCA
35	!! model are currently used depending on the envisaged applications with the
36	!! chemistry-climate model.
37	!!
38	!! This software is governed by the CeCILL license under French law and
39	!! abiding by the rules of distribution of free software. You can use,
40	!! modify and/ or redistribute the software under the terms of the CeCILL
41	!! license as circulated by CEA, CNRS and INRIA at the following URL
42	!! "http://www.cecill.info".
43	!!
44	!! As a counterpart to the access to the source code and rights to copy,
45	!! modify and redistribute granted by the license, users are provided only
46	!! with a limited warranty and the software's author, the holder of the
47	!! economic rights, and the successive licensors have only limited
48	!! liability.
49	!!
50	!! In this respect, the user's attention is drawn to the risks associated
51	!! with loading, using, modifying and/or developing or reproducing the
52	!! software by the user in light of its specific status of free software,
53	!! that may mean that it is complicated to manipulate, and that also
54	!! therefore means that it is reserved for developers and experienced
55	!! professionals having in-depth computer knowledge. Users are therefore
56	!! encouraged to load and test the software's suitability as regards their
57	!! requirements in conditions enabling the security of their systems and/or
58	!! data to be ensured and, more generally, to use and operate it in the
59	!! same conditions as regards security.
60	!!
61	!! The fact that you are presently reading this means that you have had
62	!! knowledge of the CeCILL license and that you accept its terms.
63	!! =========================================================================
64
65
66
67	SUBROUTINE CHEMMAIN( &
68	mmr , &
69	nstep , &
70	calday , &
71	ncdate , & !not used
72	ncsec , & !not used
73	lat , &
74	delt , &
75	ps , &
76	pmid , &
77	pdel , &
78	area , &
79	oro , &
80	tsurf , & ! not used
81	albs , &
82	zma , &
83	phis , &
84	cldfr , &
85	cldfr_st , &
86	cldfr_cv , &
87	cldtop , &
88	cldbot , &
89	cwat , &
90	flxrst , &
91	flxrcv , &
92	flxsst , &
93	flxscv , &
94	flxupd , &
95	flxmass_w , &
96	tfld , &
97	sh , &
98	ql , & ! variable en prevision de la chimie strato
99	rh , &
100	nx , &
101	ny , &
102	source)
103
104	!-----------------------------------------------------------------------
105	!
106	! INCA -- INteractions with Chemistry in the Atmosphere
107	!
108	! Advances the volumetric mixing ratio forward one time step via a
109	! combination of explicit, EBI, QSSA, fully implicit, and/or rodas
110	! algorithms.
111	!
112	! Didier Hauglustaine and Stacy Walters, 2000.
113	!-----------------------------------------------------------------------
114
115	USE CHEM_MODS, ONLY : nas, invariants, hrates, hrates_cv, o3_prod, o3_loss, o3_st_flx, &
116	so2_p_h2soh, so2_p_dmsoh, so2_p_dmsno3, so2_p_dmsooh, asmsam_p_dmsooh, &
117	dmso_p_dmsoh, asso4m_p_so2oh, wet3d_so2, wet3d_dms, wet3d_dmso, wet3d_hno3, wet3d_nh3, &
118	wet3d_noy, wet3d_h2o2, wet3d_hono, hno3_p_g, hno3_p_a, hno3_l_g, hno3_l_a, &
119	nh3_l_g, nh3_l_a, csno3_p_a1, csno3_p_a2, cino3_p_a, &
120	wet3d_asso4m, wet3d_asnh4m, wet3d_asno3m, wet3d_csno3m, wet3d_cino3m, &
121	hono_p_g, hono_p_a, hono_l_g, hono_l_a, &
122	asno3m_p_nh3hno3, asnh4m_p_nh3hno3, hno3_p_nh3hno3, nh3_p_nh3hno3, &
123	ASAPp1a_p, ASAPp2a_p, ASARp1a_p, ASARp2a_p, pom_p_g, &
124	co2_basprod, co2_nmhcprod, co2_radicalprod, &
125	hno3_prod, hno3_loss, co_prod, co_loss, ch4_loss, n2o_loss, nadv_mass, &
126	no_daytime, day_cnt , adv_mass, h2o, h2oc, extfrc, wetloss
127
128	USE SPECIES_NAMES
129	USE CHEM_TRACNM
130	USE CHEM_CONS, ONLY : gravit, uma
131	! USE MASS_DIAGS
132	USE AIRPLANE_SRC, ONLY : itrop
133	USE SFLX
134	USE TRANSPORT_CONTROLS
135	USE INCA_DIM
136	USE MOD_INCA_PARA
137
138	USE MOD_GRID_INCA, ONLY : PLON_GLO
139	USE RATE_INDEX_MOD
140
141	USE INCA_WRITE_FIELD_P
142	USE XIOS_INCA
143
144	USE OXYDANT_COM
145
146	USE PARAM_CHEM
147	USE LIGHTNING
148
149
150	IMPLICIT NONE
151
152	!-----------------------------------------------------------------------
153	! ... Dummy arguments
154	!-----------------------------------------------------------------------
155	INTEGER, INTENT(in) :: nstep ! time index
156	INTEGER, INTENT(in) :: lat ! latitude index (always 1 and place holder in LMDz)
157	INTEGER, INTENT(in) :: nx, ny ! dyn grid
158	INTEGER, INTENT(in) :: ncdate ! date at end present time step
159	INTEGER, INTENT(in) :: ncsec ! seconds relative to ncdate
160	INTEGER, INTENT(in) :: cldtop(PLON) ! cloud top level ( 1 ... PLEV )
161	INTEGER, INTENT(in) :: cldbot(PLON) ! cloud bot level ( 1 ... PLEV )
162	REAL, INTENT(in) :: calday ! time of year in days for midpoint time
163	REAL, INTENT(in) :: delt ! timestep in seconds
164	REAL, INTENT(inout) :: mmr(PLON,PLEV,8) ! xported species ( mmr )
165	REAL, INTENT(in) :: ps(PLON) ! surface press ( pascals )
166	REAL, INTENT(in) :: pmid(PLON,PLEV) ! midpoint press ( pascals )
167	REAL, INTENT(in) :: pdel(PLON,PLEV) ! delta press across midpoints
168	REAL, INTENT(in) :: area(PLON) ! grid cell area
169	REAL, INTENT(in) :: oro(PLON) ! surface orography flag
170	REAL, INTENT(in) :: tsurf(PLON) ! surface temperature
171	REAL, INTENT(in) :: zma(PLON,PLEV) ! abs geopot height at midpoints ( m )
172	REAL, INTENT(in) :: phis(PLON) ! surf geopot
173	REAL, INTENT(in) :: cldfr(PLON,PLEV) ! cloud fraction (all clouds)
174	REAL, INTENT(in) :: cldfr_cv(PLON,PLEV) ! cloud fraction (convective clouds)
175	REAL, INTENT(in) :: cldfr_st(PLON,PLEV) ! cloud fraction (large scale clouds)
176	REAL, INTENT(in) :: cwat(PLON,PLEV) ! total cloud water (kg/kg)
177	REAL, INTENT(in) :: flxrst(PLON,PLEVP) !liquid water flux (stratiform) kgH2O/m2/s
178	REAL, INTENT(in) :: flxrcv(PLON,PLEVP) !solid water flux (stratiform) kgH2O/m2/s
179	REAL, INTENT(in) :: flxsst(PLON,PLEVP) !liquid water flux (convection) kgH2O/m2/s
180	REAL, INTENT(in) :: flxscv(PLON,PLEVP) !solid water flux (convection) kgH2O/m2/s
181	REAL, INTENT(in) :: flxmass_w(PLON,PLEV) ! vertical mass flux
182	REAL, INTENT(in) :: flxupd(PLON,PLEV) ! entrainment flux kgAIR/m2/s
183	REAL, INTENT(in) :: tfld(PLON,PLEV) ! midpoint temperature
184	REAL, INTENT(in) :: sh(PLON,PLEV) ! specific humidity ( kg/kg )
185	REAL, INTENT(in) :: ql(PLON,PLEV) ! eau liquide
186	REAL, INTENT(in) :: rh(PLON,PLEV) ! relative humidity
187	REAL, INTENT(in) :: albs(PLON) ! surface albedo
188	REAL, INTENT(out) :: source(PLON,nbtrac)
189
190	!-----------------------------------------------------------------------
191	! ... Local variables
192	!-----------------------------------------------------------------------
193	INTEGER :: i, k, m, n, it,j
194	INTEGER :: klev
195	REAL :: col_dens(PLON,PLEV,2)
196	REAL, ALLOCATABLE, SAVE :: het_rates(:,:,:)
197	REAL, ALLOCATABLE, SAVE :: het_rates_st(:,:,:)
198	REAL, ALLOCATABLE, SAVE :: het_rates_cv(:,:,:)
199	!$OMP THREADPRIVATE(het_rates,het_rates_st,het_rates_cv)
200	REAL :: vmr(PLON,PLEV,8) ! xported species ( vmr )
201	REAL :: reaction_rates(PLON,PLEV,12)
202	REAL, DIMENSION(PLON,PLEV) :: h2ovmr ! water vapor volume mixing ratio
203	REAL, DIMENSION(PLON,PLEV) :: mbar ! mean wet atmospheric mass ( amu )
204	REAL, DIMENSION(PLON,PLEV) :: zmid ! midpoint geopotential in km
205	REAL, DIMENSION(8) :: vprodj ! volume production
206	REAL, DIMENSION(8) :: vlossj ! volume loss
207	REAL, DIMENSION(PLON) :: zen_angle
208	REAL, DIMENSION(PLON) :: loc_angle
209	REAL :: sunon(PLON)
210	REAL :: sunoff(PLON)
211	REAL :: delt_inverse
212	REAL :: zelev(PLON)
213	REAL :: zalt(PLON,PLEV)
214	LOGICAL :: polar_night(PLON)
215	LOGICAL :: polar_day(PLON)
216	LOGICAL :: zangtz(PLON)
217
218	REAL :: tfld_lim(PLON,PLEV)
219	REAL :: albs_lim(PLON)
220	REAL :: tfld_glo(PLON_GLO,PLEV)
221	REAL :: pmid_glo(PLON_GLO,PLEV)
222	REAL, DIMENSION(PLON,PLEV,12) :: reacflux
223	LOGICAL,SAVE :: first = .TRUE.
224	!$OMP THREADPRIVATE(first)
225
226
227	!-----------------------------------------------------------------------
228	! ... Function interface
229	!-----------------------------------------------------------------------
230	REAL :: TSECND
231
232	! Exemple utilisation de writefied pour le debuggage
233	! DO i=1,8
234	! CALL writefield_p('q_'//tracnam(i), mmr(:,:,i), PLEV)
235	! ENDDO
236
237	reaction_rates(:,:,:) = 0.0
238	IF (first) THEN
239	ALLOCATE(het_rates(PLON,PLEV,1))
240	ALLOCATE(het_rates_st(PLON,PLEV,1))
241	ALLOCATE(het_rates_cv(PLON,PLEV,1))
242	het_rates(:,:,:) = 0.0
243	het_rates_st(:,:,:) = 0.0
244	het_rates_cv(:,:,:) = hrates_cv(:,:,:)
245	first=.FALSE.
246
247	ENDIF
248
249	delt_inverse = 1. / delt
250
251	!-----------------------------------------------------------------------
252	! ... T and albedo adjustment for photorates
253	!-----------------------------------------------------------------------
254	tfld_lim(:,:) = MAX(tfld(:,:),226.)
255	albs_lim(:) = MIN(albs(:),0.5)
256
257	!-----------------------------------------------------------------------
258	! ... Calculate parameters for diurnal geometry
259	!-----------------------------------------------------------------------
260	CALL DIURNAL_GEOM( &
261	lat, calday, polar_night, polar_day, &
262	sunon, sunoff, loc_angle, zen_angle, &
263	zangtz)
264
265
266	!-----------------------------------------------------------------------
267	! ... Initialize xform between mass and volume mixing ratios
268	!-----------------------------------------------------------------------
269	CALL INTI_MR_XFORM( sh, mbar )
270
271	!-----------------------------------------------------------------------
272	! ... Xform from mmr to vmr
273	!-----------------------------------------------------------------------
274	CALL MMR2VMR( vmr, mmr, mbar )
275
276	!-----------------------------------------------------------------------
277	! ... Xform water vapor from mmr to vmr and adjust in stratosphere
278	!-----------------------------------------------------------------------
279	CALL ADJH2O( h2ovmr, sh,ql, h2oc, mbar, vmr )
280	!-----------------------------------------------------------------------
281	! ... Xform geopotential height from m to km
282	!-----------------------------------------------------------------------
283	zelev(:) = 1.e-3 * phis(:PLON) / gravit
284	zmid(:,:) = 1.e-3 * zma(:PLON,:) / gravit
285
286	DO i = 1, PLON
287	DO k = 1, PLEV
288	zalt(i,k) = zmid(i,k) + zelev(i)
289	END DO
290	END DO
291
292	!-----------------------------------------------------------------------
293	! ... Set the "invariants"
294	!-----------------------------------------------------------------------
295	CALL SETINV(invariants, tfld, h2ovmr, pmid)
296
297	!-----------------------------------------------------------------------
298	! ... Set the column densities at the upper boundary
299	!-----------------------------------------------------------------------
300	CALL SET_UB_COL(col_dens, vmr, invariants, pdel)
301
302
303	!-----------------------------------------------------------------------
304	! ... Set rates for "tabular" and user specified reactions
305	!-----------------------------------------------------------------------
306	CALL SETRXT( reaction_rates, tfld )
307
308	CALL USRRXT( &
309	reaction_rates(:,:,:) , &
310	tfld , &
311	invariants , &
312	rh , &
313	ps , &
314	pmid , &
315	mmr , &
316	invariants(1,1,1) , &
317	lat , &
318	zen_angle , &
319	vmr )
320
321	! changement d'unite des taux de reactions des reactions a 2 reactants
322	! ils sont calcules dans usrrxt en s-1/(mol.cm-3) et le code dans
323	! exp_sol et imp_sol attends des s-1
324	! ce n'est pas le cas des reactions a 1 reactant dont les taux
325	! sont directement calcule en s-1 par usrrxt (so4aerrxt etc.)
326	CALL ADJRXT(reaction_rates, invariants, invariants(1,1,1))
327
328	!-----------------------------------------------------------------------
329	! ... Compute the photolysis rates
330	!-----------------------------------------------------------------------
331	!-----------------------------------------------------------------------
332	! ... Set the column densities
333	!-----------------------------------------------------------------------
334	CALL SETCOL(col_dens, vmr, invariants, pdel)
335
336	!-----------------------------------------------------------------------
337	! ... Calculate the photodissociation rates
338	!-----------------------------------------------------------------------
339
340	CALL PHOTO(reaction_rates, pmid, pdel, &
341	tfld, zalt, col_dens, zen_angle, &
342	albs, cwat, cldfr, sunon, sunoff, &
343	polar_night, polar_day)
344
345	!-----------------------------------------------------------------------
346	! ... Adjust the photodissociation rates
347	!-----------------------------------------------------------------------
348	CALL PHTADJ(reaction_rates, invariants, invariants(1,1,1))
349
350
351	!-----------------------------------------------------------------------
352	! ... Compute the heterogeneous rates at time = t(n+1)
353	!-----------------------------------------------------------------------
354
355	! ... For stratiform precipitation
356
357	CALL SETHET( &
358	het_rates_st ,&
359	pmid ,&
360	pdel ,&
361	lat ,&
362	zmid ,&
363	tfld ,&
364	delt ,&
365	invariants(1,1,1) ,&
366	flxrst ,&
367	flxsst ,&
368	flxupd ,&
369	cldtop ,&
370	cldbot ,&
371	cldfr_st ,&
372	1 ,&
373	vmr )
374
375	! ... For convective precipitation
376
377	CALL SETHET( &
378	het_rates_cv ,&
379	pmid ,&
380	pdel ,&
381	lat ,&
382	zmid ,&
383	tfld ,&
384	delt ,&
385	invariants(1,1,1) ,&
386	flxrcv ,&
387	flxscv ,&
388	flxupd ,&
389	cldtop ,&
390	cldbot ,&
391	cldfr_cv ,&
392	2 ,&
393	vmr )
394
395	! ... Total removal rate
396
397	het_rates = het_rates_st + het_rates_cv
398	hrates = het_rates
399	hrates_cv = het_rates_cv
400
401	!-----------------------------------------------------------------------
402	! ... Compute the extraneous forcings
403	!-----------------------------------------------------------------------
404
405	!Initialize for this timeslice
406	DO m = 1,1
407	extfrc(:,:,m) = 0.
408	END DO
409
410	DO m = 1,1
411	extfrc(:,:,m) = 0.
412	END DO
413
414
415
416
417
418
419	!=======================================================================
420	! ... Call the class solution algorithms
421	!=======================================================================
422	!-----------------------------------------------------------------------
423	! ... Solve for "explicit" species
424	!-----------------------------------------------------------------------
425	CALL EXP_SOL(vmr, reaction_rates, het_rates, extfrc, &
426	nstep, delt)
427
428
429
430	DO i = 1, PLON
431	ENDDO
432
433	!-----------------------------------------------------------------------
434	! ... Check for negative values and reset to zero
435	!-----------------------------------------------------------------------
436	CALL NEGTRC( 'After chemistry ', vmr )
437
438
439	IF (calcul_flux) THEN
440	!-----------------------------------------------------------------------
441	! ... Compute the flux through each reaction
442	!-----------------------------------------------------------------------
443	CALL reac_flx(vmr, reaction_rates, invariants, reacflux)
444
445	CALL xios_inca_change_context("inca")
446	DO n=1, 12
447	CALL xios_inca_send_field("flux_"//trim(reacname(n)), reacflux(:,:,n))
448	ENDDO
449	CALL xios_inca_change_context("LMDZ")
450	ENDIF
451
452	!-----------------------------------------------------------------------
453	! ... Xform from vmr to mmr
454	!-----------------------------------------------------------------------
455	CALL VMR2MMR( vmr, mmr, &
456	mbar )
457
458
459
460	!DH Update the drydep flux with latest surface mixing ratio
461	CALL DRYDEP (pmid, &
462	tfld, &
463	mmr)
464
465	! calcul des sources
466	DO it=1,nbtrac
467	DO k = 1, PLON
468	source(k,it) = eflux(k,it)-dflux(k,it)
469	ENDDO
470	ENDDO
471
472
473	flux_source(:,:) = source(:,:)
474
475
476
477	END SUBROUTINE CHEMMAIN
478
479	SUBROUTINE endrun
480	!----------------------------------------------------------------------
481	! ... Abort the model
482	!----------------------------------------------------------------------
483
484	IMPLICIT NONE
485
486	EXTERNAL abort
487
488	CALL abort
489
490	END SUBROUTINE endrun
491

Note: See TracBrowser for help on using the repository browser.

Download in other formats: