Context Navigation

calchim.F90 @ 1971

Last change on this file since 1971 was 1958, checked in by jvatant, 8 years ago
Pseudo-evap was done in 2D in the photochem -> now 3D --JVO
File size: 20.5 KB

Line
1	SUBROUTINE calchim(ngrid,qy_c,declin,dtchim, &
2	ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
3
4	!---------------------------------------------------------------------------------------------------------
5	!
6	! Purpose : Interface subroutine to photochemical C model for Titan GCM.
7	! -------
8	! The subroutine computes the chemical processes for a single vertical column.
9	!
10	! - Only tendencies are returned.
11	! - With moyzon_ch=.true. and input vectors zonally averaged
12	! the calculation is done only once per lat. band
13	!
14	! Authors: + S. Lebonnois : 01/2000 \| 09/2003
15	! ------- adaptation for Titan 3D : 02/2009
16	! adaptation for // : 04/2013
17	! extension chemistry up to 1300km : 10/2013
18	!
19	! + J. Vatant d'Ollone
20	! + 02/17 - adaptation for the new generic-forked physics
21	! + 01/18 - 03/18 - Major transformations :
22	! - Upper chemistry fields are now stored in startfi
23	! and defined on a pressure grid from Vervack profile
24	! - These modifs enables to run chemistry with others resolution than 32x48x55 !
25	! - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
26	! - Chemistry can still be done in 2D
27	! -> Calcul. once per band lat and put same tendency in all longi.
28	! Check for negs in physiq_mod.
29	! -> If procs sharing a lat band, no problem, the calcul will just be done twice.
30	! -> Will not work with Dynamico, where the chemistry will have to be done in 3D.
31	! ( and there'll be work to do to get rid of averaged fields )
32	!
33	! + STILL TO DO : Replug the interaction with haze (cf titan.old) -> to see with JB.
34	!---------------------------------------------------------------------------------------------------------
35
36	! --------------------------------------------------------------------
37	! Structure :
38	! -----------
39	! 0. Declarations
40	! I. Init and firstcall
41	! 1. Read and store Vervack profile
42	! 2. Compute planetar averaged atm. properties
43	! 3. Init compound caracteristics
44	! 4. Init photodissociations rates from actinic fluxes
45	! 5. Init chemical reactions
46	! 6. Init eddy diffusion coeff
47	! II. Loop on latitudes/grid-points
48	! 0. Check on 2D chemistry
49	! 1. Compute atm. properties at grid points
50	! 2. Interpolate photodissociation rates at lat,alt,dec
51	! 3. Read composition
52	! 4. Call main solver gptitan C routine
53	! 5. Calculate output tendencies on advected tracers
54	! 6. Update upper chemistry fields
55	! 0bis. If 2D chemsitry, don't recalculate if needed
56	! -----------------------------------------------------------------
57
58	USE, INTRINSIC :: iso_c_binding
59	USE comchem_h
60	USE dimphy
61	USE datafile_mod, ONLY: datadir
62	USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
63	USE geometry_mod, ONLY: latitude
64	USE logic_mod, ONLY: moyzon_ch
65	USE moyzon_mod, ONLY: tmoy, playmoy
66
67	IMPLICIT NONE
68
69	! ------------------------------------------
70	! ********* 0. Declarations ***********
71	! ------------------------------------------
72
73	! Arguments
74	! ---------
75
76	INTEGER, INTENT(IN) :: ngrid ! Number of atmospheric columns.
77	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(IN) :: qy_c ! Chemical species on GCM layers after adv.+diss. (mol/mol).
78	REAL*8, INTENT(IN) :: declin ! Solar declination (rad).
79	REAL*8, INTENT(IN) :: dtchim ! Chemistry timsetep (s).
80	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: ctemp ! Mid-layer temperature (K).
81	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cpphi ! Mid-layer geopotential (m2.s-2).
82	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cplay ! Mid-layer pressure (Pa).
83	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: cplev ! Inter-layer pressure (Pa).
84	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: czlay ! Mid-layer altitude (m).
85	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: czlev ! Inter-layer altitude (m).
86
87	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT) :: dqyc ! Chemical species tendencies on GCM layers (mol/mol/s).
88
89	! Local variables :
90	! -----------------
91
92	INTEGER :: i , l, ic, ig, igm1
93
94	INTEGER :: dec, idec, ipres, ialt, klat
95
96	REAL*8 :: declin_c ! Declination (deg).
97	REAL*8 :: factp, factalt, factdec, factlat, krpddec, krpddecp1, krpddecm1
98	REAL*8 :: temp1, logp
99
100	! Variables sent into chemistry module (must be in double precision)
101	! ------------------------------------------------------------------
102
103	REAL*8, DIMENSION(nlaykim_tot) :: temp_c ! Temperature (K).
104	REAL*8, DIMENSION(nlaykim_tot) :: press_c ! Pressure (Pa).
105	REAL*8, DIMENSION(nlaykim_tot) :: phi_c ! Geopotential (m2.s-2) - actually not sent in chem. module but used to compute alts.
106	REAL*8, DIMENSION(nlaykim_tot) :: nb ! Density (cm-3).
107
108	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy ! Chemical species in whole column (mol/mol) sent to chem. module.
109	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy0 ! " " " " " " before modifs.
110
111	REAL*8 :: surfhaze(nlaykim_tot)
112	REAL*8 :: cprodaer(nlaykim_tot,4), cmaer(nlaykim_tot,4)
113	REAL*8 :: ccsn(nlaykim_tot,4), ccsh(nlaykim_tot,4)
114
115	REAL*8, DIMENSION(nlaykim_tot) :: rmil ! Mid-layer distance (km) to planetographic center.
116	REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
117	! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
118
119	! Saved variables initialized at firstcall
120	! ----------------------------------------
121
122	LOGICAL, SAVE :: firstcall = .TRUE.
123	!$OMP THREADPRIVATE(firstcall)
124
125	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: kedd ! Eddy mixing coefficient for upper chemistry (cm^2.s-1)
126	!$OMP THREADPRIVATE(kedd)
127
128	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: mass ! Molar mass of the compunds (g.mol-1)
129	REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: md
130	!$OMP THREADPRIVATE(mass,md)
131
132	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: r1d, ct1d, p1d, t1d ! Vervack profile
133	! JVO 18 : No threadprivate for those as they'll be read in tcp.ver by master
134
135	REAL*8, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: krpd ! Photodissociations rate table
136	REAL*8, DIMENSION(:,:) , ALLOCATABLE, SAVE :: krate ! Reactions rate ( photo + chem )
137	!$OMP THREADPRIVATE(krpd,krate)
138
139	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nom_prod, nom_perte
140	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: reactif
141	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: prod
142	INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: perte
143	!$OMP THREADPRIVATE(nom_prod,nom_perte,reactif,prod,perte)
144
145	! TEMPORARY : Dummy parameters without microphysics
146	! Here to keep the whole stuff running without modif chem. module
147	! ---------------------------------------------------------------
148
149	INTEGER :: utilaer(16)
150	INTEGER :: aerprod = 0
151	INTEGER :: htoh2 = 0
152
153	! -----------------------------------------------------------------------
154	! *************** I. Initialisations and Firstcall ******************
155	! -----------------------------------------------------------------------
156
157	IF (firstcall) THEN
158
159	PRINT*, 'CHIMIE, premier appel'
160
161	call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
162
163	ALLOCATE(r1d(131))
164	ALLOCATE(ct1d(131))
165	ALLOCATE(p1d(131))
166	ALLOCATE(t1d(131))
167
168	ALLOCATE(md(nlaykim_tot,nkim))
169	ALLOCATE(mass(nkim))
170
171	ALLOCATE(kedd(nlaykim_tot))
172
173	ALLOCATE(krate(nlaykim_tot,nr_kim))
174	ALLOCATE(krpd(nd_kim+1,nlrt_kim,15,nlat_actfluxes))
175
176	ALLOCATE(nom_prod(nkim))
177	ALLOCATE(nom_perte(nkim))
178	ALLOCATE(reactif(5,nr_kim))
179	ALLOCATE(prod(200,nkim))
180	ALLOCATE(perte(2,200,nkim))
181
182	! 1. Read Vervack profile "tcp.ver", once for all
183	! -----------------------------------------------
184
185	!$OMP MASTER
186	OPEN(11,FILE=TRIM(datadir)//'/tcp.ver',STATUS='old')
187	READ(11,*)
188	DO i=1,131
189	READ(11,*) r1d(i), t1d(i), ct1d(i), p1d(i)
190	! For debug :
191	! -----------
192	! PRINT*, "tcp.ver", r1d(i), t1d(i), ct1d(i), p1d(i)
193	ENDDO
194	CLOSE(11)
195	!$OMP END MASTER
196	!$OMP BARRIER
197
198	! 2. Calculation of temp_c, densities and altitudes in planetary average
199	! ----------------------------------------------------------------------
200
201	! JVO18 : altitudes are no more calculated in firstcall, as I set kedd in pressure grid
202
203	! a. For GCM layers we just copy-paste
204
205	PRINT*,'Init chemistry : pressure, density, temperature ... :'
206	PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'
207
208	IF (ngrid.NE.1) THEN
209	DO l=1,klev
210	temp_c(l) = tmoy(l) ! K
211	press_c(l) = playmoy(l)/100. ! mbar
212	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
213	PRINT*, l, press_c(l), nb(l), temp_c(l), phi_c(l)
214	ENDDO
215	ELSE
216	DO l=1,klev
217	temp_c(l) = ctemp(1,l) ! K
218	press_c(l) = cplay(1,l)/100. ! mbar
219	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
220	PRINT*, l, press_c(l), nb(l), temp_c(l)
221	ENDDO
222	ENDIF
223
224	! b. Extension in upper atmosphere with Vervack profile
225	! NB : Maybe the transition klev/klev+1 is harsh ...
226
227	ipres=1
228	DO l=klev+1,nlaykim_tot
229	press_c(l) = preskim(l-klev) / 100.0
230	DO i=ipres,130
231	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
232	ipres=i
233	ENDIF
234	ENDDO
235	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
236
237	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
238	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
239	PRINT*, l , press_c(l), nb(l), temp_c(l)
240	ENDDO
241
242	! 3. Compounds caracteristics
243	! ---------------------------
244	mass(:) = 0.0
245	call comp(nomqy_c,nb,temp_c,mass,md)
246	PRINT*,' Mass'
247	DO ic=1,nkim
248	PRINT*, nomqy_c(ic), mass(ic)
249	ENDDO
250
251	! 4. Photodissociation rates
252	! --------------------------
253	call disso(krpd,nlat_actfluxes)
254
255	! 5. Init. chemical reactions with planetary average T prof.
256	! ----------------------------------------------------------
257
258	! NB : Chemical reactions rate are assumed to be constant within the T range of Titan's atm
259	! so we fill their krate once for all but krate for photodiss will be filled at each timestep
260
261	call chimie(nomqy_c,nb,temp_c,krate,reactif, &
262	nom_perte,nom_prod,perte,prod)
263
264	! 6. Eddy mixing coefficients (constant with time and space)
265	! ----------------------------------------------------------
266
267	kedd(:) = 1.e3 ! Default value =/= zero
268
269	! NB : Eddy coeffs (e.g. Lavvas et al 08, Yelle et al 08) in altitude but they're rather linked to pressure
270	! Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
271
272	! First calculate kedd for upper chemistry layers
273	DO l=klev+1,nlaykim_tot
274	logp=-log10(press_c(l))
275	! 1E6 at 300 km ~ 10-1 mbar
276	IF ( logp.ge.1.0 .and. logp.le.4.0 ) THEN
277	kedd(l) = 10.*(6.0+1.3(logp-1.0)/3.0)
278	! 2E7 above 600 km ~ 10-4 mbar
279	ELSEIF ( logp.gt.4.0 ) THEN
280	kedd(l) = 2.e7
281	ENDIF
282	ENDDO
283
284	! Then adjust 15 last layers profile fading to default value depending on kedd(ptop)
285	DO l=klev-15,klev
286	temp1 = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev+1)/press_c(klev-15)) )
287	kedd(l) = 10.*( 3.0 + log10(kedd(klev+1)/1.e3) temp1 )
288	ENDDO
289
290	firstcall = .FALSE.
291	ENDIF ! firstcall
292
293	declin_c = declin*180./pi
294
295	! -----------------------------------------------------------------------
296	! ********************* II. Loop on latitudes ***********************
297	! -----------------------------------------------------------------------
298
299	DO ig=1,ngrid
300
301	IF (ig.eq.1) THEN
302	igm1=1
303	ELSE
304	igm1=ig-1
305	ENDIF
306
307	! If 2D chemistry, trick to do the calculation only once per latitude band within the chunk
308	! NB1 : Will be obsolete with DYNAMICO, the chemistry will necessarly be 3D
309	! NB2 : Test of same latitude with dlat=0.5 : I think that if you run sims better than 1/10th degree then
310	! either it's with Dynamico and doesn't apply OR it is more than enough in terms of "preco / calc time" !
311	! -------------------------------------------------------------------------------------------------------
312
313	IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.)) ) .OR. (.NOT. moyzon_ch) ) THEN
314
315	! 1. Compute altitude for the grid point with hydrostat. equilib.
316	! ---------------------------------------------------------------
317
318	! a. For GCM layers we just copy-paste
319
320	DO l=1,klev
321	rinter(l) = (czlev(ig,l)+rad)/1000.0 ! km
322	rmil(l) = (czlay(ig,l)+rad)/1000.0 ! km
323	temp_c(l) = ctemp(ig,l) ! K
324	phi_c(l) = cpphi(ig,l) ! m2.s-2
325	press_c(l) = cplay(ig,l)/100. ! mbar
326	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
327	ENDDO
328	rinter(klev+1)=(czlev(ig,klev+1)+rad)/1000.
329
330	! b. Extension in upper atmosphere with Vervack profile
331	! NB : The transition klev/klev+1 is maybe harsh if we have correct g above but not under
332
333	ipres=1
334	DO l=klev+1,nlaykim_tot
335	press_c(l) = preskim(l-klev) / 100.0
336	DO i=ipres,130
337	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
338	ipres=i
339	ENDIF
340	ENDDO
341	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
342
343	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
344	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
345	ENDDO
346
347	! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
348	! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
349
350	DO l=klev+1,nlaykim_tot
351
352	! Compute geopotential on the upper grid, to have correct altitude
353	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
354	! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
355	! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
356
357	temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
358	phi_c(l) = phi_c(l-1) + rtemp1log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
359
360	rmil(l) = ( gradrad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
361	ENDDO
362
363	DO l=klev+2,nlaykim_tot
364	rinter(l) = 0.5*(rmil(l-1) + rmil(l))
365	ENDDO
366
367	! 2. From krpd, compute krate for dissociations (declination-latitude-altitude interpolation)
368	! -------------------------------------------------------------------------------------------
369
370	! a. Calculate declination dependence
371
372	if ((declin_c*10+267).lt.14.) then
373	idec = 0
374	dec = 0
375	else
376	if ((declin_c*10+267).gt.520.) then
377	idec = 14
378	dec = 534
379	else
380	idec = 1
381	dec = 27
382	do while( (declin_c*10+267).ge.real(dec+20) )
383	dec = dec+40
384	idec = idec+1
385	enddo
386	endif
387	endif
388	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
389	factdec = ( declin_c - (dec-267)/10. ) / 4.
390	else
391	factdec = ( declin_c - (dec-267)/10. ) / 2.7
392	endif
393
394	! b. Calculate klat for interpolation on fixed latitudes of actinic fluxes input
395
396	klat=1
397	DO i=1,nlat_actfluxes
398	IF (latitude(ig).LT.lat_actfluxes(i)) klat=i
399	ENDDO
400	IF (klat==nlat_actfluxes) THEN ! avoid rounding problems
401	klat = nlat_actfluxes-1
402	factlat = 1.0
403	ELSE
404	factlat = (latitude(ig)-lat_actfluxes(klat))/(lat_actfluxes(klat+1)-lat_actfluxes(klat))
405	ENDIF
406
407	! c. Altitude loop
408
409	DO l=1,nlaykim_tot
410
411	! Calculate ialt for interpolation in altitude (krpd every 2 km)
412	ialt = int((rmil(l)-rad/1000.)/2.)+1
413	factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
414
415	! Altitude can go above top limit of UV levels - in this case we keep the 1310km top fluxes
416	IF (ialt.GT.nlrt_kim-1) THEN
417	ialt = nlrt_kim-1 ! avoid out-of-bound array
418	factalt = 1.0
419	ENDIF
420
421	DO i=1,nd_kim+1 ! nd_kim+1 is dissociation of N2 by GCR
422
423	krpddec = ( krpd(i,ialt ,idec+1,klat) * (1.0-factalt) &
424	+ krpd(i,ialt+1,idec+1,klat) * factalt ) * (1.0-factlat) &
425	+ ( krpd(i,ialt ,idec+1,klat+1) * (1.0-factalt) &
426	+ krpd(i,ialt+1,idec+1,klat+1) * factalt ) * factlat
427
428	if ( factdec.lt.0. ) then
429	krpddecm1 = ( krpd(i,ialt ,idec ,klat) * (1.0-factalt) &
430	+ krpd(i,ialt+1,idec ,klat) * factalt ) * (1.0-factlat) &
431	+ ( krpd(i,ialt ,idec ,klat+1) * (1.0-factalt) &
432	+ krpd(i,ialt+1,idec ,klat+1) * factalt ) * factlat
433	krate(l,i) = krpddecm1 * abs(factdec) + krpddec * ( 1.0 + factdec)
434	else if ( factdec.gt.0. ) then
435	krpddecp1 = ( krpd(i,ialt ,idec+2,klat) * (1.0-factalt) &
436	+ krpd(i,ialt+1,idec+2,klat) * factalt ) * (1.0-factlat) &
437	+ ( krpd(i,ialt ,idec+2,klat+1) * (1.0-factalt) &
438	+ krpd(i,ialt+1,idec+2,klat+1) * factalt ) * factlat
439	krate(l,i) = krpddecp1 * factdec + krpddec * ( 1.0 - factdec)
440	else if ( factdec.eq.0. ) then
441	krate(l,i) = krpddec
442	endif
443
444	ENDDO ! i=1,nd_kim+1
445	ENDDO ! l=1,nlaykim_tot
446
447	! 3. Read composition
448	! -------------------
449
450	DO ic=1,nkim
451	DO l=1,klev
452	cqy(l,ic) = qy_c(ig,l,ic) ! advected tracers for the GCM part converted to molar frac.
453	ENDDO
454
455	DO l=1,nlaykim_up
456	cqy(klev+l,ic) = ykim_up(ic,ig,l) ! ykim_up for the upper atm.
457	ENDDO
458	ENDDO
459
460	cqy0(:,:) = cqy(:,:) ! Stores compo. before modifs
461
462	! 4. Call main Titan chemistry C routine
463	! --------------------------------------
464
465	call gptitan(rinter,temp_c,nb, &
466	nomqy_c,cqy, &
467	dtchim,latitude(ig)*180./pi,mass,md, &
468	kedd,krate,reactif, &
469	nom_prod,nom_perte,prod,perte, &
470	aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, &
471	htoh2,surfhaze)
472
473	! 5. Calculates tendencies on composition for advected tracers
474	! ------------------------------------------------------------
475	DO ic=1,nkim
476	DO l=1,klev
477	dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
478	ENDDO
479	ENDDO
480
481	! 6. Update ykim_up
482	! -----------------
483	DO ic=1,nkim
484	DO l=1,nlaykim_up
485	ykim_up(ic,ig,l) = cqy(klev+l,ic)
486	ENDDO
487	ENDDO
488	! NB: The full vertical composition grid will be created only for the outputs
489
490
491	ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
492	dqyc(ig,:,:) = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
493	ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations
494	ENDIF
495
496	ENDDO
497
498	END SUBROUTINE calchim

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

Download in other formats:

Original Format