Context Navigation

vdifc_pluto_mod.F90 @ 3853

Last change on this file since 3853 was 3853, checked in by tbertrand, 5 months ago
Pluto PCM: Implementing the CO cycle in the mode no_n2frost in the VTM and GCM TB
File size: 24.9 KB

Line
1	module vdifc_pluto_mod
2
3	implicit none
4
5	contains
6
7	SUBROUTINE vdifc_pluto(ngrid,nlay,nq,ppopsk, &
8	ptimestep,pcapcal,lecrit, &
9	pplay,pplev,pzlay,pzlev,pz0, &
10	pu,pv,ph,pq,pt,ptsrf,pemis,pqsurf, &
11	pdufi,pdvfi,pdhfi,pdqfi,pdtfi,pfluxsrf, &
12	pdudif,pdvdif,pdhdif,pdtsrf,pq2, &
13	pdqdif,pdqsdif,qsat_ch4,qsat_ch4_l1)
14
15	use comgeomfi_h
16	use callkeys_mod, only: carbox, methane, condcosurf, condensn2, condmetsurf,&
17	kmix_proffix, vertdiff, tracer, kmixmin, no_n2frost
18	use datafile_mod, only: datadir
19	use surfdat_h, only: phisfi
20	use comcstfi_mod, only: g, r, rcp, cpp
21	USE tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, igcm_co_gas, igcm_co_ice,&
22	igcm_n2, lw_ch4, lw_co, lw_n2
23
24	implicit none
25
26	!=======================================================================
27	!
28	! subject:
29	! --------
30	! Turbulent diffusion (mixing) for potential T, U, V and tracer
31	!
32	! Shema implicite
33	! On commence par rajouter au variables x la tendance physique
34	! et on resoult en fait:
35	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
36	!
37	! author:
38	! ------
39	! Hourdin/Forget/Fournier
40	!=======================================================================
41
42	!-----------------------------------------------------------------------
43	! declarations:
44	! -------------
45
46	#include "dimensions.h"
47
48	!
49	! arguments:
50	! ----------
51
52	INTEGER ngrid,nlay
53	REAL ptimestep
54	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
55	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
56	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
57	REAL ptsrf(ngrid),pemis(ngrid)
58	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
59	REAL pdtfi(ngrid,nlay)
60	REAL pt(ngrid,nlay)
61	REAL pfluxsrf(ngrid)
62	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
63	REAL pdtsrf(ngrid),pcapcal(ngrid)
64	REAL pq2(ngrid,nlay+1)
65	REAL qsat_ch4(ngrid)
66	REAL qsat_co(ngrid)
67	REAL qsat_ch4_l1(ngrid)
68	REAL zq1temp_ch4(ngrid)
69
70	! Argument added for condensation:
71	! REAL n2ice (ngrid)
72	REAL ppopsk(ngrid,nlay)
73	logical lecrit
74	REAL pz0
75
76	! Traceurs :
77	integer nq
78	REAL pqsurf(ngrid,nq)
79	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
80	real pdqdif(ngrid,nlay,nq)
81	real pdqdifeddy(ngrid,nlay,nq)
82	real pdqsdif(ngrid,nq),pdqsdif1(ngrid,nq)
83
84	! local:
85	! ------
86
87	INTEGER ilev,ig,ilay,nlev
88
89	REAL z4st,zdplanck(ngrid)
90	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
91	REAL zcdv(ngrid),zcdh(ngrid),sat2(ngrid)
92	REAL zcdv_true(ngrid),zcdh_true(ngrid)
93	REAL zu(ngrid,nlay),zv(ngrid,nlay)
94	REAL zh(ngrid,nlay),zt(ngrid,nlay)
95	REAL ztsrf2(ngrid),sat(ngrid),sat1(ngrid)
96	REAL z1(ngrid),z2(ngrid)
97	REAL za(ngrid,nlay),zb(ngrid,nlay)
98	REAL zb0(ngrid,nlay)
99	REAL zc(ngrid,nlay),zd(ngrid,nlay)
100	REAL zcst1
101	REAL zu2
102	EXTERNAL SSUM,SCOPY
103	REAL SSUM
104	LOGICAL firstcall
105	SAVE firstcall
106	real,dimension(:),save,allocatable :: qsat_co_factor ! factor to prevent co frost formation if no n2 frost
107	!$OMP THREADPRIVATE(qsat_co_factor)
108
109	!!read fixed profile for kmix
110	integer Nfine,ifine
111	parameter(Nfine=701)
112	character(len=100) :: file_path
113	real,save :: levdat(Nfine),kmixdat(Nfine)
114	real :: kmix_z(nlay) ! kmix from kmix_proffix
115
116	! variable added for N2 condensation:
117	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118	REAL hh , zhcond(ngrid,nlay)
119	! REAL latcond,tcond1p4Pa
120	REAL tcond1p4Pa
121	REAL acond,bcond
122	SAVE acond,bcond
123	! DATA latcond,tcond1p4Pa/2.5e5,38/
124	DATA tcond1p4Pa/38/
125
126	! Tracers :
127	! ~~~~~~~
128	INTEGER iq
129	REAL zq(ngrid,nlay,nq)
130	REAL zq1temp_co(ngrid)
131	REAL rho(ngrid) ! near surface air density
132	DATA firstcall/.true./
133
134	! ** un petit test de coherence
135	! --------------------------
136
137	IF (firstcall) THEN
138	IF(ngrid.NE.ngrid) THEN
139	write(,) 'STOP dans vdifc'
140	write(,) 'probleme de dimensions :'
141	write(,) 'ngrid =',ngrid
142	write(,) 'ngrid =',ngrid
143	STOP
144	ENDIF
145	! To compute: Tcond= 1./(bcond-acondlog(.7143p)) (p in pascal)
146	write(,) 'In vdifc: Tcond(P=1.4Pa)=',tcond1p4Pa,' Lcond=',lw_n2
147	bcond=1./tcond1p4Pa
148	acond=r/lw_n2
149	write(,) ' acond,bcond',acond,bcond
150
151	firstcall=.false.
152
153	! If fixed profile of kmix
154	IF (kmix_proffix) then
155	file_path=trim(datadir)//'/gas_prop/kmix.txt'
156	open(114,file=file_path,form='formatted')
157	do ifine=1,Nfine
158	read(114,*) levdat(ifine), kmixdat(ifine)
159	enddo
160	close(114)
161	ENDIF
162
163	! If fixed distribution of N2, then no CO frost either
164	ALLOCATE(qsat_co_factor(ngrid))
165	qsat_co_factor(:)=1.
166	IF (no_n2frost) then
167	DO ig=1,ngrid
168	if (pqsurf(ig,igcm_n2).eq.0.) then
169	qsat_co_factor(ig) = 1.e6
170	endif
171	ENDDO
172	ENDIF
173
174	ENDIF
175
176	!-----------------------------------------------------------------------
177	! 1. initialisation
178	! -----------------
179
180	nlev=nlay+1
181
182	! ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
183	! avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
184	! ----------------------------------------
185
186	DO ilay=1,nlay
187	DO ig=1,ngrid
188	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
189	ENDDO
190	ENDDO
191
192	zcst1=4.gptimestep/(r*r)
193	DO ilev=2,nlev-1
194	DO ig=1,ngrid
195	zb0(ig,ilev)=pplev(ig,ilev)* &
196	(pplev(ig,1)/pplev(ig,ilev))**rcp / &
197	(ph(ig,ilev-1)+ph(ig,ilev))
198	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/ &
199	(pplay(ig,ilev-1)-pplay(ig,ilev))
200	! write(300,*)'zb0',zb0(ig,ilev)
201	ENDDO
202	ENDDO
203	DO ig=1,ngrid
204	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
205	ENDDO
206
207	! ** diagnostique pour l'initialisation
208	! ----------------------------------
209
210	! IF(lecrit) THEN
211	! ig=ngrid/2+1
212	! write(,) 'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
213	! DO ilay=1,nlay
214	! WRITE(*,'(6f11.5)') &
215	! .01pplay(ig,ilay),.001pzlay(ig,ilay), &
216	! pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
217	! ENDDO
218	! write(,) 'Pression (mbar) ,altitude (km),zb'
219	! DO ilev=1,nlay
220	! WRITE(*,'(3f15.7)') &
221	! .01pplev(ig,ilev),.001pzlev(ig,ilev), &
222	! zb0(ig,ilev)
223	! ENDDO
224	! ENDIF
225
226	! Potential Condensation temperature:
227	! -----------------------------------
228
229	if (condensn2) then
230	DO ilev=1,nlay
231	DO ig=1,ngrid
232	zhcond(ig,ilev) = &
233	(1./(bcond-acondlog(.7143pplay(ig,ilev))))/ppopsk(ig,ilev)
234	END DO
235	END DO
236	else
237	DO ilev=1,nlay
238	DO ig=1,ngrid
239	zhcond(ig,ilev) = 0.
240	END DO
241	END DO
242	end if
243
244
245	!-----------------------------------------------------------------------
246	! 2. ajout des tendances physiques
247	! -----------------------------
248
249	DO ilev=1,nlay
250	DO ig=1,ngrid
251	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
252	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
253	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
254	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
255	zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
256	ENDDO
257	ENDDO
258	if(tracer) then
259	DO iq =1, nq
260	DO ilev=1,nlay
261	DO ig=1,ngrid
262	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
263	ENDDO
264	ENDDO
265	ENDDO
266	end if
267
268	!-----------------------------------------------------------------------
269	! 3. schema de turbulence
270	! --------------------
271
272	! ** source d'energie cinetique turbulente a la surface
273	! (condition aux limites du schema de diffusion turbulente
274	! dans la couche limite
275	! ---------------------
276	CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph &
277	,zcdv_true,zcdh_true)
278	DO ig=1,ngrid
279	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
280	!TB16: GCM wind for flat hemisphere
281	IF (phisfi(ig).eq.0.) zu2=max(2.,zu2)
282
283	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
284	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
285	ENDDO
286
287	! ** schema de diffusion turbulente dans la couche limite
288	! ----------------------------------------------------
289
290	CALL vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay &
291	,pu,pv,ph,zcdv_true &
292	,pq2,zkv,zkh)
293
294
295	! Adding eddy mixing to mimic 3D general circulation in 1D
296	! RW FF 2010
297	if ((ngrid.eq.1)) then
298	!kmixmin is the minimum eddy mix coeff in 1D
299
300	! If fixed profile of kmix
301	IF (kmix_proffix) then
302	!! Interpolate on the model vertical grid
303	CALL interp_line(levdat,kmixdat,Nfine,&
304	pzlay(1,:)/1000.,kmix_z(:),nlay)
305
306	do ilev=1,nlay
307	zkh(1,ilev) = max(kmix_z(ilev),zkh(1,ilev))
308	zkv(1,ilev) = max(kmix_z(ilev),zkv(1,ilev))
309	!zkh(1,ilev) = kmixmin
310	!zkv(1,ilev) = kmixmin
311	end do
312	ELSE
313	do ilev=1,nlay
314	zkh(1,ilev) = max(kmixmin,zkh(1,ilev))
315	zkv(1,ilev) = max(kmixmin,zkv(1,ilev))
316	!zkh(1,ilev) = kmixmin
317	!zkv(1,ilev) = kmixmin
318	end do
319	ENDIF
320	endif ! ngrid.eq.1
321
322	!! Temporary:
323	! zkh = zkh*0.1
324	! zkv = zkv*0.1
325
326	! ** diagnostique pour le schema de turbulence
327	! -----------------------------------------
328
329	! IF(lecrit) THEN
330	! write(,) ! write(,) 'Diagnostic for the vertical turbulent mixing'
331	! write(,) 'Cd for momentum and potential temperature'
332
333	! write(,) zcdv(ngrid/2+1),zcdh(ngrid/2+1)
334	! write(,) 'Mixing coefficient for momentum and pot.temp.'
335	! DO ilev=1,nlay
336	! write(,) zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
337	! ENDDO
338	! ENDIF
339
340	!-----------------------------------------------------------------------
341	! 4. inversion pour l'implicite sur u
342	! --------------------------------
343
344	! ** l'equation est
345	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
346	! avec
347	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
348	! et
349	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
350	! donc les entrees sont /zcdv/ pour la condition a la limite sol
351	! et /zkv/ = Ku
352
353	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
354	CALL multipl(ngrid,zcdv,zb0,zb)
355
356	DO ig=1,ngrid
357	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
358	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
359	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
360	ENDDO
361
362	DO ilay=nlay-1,1,-1
363	DO ig=1,ngrid
364	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
365	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
366	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ &
367	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
368	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
369	ENDDO
370	ENDDO
371
372	DO ig=1,ngrid
373	zu(ig,1)=zc(ig,1)
374	ENDDO
375	DO ilay=2,nlay
376	DO ig=1,ngrid
377	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
378	ENDDO
379	ENDDO
380
381	!-----------------------------------------------------------------------
382	! 5. inversion pour l'implicite sur v
383	! --------------------------------
384
385	! ** l'equation est
386	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
387	! avec
388	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
389	! et
390	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
391	! donc les entrees sont /zcdv/ pour la condition a la limite sol
392	! et /zkv/ = Kv
393
394	DO ig=1,ngrid
395	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
396	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
397	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
398	ENDDO
399
400	DO ilay=nlay-1,1,-1
401	DO ig=1,ngrid
402	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
403	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
404	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ &
405	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
406	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
407	ENDDO
408	ENDDO
409
410	DO ig=1,ngrid
411	zv(ig,1)=zc(ig,1)
412	ENDDO
413	DO ilay=2,nlay
414	DO ig=1,ngrid
415	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
416	ENDDO
417	ENDDO
418
419	!-----------------------------------------------------------------------
420	! 6. inversion pour l'implicite sur h sans oublier le couplage
421	! avec le sol (conduction)
422	! ------------------------
423
424	! ** l'equation est
425	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
426	! avec
427	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
428	! et
429	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
430	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
431	! et /zkh/ = Kh
432	! -------------
433
434	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
435	CALL multipl(ngrid,zcdh,zb0,zb)
436
437	DO ig=1,ngrid
438	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
439	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
440	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
441	ENDDO
442
443	DO ilay=nlay-1,1,-1
444	DO ig=1,ngrid
445	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
446	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
447	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ &
448	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
449	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
450	ENDDO
451	ENDDO
452
453	! ** calcul de (d Planck / dT) a la temperature d'interface
454	! ------------------------------------------------------
455
456	z4st=4.5.67e-8ptimestep
457	DO ig=1,ngrid
458	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
459	ENDDO
460
461	! ** calcul de la temperature_d'interface et de sa tendance.
462	! on ecrit que la somme des flux est nulle a l'interface
463	! a t + \delta t,
464	! !'est a dire le flux radiatif a {t + \delta t}
465	! + le flux turbulent a {t + \delta t} &
466	! qui 'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
467	! (notation K dt = /cpp*b/)
468	! + le flux dans le sol a t
469	! + l'evolution du flux dans le sol lorsque la temperature d'interface
470	! passe de sa valeur a t a sa valeur a {t + \delta t}.
471	! ----------------------------------------------------
472
473	DO ig=1,ngrid
474	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1) &
475	+zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
476	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
477	ztsrf2(ig)=z1(ig)/z2(ig)
478	pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
479
480	! Modif speciale N2 condensation:
481	! tconds = 1./(bcond-acondlog(.7143pplev(ig,1)))
482	! if ((condensn2).and. &
483	! ((n2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
484	! zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
485	! else
486	zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
487	! end if
488	ENDDO
489
490	! ** et a partir de la temperature au sol on remonte
491	! -----------------------------------------------
492
493	DO ilay=2,nlay
494	DO ig=1,ngrid
495	hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif n2cond
496	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
497	ENDDO
498	ENDDO
499
500
501	!-----------------------------------------------------------------------
502	! TRACERS
503	! -------
504
505	if(tracer) then
506
507	! Using the wind modified by friction for lifting and sublimation
508	! ----------------------------------------------------------------
509	! This is computed above
510
511	! DO ig=1,ngrid
512	! zu2=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
513	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
514	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
515	! ENDDO
516
517	! Calcul flux vertical au bas de la premiere couche (cf dust on Mars)
518	! -----------------------------------------------------------
519	do ig=1,ngrid
520	rho(ig) = zb0(ig,1) /ptimestep
521	! zb(ig,1) = 0.
522	end do
523
524	pdqsdif(:,:) = 0
525	pdqdif(:,:,:)=0.
526
527
528	! TB: Eddy lifting of tracers :
529	! ****************************************************************
530	! DO ig=1,ngrid
531	!! option : injection only on an equatorial band
532	!! if (abs(lati(ig))*180./pi.le.25.) then
533	! pdqsdif(ig,igcm_eddy1e6) =-1.e-12
534	! pdqsdif(ig,igcm_eddy1e7) =-1.e-12
535	! pdqsdif(ig,igcm_eddy5e7) =-1.e-12
536	! pdqsdif(ig,igcm_eddy1e8) =-1.e-12
537	! pdqsdif(ig,igcm_eddy5e8) =-1.e-12
538	! endif
539	! ENDDO
540
541	! pdqdifeddy(:,:,:)=0.
542	! injection a 50 km
543	! DO ig=1,ngrid
544	! pdqdifeddy(ig,17,igcm_eddy1e6)=1e-12
545	! pdqdifeddy(ig,17,igcm_eddy1e7)=1e-12
546	! pdqdifeddy(ig,17,igcm_eddy5e7)=1e-12
547	! pdqdifeddy(ig,17,igcm_eddy1e8)=1e-12
548	! pdqdifeddy(ig,17,igcm_eddy5e8)=1e-12
549	! ENDDO
550
551	! Inversion pour l'implicite sur q
552	! --------------------------------
553	do iq=1,nq
554	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
555
556	if ((methane).and.(iq.eq.igcm_ch4_gas)) then
557	! This line is required to account for turbulent transport
558	! from surface (e.g. ice) to mid-layer of atmosphere:
559	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
560	else if ((carbox).and.(iq.eq.igcm_co_gas)) then
561	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
562	else ! (re)-initialize zb(:,1)
563	zb(1:ngrid,1)=0
564	end if
565
566	DO ig=1,ngrid
567	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
568	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
569	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
570	ENDDO
571
572	DO ilay=nlay-1,2,-1
573	DO ig=1,ngrid
574	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
575	zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
576	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ &
577	zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
578	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
579	ENDDO
580	ENDDO
581
582	! special case for methane and CO ice tracer: do not include
583	! ice tracer from surface (which is set when handling
584	! vapour case (see further down).
585	if (methane.and.(iq.eq.igcm_ch4_ice)) then
586	DO ig=1,ngrid
587	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
588	zb(ig,2)*(1.-zd(ig,2)))
589	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
590	zb(ig,2)zc(ig,2)) z1(ig)
591	ENDDO
592	else if (carbox.and.(iq.eq.igcm_co_ice)) then
593	DO ig=1,ngrid
594	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
595	zb(ig,2)*(1.-zd(ig,2)))
596	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
597	zb(ig,2)zc(ig,2)) z1(ig)
598	ENDDO
599
600	else ! general case
601	DO ig=1,ngrid
602	z1(ig)=1./(za(ig,1)+zb(ig,1)+ &
603	zb(ig,2)*(1.-zd(ig,2)))
604	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
605	zb(ig,2)*zc(ig,2) + &
606	(-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
607	ENDDO
608	endif ! of if (methane.and.(iq.eq.igcm_ch4_ice))
609
610	! Calculation for turbulent exchange with the surface (for ice)
611	IF (methane.and.(iq.eq.igcm_ch4_gas)) then
612
613	!! calcul de la valeur de q a la surface :
614	call methanesat(ngrid,ptsrf,pplev(1,1), &
615	qsat_ch4(:),pqsurf(:,igcm_n2))
616
617	!! For output:
618	call methanesat(ngrid,zt(:,1),pplev(1,1), &
619	qsat_ch4_l1(:),pqsurf(:,igcm_n2))
620
621	!! Prevent CH4 condensation at the surface
622	if (.not.condmetsurf) then
623	qsat_ch4=qsat_ch4*1.e6
624	endif
625
626	DO ig=1,ngrid
627	zd(ig,1)=zb(ig,1)*z1(ig)
628	zq1temp_ch4(ig)=zc(ig,1)+ zd(ig,1)*qsat_ch4(ig)
629	pdqsdif(ig,igcm_ch4_ice)=rho(ig)*zcdv(ig) &
630	*(zq1temp_ch4(ig)-qsat_ch4(ig))
631	END DO
632
633	DO ig=1,ngrid
634	if ((-pdqsdif(ig,igcm_ch4_ice)*ptimestep) &
635	.gt.(pqsurf(ig,igcm_ch4_ice))) then
636
637	pdqsdif(ig,igcm_ch4_ice)= &
638	-pqsurf(ig,igcm_ch4_ice)/ptimestep
639
640	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
641
642	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_ch4_gas)+ &
643	zb(ig,2)*zc(ig,2) + &
644	(-pdqsdif(ig,igcm_ch4_ice)) ptimestep) z1(ig)
645
646	zq1temp_ch4(ig)=zc(ig,1)
647	endif
648
649	zq(ig,1,igcm_ch4_gas)=zq1temp_ch4(ig)
650
651	! TB: MODIF speciale pour CH4
652	pdtsrf(ig)=pdtsrf(ig)+ &
653	lw_ch4*pdqsdif(ig,igcm_ch4_ice)/pcapcal(ig)
654
655
656	ENDDO ! of DO ig=1,ngrid
657
658	ELSE IF (carbox.and.(iq.eq.igcm_co_gas)) then
659
660	!! Calculating saturation mixing ratio at surface
661	call cosat(ngrid,ptsrf,pplev(1,1),qsat_co, &
662	pqsurf(:,igcm_n2))
663
664	!! Prevent CO condensation at the surface
665	if (.not.condcosurf) then
666	qsat_co(:)=qsat_co(:)*1.e6
667	endif
668	if (no_n2frost) then
669	qsat_co(:)=qsat_co(:)*qsat_co_factor(:)
670	endif
671
672	DO ig=1,ngrid
673	zd(ig,1)=zb(ig,1)*z1(ig)
674	zq1temp_co(ig)=zc(ig,1)+ zd(ig,1)*qsat_co(ig)
675	pdqsdif(ig,igcm_co_ice)=rho(ig)*zcdv(ig) &
676	*(zq1temp_co(ig)-qsat_co(ig))
677	END DO
678
679
680	DO ig=1,ngrid
681	! -------------------------------------------------------------
682	! TEMPORAIRE : pour initialiser CO si glacier N2
683	! meme si il n'y a pas de CO disponible
684	! if (pqsurf(ig,igcm_n2).le.10.) then
685	! -------------------------------------------------------------
686	!
687	if ((-pdqsdif(ig,igcm_co_ice)*ptimestep) &
688	.gt.(pqsurf(ig,igcm_co_ice))) then
689	pdqsdif(ig,igcm_co_ice)= &
690	-pqsurf(ig,igcm_co_ice)/ptimestep
691	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
692	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_co_gas)+ &
693	zb(ig,2)*zc(ig,2) + &
694	(-pdqsdif(ig,igcm_co_ice)) ptimestep) z1(ig)
695	zq1temp_co(ig)=zc(ig,1)
696	endif
697	! endif
698
699	zq(ig,1,igcm_co_gas)=zq1temp_co(ig)
700
701	! MODIF speciale pour CO / corrected by FF 2016
702	pdtsrf(ig)=pdtsrf(ig)+ &
703	lw_co*pdqsdif(ig,igcm_co_ice)/pcapcal(ig)
704
705	ENDDO ! of DO ig=1,ngrid
706
707	ELSE ! if (methane)
708
709	DO ig=1,ngrid
710	zq(ig,1,iq)=zc(ig,1)
711	ENDDO
712
713	END IF ! of IF ((methane).and.(iq.eq.igcm_ch4_gas))
714
715	!! Diffusion verticale : shut down vertical transport if vertdiff = false
716	if (vertdiff) then
717	DO ilay=2,nlay
718	DO ig=1,ngrid
719	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
720	ENDDO
721	ENDDO
722	endif
723
724	enddo ! of do iq=1,nq
725	end if ! of if(tracer)
726
727	!-----------------------------------------------------------------------
728	! 8. calcul final des tendances de la diffusion verticale
729	! ----------------------------------------------------
730	DO ilev = 1, nlay
731	DO ig=1,ngrid
732	pdudif(ig,ilev)=( zu(ig,ilev)- &
733	(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
734	pdvdif(ig,ilev)=( zv(ig,ilev)- &
735	(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
736	hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep , &
737	zhcond(ig,ilev)) ! modif n2cond
738	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
739	ENDDO
740	ENDDO
741
742	if (tracer) then
743	DO iq = 1, nq
744	DO ilev = 1, nlay
745	DO ig=1,ngrid
746	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)- &
747	(pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
748	! pdqdif(ig,ilev,iq)=pdqdifeddy(ig,ilev,iq)+(zq(ig,ilev,iq)- &
749	! (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
750	ENDDO
751	ENDDO
752	ENDDO
753	end if
754
755	! ** diagnostique final
756	! ------------------
757
758	IF(lecrit) THEN
759	write(,) 'In vdif'
760	write(,) 'Ts (t) and Ts (t+st)'
761	WRITE(*,'(a10,3a15)') &
762	'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
763	write(,) ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
764	DO ilev=1,nlay
765	WRITE(*,'(4f15.7)') &
766	ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev), &
767	pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
768
769	ENDDO
770	ENDIF
771
772	RETURN
773	! END
774
775	end subroutine vdifc_pluto
776
777	end module vdifc_pluto_mod

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

Download in other formats:

Original Format