Context Navigation

vdifc.F @ 3436

Last change on this file since 3436 was 3237, checked in by afalco, 9 months ago
Pluto PCM: Import methane and ch4 clouds from pluto.old AF
File size: 23.9 KB

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

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

Download in other formats:

Original Format