Context Navigation

vdifc.F @ 469

Last change on this file since 469 was 325, checked in by acolaitis, 14 years ago
Included variation of mean molar mass in potential temperature definition: this modification realistically accounts for mixing in the polar night when the atm is enriched in Ar. This follows present computations in convadj for the polar night case.
File size: 22.0 KB

Line
1	SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk,
2	$ ptimestep,pcapcal,lecrit,
3	$ pplay,pplev,pzlay,pzlev,pz0,
4	$ pu,pv,ph,pq,ptsrf,pemis,pqsurf,
5	$ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
6	$ pdudif,pdvdif,pdhdif,pdtsrf,pq2,
7	$ pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
8	IMPLICIT NONE
9
10	c=======================================================================
11	c
12	c subject:
13	c --------
14	c Turbulent diffusion (mixing) for potential T, U, V and tracer
15	c
16	c Shema implicite
17	c On commence par rajouter au variables x la tendance physique
18	c et on resoult en fait:
19	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
20	c
21	c author:
22	c ------
23	c Hourdin/Forget/Fournier
24	c=======================================================================
25
26	c-----------------------------------------------------------------------
27	c declarations:
28	c -------------
29
30	#include "dimensions.h"
31	#include "dimphys.h"
32	#include "comcstfi.h"
33	#include "callkeys.h"
34	#include "surfdat.h"
35	#include "comgeomfi.h"
36	#include "tracer.h"
37
38	c
39	c arguments:
40	c ----------
41
42	INTEGER ngrid,nlay
43	REAL ptimestep
44	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
45	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
46	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
47	REAL ptsrf(ngrid),pemis(ngrid)
48	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
49	REAL pfluxsrf(ngrid)
50	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
51	REAL pdtsrf(ngrid),pcapcal(ngrid)
52	REAL pq2(ngrid,nlay+1)
53
54	c Argument added for condensation:
55	REAL co2ice (ngrid), ppopsk(ngrid,nlay)
56	logical lecrit
57
58	REAL pz0(ngridmx) ! surface roughness length (m)
59
60	c Argument added to account for subgrid gustiness :
61
62	REAL wmax(ngridmx)
63
64	c Traceurs :
65	integer nq
66	REAL pqsurf(ngrid,nq)
67	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
68	real pdqdif(ngrid,nlay,nq)
69	real pdqsdif(ngrid,nq)
70
71	c local:
72	c ------
73
74	INTEGER ilev,ig,ilay,nlev
75
76	REAL z4st,zdplanck(ngridmx)
77	REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
78	REAL zcdv(ngridmx),zcdh(ngridmx)
79	REAL zcdv_true(ngridmx),zcdh_true(ngridmx) ! drag coeff are used by the LES to recompute u* and hfx
80	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
81	REAL zh(ngridmx,nlayermx)
82	REAL ztsrf2(ngridmx)
83	REAL z1(ngridmx),z2(ngridmx)
84	REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
85	REAL zb0(ngridmx,nlayermx)
86	REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
87	REAL zcst1
88	REAL zu2(ngridmx)
89
90	EXTERNAL SSUM,SCOPY
91	REAL SSUM
92	LOGICAL firstcall
93	SAVE firstcall
94
95	c variable added for CO2 condensation:
96	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97	REAL hh , zhcond(ngridmx,nlayermx)
98	REAL latcond,tcond1mb
99	REAL acond,bcond
100	SAVE acond,bcond
101	DATA latcond,tcond1mb/5.9e5,136.27/
102
103	c Tracers :
104	c ~~~~~~~
105	INTEGER iq
106	REAL zq(ngridmx,nlayermx,nqmx)
107	REAL zq1temp(ngridmx)
108	REAL rho(ngridmx) ! near surface air density
109	REAL qsat(ngridmx)
110	DATA firstcall/.true./
111
112	REAL kmixmin
113
114	c ** un petit test de coherence
115	c --------------------------
116
117	IF (firstcall) THEN
118	IF(ngrid.NE.ngridmx) THEN
119	PRINT*,'STOP dans vdifc'
120	PRINT*,'probleme de dimensions :'
121	PRINT*,'ngrid =',ngrid
122	PRINT*,'ngridmx =',ngridmx
123	STOP
124	ENDIF
125	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
126	bcond=1./tcond1mb
127	acond=r/latcond
128	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
129	PRINT*,' acond,bcond',acond,bcond
130
131	firstcall=.false.
132	ENDIF
133
134
135
136
137
138	c-----------------------------------------------------------------------
139	c 1. initialisation
140	c -----------------
141
142	nlev=nlay+1
143
144	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
145	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
146	c ----------------------------------------
147
148	DO ilay=1,nlay
149	DO ig=1,ngrid
150	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
151	ENDDO
152	ENDDO
153
154	zcst1=4.gptimestep/(r*r)
155	DO ilev=2,nlev-1
156	DO ig=1,ngrid
157	zb0(ig,ilev)=pplev(ig,ilev)*
158	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
159	s (ph(ig,ilev-1)+ph(ig,ilev))
160	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
161	s (pplay(ig,ilev-1)-pplay(ig,ilev))
162	ENDDO
163	ENDDO
164	DO ig=1,ngrid
165	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
166	ENDDO
167
168	c ** diagnostique pour l'initialisation
169	c ----------------------------------
170
171	IF(lecrit) THEN
172	ig=ngrid/2+1
173	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
174	DO ilay=1,nlay
175	WRITE(*,'(6f11.5)')
176	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
177	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
178	ENDDO
179	PRINT*,'Pression (mbar) ,altitude (km),zb'
180	DO ilev=1,nlay
181	WRITE(*,'(3f15.7)')
182	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
183	s zb0(ig,ilev)
184	ENDDO
185	ENDIF
186
187	c Potential Condensation temperature:
188	c -----------------------------------
189
190	c if (callcond) then
191	c DO ilev=1,nlay
192	c DO ig=1,ngrid
193	c zhcond(ig,ilev) =
194	c & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
195	c END DO
196	c END DO
197	c else
198	call zerophys(ngrid*nlay,zhcond)
199	c end if
200
201
202	c-----------------------------------------------------------------------
203	c 2. ajout des tendances physiques
204	c -----------------------------
205
206	DO ilev=1,nlay
207	DO ig=1,ngrid
208	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
209	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
210	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
211	zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
212	ENDDO
213	ENDDO
214	if(tracer) then
215	DO iq =1, nq
216	DO ilev=1,nlay
217	DO ig=1,ngrid
218	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
219	ENDDO
220	ENDDO
221	ENDDO
222	end if
223
224	c-----------------------------------------------------------------------
225	c 3. schema de turbulence
226	c --------------------
227
228	c ** source d'energie cinetique turbulente a la surface
229	c (condition aux limites du schema de diffusion turbulente
230	c dans la couche limite
231	c ---------------------
232
233	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wmax,ptsrf,ph
234	& ,zcdv_true,zcdh_true)
235
236	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
237
238	IF (callrichsl) THEN
239	zcdv(:)=zcdv_true(:)sqrt(zu2(:)+(0.7wmax(:))**2) ! subgrid gustiness added by quadratic interpolation with a coeff beta
240	zcdh(:)=zcdh_true(:)sqrt(zu2(:)+(1.2wmax(:))**2) ! LES comparisons. This parameter is linked to the thermals settings)
241	ELSE
242	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
243	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
244	ENDIF
245
246
247	! Some usefull diagnostics for the new surface layer parametrization :
248
249	! call WRITEDIAGFI(ngridmx,'pcdv',
250	! & 'momentum drag','no units',
251	! & 2,zcdv_true)
252	! call WRITEDIAGFI(ngridmx,'pcdh',
253	! & 'heat drag','no units',
254	! & 2,zcdh_true)
255	! call WRITEDIAGFI(ngridmx,'u*',
256	! & 'friction velocity','m/s',
257	! & 2,sqrt(zcdv_true(:))sqrt(zu2(:)+(0.7wmax(:))**2))
258	! call WRITEDIAGFI(ngridmx,'T*',
259	! & 'friction temperature','K',
260	! & 2,zcdh_true(:)sqrt(zu2(:)+(1.2wmax(:))*2)
261	! & -(ptsrf(:)-ph(:,1))
262	! & /(sqrt(zcdv_true(:))sqrt(zu2(:)+(0.7wmax(:))**2)))
263	! call WRITEDIAGFI(ngridmx,'rm-1',
264	! & 'aerodyn momentum conductance','m/s',
265	! & 2,zcdv)
266	! call WRITEDIAGFI(ngridmx,'rh-1',
267	! & 'aerodyn heat conductance','m/s',
268	! & 2,zcdh)
269
270	c ** schema de diffusion turbulente dans la couche limite
271	c ----------------------------------------------------
272
273	CALL vdif_kc(ptimestep,g,pzlev,pzlay
274	& ,pu,pv,ph,zcdv_true
275	& ,pq2,zkv,zkh,zq)
276
277	if ((doubleq).and.(ngrid.eq.1)) then
278	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
279	do ilev=1,nlay
280	do ig=1,ngrid
281	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
282	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
283	end do
284	end do
285	end if
286
287	c ** diagnostique pour le schema de turbulence
288	c -----------------------------------------
289
290	IF(lecrit) THEN
291	PRINT*
292	PRINT*,'Diagnostic for the vertical turbulent mixing'
293	PRINT*,'Cd for momentum and potential temperature'
294
295	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
296	PRINT*,'Mixing coefficient for momentum and pot.temp.'
297	DO ilev=1,nlay
298	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
299	ENDDO
300	ENDIF
301
302
303
304
305	c-----------------------------------------------------------------------
306	c 4. inversion pour l'implicite sur u
307	c --------------------------------
308
309	c ** l'equation est
310	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
311	c avec
312	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
313	c et
314	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
315	c donc les entrees sont /zcdv/ pour la condition a la limite sol
316	c et /zkv/ = Ku
317
318	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
319	CALL multipl(ngrid,zcdv,zb0,zb)
320
321	DO ig=1,ngrid
322	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
323	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
324	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
325	ENDDO
326
327	DO ilay=nlay-1,1,-1
328	DO ig=1,ngrid
329	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
330	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
331	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
332	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
333	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
334	ENDDO
335	ENDDO
336
337	DO ig=1,ngrid
338	zu(ig,1)=zc(ig,1)
339	ENDDO
340	DO ilay=2,nlay
341	DO ig=1,ngrid
342	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
343	ENDDO
344	ENDDO
345
346
347
348
349
350	c-----------------------------------------------------------------------
351	c 5. inversion pour l'implicite sur v
352	c --------------------------------
353
354	c ** l'equation est
355	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
356	c avec
357	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
358	c et
359	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
360	c donc les entrees sont /zcdv/ pour la condition a la limite sol
361	c et /zkv/ = Kv
362
363	DO ig=1,ngrid
364	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
365	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
366	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
367	ENDDO
368
369	DO ilay=nlay-1,1,-1
370	DO ig=1,ngrid
371	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
372	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
373	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
374	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
375	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
376	ENDDO
377	ENDDO
378
379	DO ig=1,ngrid
380	zv(ig,1)=zc(ig,1)
381	ENDDO
382	DO ilay=2,nlay
383	DO ig=1,ngrid
384	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
385	ENDDO
386	ENDDO
387
388
389
390
391
392	c-----------------------------------------------------------------------
393	c 6. inversion pour l'implicite sur h sans oublier le couplage
394	c avec le sol (conduction)
395	c ------------------------
396
397	c ** l'equation est
398	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
399	c avec
400	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
401	c et
402	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
403	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
404	c et /zkh/ = Kh
405	c -------------
406
407	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
408	CALL multipl(ngrid,zcdh,zb0,zb)
409
410	DO ig=1,ngrid
411	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
412	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
413	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
414	ENDDO
415
416	DO ilay=nlay-1,1,-1
417	DO ig=1,ngrid
418	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
419	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
420	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
421	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
422	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
423	ENDDO
424	ENDDO
425
426	c ** calcul de (d Planck / dT) a la temperature d'interface
427	c ------------------------------------------------------
428
429	z4st=4.5.67e-8ptimestep
430	DO ig=1,ngrid
431	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
432	ENDDO
433
434	c ** calcul de la temperature_d'interface et de sa tendance.
435	c on ecrit que la somme des flux est nulle a l'interface
436	c a t + \delta t,
437	c c'est a dire le flux radiatif a {t + \delta t}
438	c + le flux turbulent a {t + \delta t}
439	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
440	c (notation K dt = /cpp*b/)
441	c + le flux dans le sol a t
442	c + l'evolution du flux dans le sol lorsque la temperature d'interface
443	c passe de sa valeur a t a sa valeur a {t + \delta t}.
444	c ----------------------------------------------------
445
446	DO ig=1,ngrid
447	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
448	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
449	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
450	ztsrf2(ig)=z1(ig)/z2(ig)
451	pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
452
453	c Modif speciale CO2 condensation:
454	c tconds = 1./(bcond-acondlog(.0095pplev(ig,1)))
455	c if ((callcond).and.
456	c & ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
457	c zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
458	c else
459	zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
460	c end if
461	ENDDO
462
463	c ** et a partir de la temperature au sol on remonte
464	c -----------------------------------------------
465
466	DO ilay=2,nlay
467	DO ig=1,ngrid
468	hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond
469	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
470	ENDDO
471	ENDDO
472
473
474	c-----------------------------------------------------------------------
475	c TRACERS
476	c -------
477
478	if(tracer) then
479
480	c Using the wind modified by friction for lifting and sublimation
481	c ----------------------------------------------------------------
482	DO ig=1,ngrid
483	zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
484	zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
485	zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
486	ENDDO
487
488	c Calcul du flux vertical au bas de la premiere couche (dust) :
489	c -----------------------------------------------------------
490	do ig=1,ngridmx
491	rho(ig) = zb0(ig,1) /ptimestep
492	c zb(ig,1) = 0.
493	end do
494	c Dust lifting:
495	if (lifting) then
496	#ifndef MESOSCALE
497	if (doubleq.AND.submicron) then
498	do ig=1,ngrid
499	c if(co2ice(ig).lt.1) then
500	pdqsdif(ig,igcm_dust_mass) =
501	& -alpha_lift(igcm_dust_mass)
502	pdqsdif(ig,igcm_dust_number) =
503	& -alpha_lift(igcm_dust_number)
504	pdqsdif(ig,igcm_dust_submicron) =
505	& -alpha_lift(igcm_dust_submicron)
506	c end if
507	end do
508	else if (doubleq) then
509	do ig=1,ngrid
510	!!! soulevement constant
511	pdqsdif(ig,igcm_dust_mass) =
512	& -alpha_lift(igcm_dust_mass)
513	pdqsdif(ig,igcm_dust_number) =
514	& -alpha_lift(igcm_dust_number)
515	end do
516	else if (submicron) then
517	do ig=1,ngrid
518	pdqsdif(ig,igcm_dust_submicron) =
519	& -alpha_lift(igcm_dust_submicron)
520	end do
521	else
522	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
523	& pdqsdif)
524	endif !doubleq.AND.submicron
525	#else
526	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
527	& pdqsdif)
528	#endif
529	else
530	pdqsdif(1:ngrid,1:nq) = 0.
531	end if
532
533	c OU calcul de la valeur de q a la surface (water) :
534	c ----------------------------------------
535	if (water) then
536	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
537	end if
538
539	c Inversion pour l'implicite sur q
540	c --------------------------------
541	do iq=1,nq
542	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
543
544	if ((water).and.(iq.eq.igcm_h2o_vap)) then
545	c This line is required to account for turbulent transport
546	c from surface (e.g. ice) to mid-layer of atmosphere:
547	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
548	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
549	else ! (re)-initialize zb(:,1)
550	zb(1:ngrid,1)=0
551	end if
552
553	DO ig=1,ngrid
554	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
555	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
556	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
557	ENDDO
558
559	DO ilay=nlay-1,2,-1
560	DO ig=1,ngrid
561	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
562	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
563	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
564	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
565	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
566	ENDDO
567	ENDDO
568
569	if (water.and.(iq.eq.igcm_h2o_ice)) then
570	! special case for water ice tracer: do not include
571	! h2o ice tracer from surface (which is set when handling
572	! h2o vapour case (see further down).
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	else ! general case
580	DO ig=1,ngrid
581	z1(ig)=1./(za(ig,1)+zb(ig,1)+
582	$ zb(ig,2)*(1.-zd(ig,2)))
583	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
584	$ zb(ig,2)*zc(ig,2) +
585	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
586	ENDDO
587	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
588
589	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
590	c Calculation for turbulent exchange with the surface (for ice)
591	DO ig=1,ngrid
592	zd(ig,1)=zb(ig,1)*z1(ig)
593	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
594
595	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
596	& *(zq1temp(ig)-qsat(ig))
597	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
598	END DO
599
600	DO ig=1,ngrid
601	if(.not.watercaptag(ig)) then
602	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
603	& .gt.pqsurf(ig,igcm_h2o_ice)) then
604	c write(,)'on sublime plus que qsurf!'
605	pdqsdif(ig,igcm_h2o_ice)=
606	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
607	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
608	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
609	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
610	$ zb(ig,2)*zc(ig,2) +
611	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
612	zq1temp(ig)=zc(ig,1)
613	endif
614	endif ! if (.not.watercaptag(ig))
615	c Starting upward calculations for water :
616	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
617	ENDDO ! of DO ig=1,ngrid
618	ELSE
619	c Starting upward calculations for simple mixing of tracer (dust)
620	DO ig=1,ngrid
621	zq(ig,1,iq)=zc(ig,1)
622	ENDDO
623	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
624
625	DO ilay=2,nlay
626	DO ig=1,ngrid
627	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
628	ENDDO
629	ENDDO
630	enddo ! of do iq=1,nq
631	end if ! of if(tracer)
632
633	c-----------------------------------------------------------------------
634	c 8. calcul final des tendances de la diffusion verticale
635	c ----------------------------------------------------
636
637	DO ilev = 1, nlay
638	DO ig=1,ngrid
639	pdudif(ig,ilev)=( zu(ig,ilev)-
640	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
641	pdvdif(ig,ilev)=( zv(ig,ilev)-
642	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
643	hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
644	$ zhcond(ig,ilev)) ! modif co2cond
645	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
646	ENDDO
647	ENDDO
648
649
650	if (tracer) then
651	DO iq = 1, nq
652	DO ilev = 1, nlay
653	DO ig=1,ngrid
654	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
655	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
656	ENDDO
657	ENDDO
658	ENDDO
659	end if
660
661	c ** diagnostique final
662	c ------------------
663
664	IF(lecrit) THEN
665	PRINT*,'In vdif'
666	PRINT*,'Ts (t) and Ts (t+st)'
667	WRITE(*,'(a10,3a15)')
668	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
669	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
670	DO ilev=1,nlay
671	WRITE(*,'(4f15.7)')
672	s ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
673	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
674
675	ENDDO
676	ENDIF
677
678	if (calltherm .and. outptherm) then
679	if (ngrid .eq. 1) then
680	call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
681	& 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
682	call WRITEDIAGFI(ngrid,'zkh','zkh',
683	& 'SI',1,zkh)
684	endif
685	endif
686
687	RETURN
688	END

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

Download in other formats:

Original Format