Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

calchim_asis.F90

Last change on this file was 3312, checked in by yjaziri, 7 months ago

Generic PCM:

Photochemistry: correction for using diurnal equal true in 1D
No correction factor with diurnal equal true for photolysis rate
+ writediagspecUV change name to flux_surf for clarity

Property svn:executable set to *

File size: 21.1 KB

Line
1	subroutine calchim_asis(ngrid,nlayer,nq, &
2	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,fract, &
3	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim, &
4	tauref,co2ice, &
5	pu,pdu,pv,pdv,surfdust,surfice,icount,zdtchim)
6
7	use tracer_h, only: mmol, noms, nesp, is_chim
8
9	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
10	USE comcstfi_mod
11	use callkeys_mod
12	use time_phylmdz_mod, only: ecritphy, iphysiq ! output rate set by ecritphy
13	use types_asis, only: v_phot_3d, jlabel, nb_phot_hv_max
14	use chimiedata_h
15	use radcommon_h, only: glat
16	use wstats_mod, only: wstats
17
18	implicit none
19
20	!=======================================================================
21	!
22	! subject:
23	! --------
24	!
25	! Prepare the call for the photochemical module, and send back the
26	! tendencies from photochemistry in the chemical species mass mixing ratios
27	!
28	! Author: Sebastien Lebonnois (08/11/2002)
29	! -------
30	! update 12/06/2003 for water ice clouds and compatibility with dust
31	! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
32	! update 03/05/2005 cosmetic changes (Franck Lefevre)
33	! update sept. 2008 identify tracers by their names (Ehouarn Millour)
34	! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
35	! update 16/03/2012 optimization (Franck Lefevre)
36	! update 11/12/2013 optimization (Franck Lefevre)
37	! update 20/10/2017 adapted to LMDZ GENERIC+cosmetic changes (Benjamin Charnay)
38	! update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
39	!
40	! Arguments:
41	! ----------
42	!
43	! Input:
44	!
45	! ptimestep timestep (s)
46	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
47	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
48	! pt(ngrid,nlayer) Temperature (K)
49	! pdt(ngrid,nlayer) Temperature tendency (K)
50	! pu(ngrid,nlayer) u component of the wind (ms-1)
51	! pdu(ngrid,nlayer) u component tendency (K)
52	! pv(ngrid,nlayer) v component of the wind (ms-1)
53	! pdv(ngrid,nlayer) v component tendency (K)
54	! dist_sol distance of the sun (AU)
55	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
56	! fract(ngrid) day fraction (for diurnal=false)
57	! pq(ngrid,nlayer,nq) Advected fields, ie chemical species here
58	! pdq(ngrid,nlayer,nq) Previous tendencies on pq
59	! tauref(ngrid) Optical depth at 7 hPa
60	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
61	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
62	! surfice(ngrid,nlayer) ice surface area (m2/m3)
63	!
64	! Output:
65	!
66	! dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry
67	! dqschim(ngrid,nesp) ! tendencies on qsurf
68	!
69	!=======================================================================
70
71	! input:
72
73	integer,intent(in) :: ngrid ! number of atmospheric columns
74	integer,intent(in) :: nlayer ! number of atmospheric layers
75	integer,intent(in) :: nq ! number of tracers
76	real :: ptimestep
77	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
78	real :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers
79	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
80	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
81	real :: pt(ngrid,nlayer) ! temperature
82	real :: pdt(ngrid,nlayer) ! temperature tendency
83	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
84	real :: pdu(ngrid,nlayer) ! u component tendency
85	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
86	real :: pdv(ngrid,nlayer) ! v component tendency
87	real :: dist_sol ! distance of the sun (AU)
88	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
89	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
90	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
91	real :: zday ! date (time since Ls=0, in martian days)
92	real :: tauref(ngrid) ! optical depth at 7 hPa
93	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
94	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
95	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
96
97	! output:
98
99	real :: dqchim(ngrid,nlayer,nesp) ! tendencies on pq due to chemistry
100	real :: dqschim(ngrid,nesp) ! tendencies on qsurf
101
102	! local variables:
103
104	integer :: iloc(1) ! index of major species
105	integer :: ig,l,i,iq,iqmax,iesp
106	integer :: foundswitch, lswitch ! to switch between photochem and thermochem ? (YJ)
107	integer,save :: chemthermod
108
109	real :: zq(ngrid,nlayer,nesp) ! pq+pdq*ptimestep before chemistry
110	! new mole fraction after
111	real :: zt(ngrid,nlayer) ! temperature
112	real :: zu(ngrid,nlayer) ! u component of the wind
113	real :: zv(ngrid,nlayer) ! v component of the wind
114	real :: taucol ! optical depth at 7 hPa
115	real :: xmmol(nlayer,nesp) ! mmol/mmean but only for chemical species
116
117	logical,save :: firstcall = .true.
118
119	! for each column of atmosphere:
120
121	real :: zpress(nlayer) ! Pressure (mbar)
122	real :: zdens(nlayer) ! Density (cm-3)
123	real :: ztemp(nlayer) ! Temperature (K)
124	real :: zlocal(nlayer) ! Altitude (km)
125	real :: zycol(nlayer,nesp) ! Composition (mole fractions)
126	real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1)
127	real :: szacol ! Solar zenith angle
128	real :: fract(ngrid) ! day fraction
129	real :: fractcol ! day fraction
130	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
131	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
132
133	integer :: iter(nlayer) ! Number of chemical iterations
134	! within one physical timestep
135	integer :: icount
136	! for output:
137
138	logical :: output ! to issue calls to writediagfi and stats
139	parameter (output = .true.)
140	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
141	! within one physical timestep
142	integer :: ierr
143	real :: zdtchim(ngrid,nlayer) ! temperature modification by chemistry
144	real :: dEzchim(ngrid,nlayer) ! energy modification by chemistry
145	real :: zdtchim_output(ngrid) ! flux modification by chemistry in W.m-2
146
147	!=======================================================================
148	! initialization of the chemistry (first call only)
149	!=======================================================================
150
151	if (firstcall) then
152
153	!! Moved to routine indice in photochemistry_asis
154	!! because nb_phot_hv_max value needed in order
155	!! to choose if we call read_phototable or not.
156	!! A cleaner solution need to be find.
157	! if (photochem .and. .not. jonline) then
158	! print*,'calchim: Read photolysis lookup table'
159	! call read_phototable
160	! end if
161
162	if (.not.allocated(SF_mode)) allocate(SF_mode(nesp))
163	if (.not.allocated(SF_value)) allocate(SF_value(nesp))
164	if (.not.allocated(prod_rate)) allocate(prod_rate(nesp))
165	if (.not.allocated(surface_flux)) allocate(surface_flux(ngrid,nesp))
166	if (.not.allocated(surface_flux2)) allocate(surface_flux2(ngrid,nesp))
167	if (.not.allocated(escape)) allocate(escape(ngrid,nesp))
168	if (.not.allocated(chemnoms)) allocate(chemnoms(nesp))
169
170	surface_flux(:,:) = 0.
171	surface_flux2(:,:) = 0.
172	escape(:,:) = 0.
173	SF_mode(:) = 2
174	SF_value(:) = 0.
175	prod_rate(:) = 0.
176	iter_3d(:,:) = 0.
177	iter(:) = 0.
178
179	call ini_tracchim
180
181	! Sanity check mmol /= 0. in chemistry
182	do iq = 1,nq
183	if (is_chim(iq).eq.1 .and. mmol(iq).eq.0.) then
184	write(,) 'Error in calchim:'
185	write(,) 'Mmol cannot be equal to 0 for chemical species'
186	stop
187	end if
188	end do
189
190	firstcall = .false.
191	end if ! if (firstcall)
192
193	! Initializations
194
195	zycol(:,:) = 0.
196	dqchim(:,:,:) = 0.
197	dqschim(:,:) = 0.
198
199	!=======================================================================
200	! loop over grid
201	!=======================================================================
202
203	do ig = 1,ngrid
204
205	foundswitch = 0
206	do l = 1,nlayer
207	iesp = 0
208	do iq = 1,nq
209	if (is_chim(iq).eq.1) then
210	iesp = iesp + 1
211	zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
212	xmmol(l,iesp) = mmol(iq)/mmean(ig,l)
213	zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp)
214	end if
215	end do
216
217	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
218	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
219	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
220	zpress(l) = pplay(ig,l)/100.
221	ztemp(l) = zt(ig,l)
222	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
223	zlocal(l) = zzlay(ig,l)/1000.
224	zmmean(l) = mmean(ig,l)
225
226	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
227
228	surfdust1d(l) = surfdust(ig,l)*1.e-2
229	surfice1d(l) = surfice(ig,l)*1.e-2
230
231	! search for switch index between regions
232
233	! if (photochem .and. thermochem) then
234	! if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
235	! lswitch = l
236	! foundswitch = 1
237	! end if
238	! end if
239	if (.not. photochem) then
240	lswitch = 22
241	end if
242	! if (.not. thermochem) then
243	lswitch = min(50,nlayer+1)
244	! end if
245
246	end do ! of do l=1,nlayer
247
248	szacol = acos(mu0(ig))*180./pi
249	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
250	fractcol = fract(ig)
251
252	!=======================================================================
253	! call chemical subroutines
254	!=======================================================================
255
256	! chemistry in lower atmosphere
257
258	call photochemistry_asis(nlayer,ngrid, &
259	ig,lswitch,zycol,szacol,fractcol,ptimestep,&
260	zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
261	surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:))
262
263	! diagnostic photochemical heating
264	zdtchim_output(ig) = 0.
265	do l = 1,nlayer
266	dEzchim(ig,l) = zdtchim(ig,l)cpp(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
267	zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)cpp(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
268	end do
269
270	do l = 1,nlayer
271	iter_3d(ig,l) = iter(l)
272	end do
273
274	! chemistry in upper atmosphere
275
276	! if (thermochem) then
277	! call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
278	! zdens,zpress,zlocal,szacol,ptimestep,zday)
279	! end if
280
281	! dry deposition
282
283	if (depos) then
284	call deposition_source(ngrid, nlayer, nq, &
285	ig, zzlay, zzlev, zdens, zycol, ptimestep)
286	surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:)
287	if (ngrid==1) then
288	if(mod(icount,ecritphy).eq.0) then
289	surface_flux2(ig,:) = surface_flux2(ig,:)/float(ecritphy)
290	endif
291	else
292	if(mod(icount*iphysiq,ecritphy).eq.0) then
293	surface_flux2(ig,:) = surface_flux2(ig,:)*iphysiq/float(ecritphy)
294	endif
295	endif
296	end if
297
298	!=======================================================================
299	! tendencies
300	!=======================================================================
301
302	! index of the most abundant species at each level
303
304	! major(:) = maxloc(zycol, dim = 2)
305
306	! tendency for the most abundant species = - sum of others
307
308	do l = 1,nlayer
309	iloc=maxloc(zycol(l,:))
310	iqmax=iloc(1)
311	do iq = 1,nesp
312	if (iq /= iqmax) then
313	dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep
314	dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep)
315	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq)
316	end if
317	end do
318
319	end do ! of do l = 1,nlayer
320
321
322
323
324
325
326	!=======================================================================
327	! end of loop over grid
328	!=======================================================================
329
330	end do ! of do ig=1,ngridbidon(ig,:) = 1
331
332	!=======================================================================
333	! write outputs
334	!=======================================================================
335
336	! value of parameter 'output' to trigger writting of outputs
337	! is set above at the declaration of the variable.
338	if (output) then
339
340	! photoloysis
341	do i=1,nb_phot_hv_max
342	call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1), &
343	's-1',3,v_phot_3d(1,1,i))
344	end do
345	! iter
346	call wstats(ngrid,'iter','iterations', &
347	' ',3,iter_3d(1,1))
348	! phothcemical heating
349	call wstats(ngrid,'zdtchim','dT chim', &
350	'K.s-1',3,zdtchim)
351	call wstats(ngrid,'dEzchim','dE chim', &
352	'W m-2',3,dEzchim)
353	call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
354	! Mean molecular mass
355	call wstats(ngrid,'mmean','mean molecular mass', &
356	'g.mole-1',3,mmean(1,1))
357	! Chemical tendencies
358	do iesp=1,nesp
359	call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', &
360	'kg/kg s^-1',3,dqchim(1,1,iesp) )
361	end do
362	if (depos) then
363	! Surface fluxes and escape
364	do iesp=1,nesp
365	call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', &
366	'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
367	call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', &
368	'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
369	call wstats(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape', &
370	'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
371	end do
372	endif ! end depos
373
374	! photoloysis
375	do i=1,nb_phot_hv_max
376	call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1), &
377	's-1',3,v_phot_3d(1,1,i))
378	end do
379	! iter
380	call writediagfi(ngrid,'iter','iterations', &
381	' ',3,iter_3d(1,1))
382	! phothcemical heating
383	call writediagfi(ngrid,'zdtchim','dT chim', &
384	'K.s-1',3,zdtchim)
385	call writediagfi(ngrid,'dEzchim','dE chim', &
386	'W m-2',3,dEzchim)
387	call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
388	! Mean molecular mass
389	call writediagfi(ngrid,'mmean','mean molecular mass', &
390	'g.mole-1',3,mmean(1,1))
391	! Chemical tendencies
392	do iesp=1,nesp
393	call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq', &
394	'kg/kg s^-1',3,dqchim(1,1,iesp) )
395	end do
396	if (depos) then
397	! Surface fluxes and escape
398	do iesp=1,nesp
399	call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux', &
400	'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
401	call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux', &
402	'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
403	call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape', &
404	'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
405	end do
406	! Restart flux average
407	if (ngrid == 1) then
408	if(mod(icount,ecritphy).eq.0) then
409	surface_flux2(:,:) = 0.0
410	endif
411	else
412	if(mod(icount*iphysiq,ecritphy).eq.0) then
413	surface_flux2(:,:) = 0.0
414	endif
415	endif
416	endif ! end depos
417
418	if (jonline .and. output_writediagspecUV) then
419
420	!flux at the surface
421	call writediagspecUV(ngrid,"flux_surf","Surface flux(lon,lat,band)","photon.s-1.nm-1.cm-2",3,fluxUV(:,:,1))
422
423	endif ! end jonline
424
425	end if ! of if (output)
426
427	return
428
429
430	contains
431
432
433	!======================================================================
434
435	subroutine ini_tracchim
436
437	!======================================================================
438	! YJ Modern tracer.def
439	! Get in tracer.def chemical parameters
440	!======================================================================
441
442	use chimiedata_h
443	use tracer_h, only: is_chim
444	implicit none
445
446	! local
447	character(len=500) :: tracline ! to store a line of text
448	integer :: nq ! number of tracers
449	integer :: iesp, iq
450
451	! Get nq
452	open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
453	if (ierr.eq.0) then
454	READ(90,'(A)') tracline
455	IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
456	READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def
457	ELSE
458	DO
459	READ(90,'(A)',iostat=ierr) tracline
460	IF (ierr.eq.0) THEN
461	IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
462	READ(tracline,*,iostat=ierr) nq
463	EXIT
464	ENDIF
465	ELSE ! If pb, or if reached EOF without having found number of tracer
466	write(,) "calchim: error reading line of tracers"
467	write(,) " (first lines of traceur.def) "
468	stop
469	ENDIF
470	ENDDO
471	ENDIF ! if modern or standard traceur.def
472	else
473	write(,) "calchim: error opening traceur.def"
474	stop
475	endif
476
477	! Get data of tracers
478	iesp = 0
479	do iq=1,nq
480	read(90,'(A)') tracline
481	if (is_chim(iq).eq.1) then
482	iesp = iesp + 1
483	read(tracline,*) chemnoms(iesp)
484	write(,)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp))
485	if (index(tracline,'SF_mode=' ) /= 0) then
486	read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp)
487	write(,) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp)
488	else
489	write(,) ' Parameter value (default) : SF_mode=', SF_mode(iesp)
490	end if
491	if (index(tracline,'SF_value=' ) /= 0) then
492	read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp)
493	write(,) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp)
494	else
495	write(,) ' Parameter value (default) : SF_value=', SF_value(iesp)
496	end if
497	if (index(tracline,'prod_rate=' ) /= 0) then
498	read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp)
499	write(,) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp)
500	else
501	write(,) ' Parameter value (default) : prod_rate=', prod_rate(iesp)
502	end if
503	if (index(tracline,'escape=' ) /= 0) then
504	read(tracline(index(tracline,'escape=')+len('escape='):),*) escape(:,iesp)
505	write(,) ' Parameter value (traceur.def) : escape=', escape(:,iesp)
506	else
507	write(,) ' Parameter value (default) : escape=', escape(:,iesp)
508	end if
509	end if
510	end do
511
512	close(90)
513
514	end subroutine ini_tracchim
515
516	end subroutine calchim_asis
517

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

Download in other formats:

Original Format