Context Navigation

vdifc.F @ 132

Last change on this file since 132 was 77, checked in by aslmd, 14 years ago

LMD_MM_MARS: nuages et poussiere radiativement actifs

--> tests concluants
VERSION DE REFERENCE pour NOUVELLE PHYSIQUE
-- il reste le nest a implementer (cf. lignes trop longues)

M 76 mars/libf/phymars/dimradmars.h
M 76 mars/libf/phymars/callradite.F
M 76 mesoscale/LMDZ.MARS.new/myGCM/callphys.def
M 76 mesoscale/TESTS/newphys_tracers/namelist.input
reglages pour nuages et poussiere radiativement actifs

M 76 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
M 76 mars/libf/phymars/meso_physiq.F
correction pour bien sortir les depots de glace d'eau au sol [QSURFICE]
correction pour que pd_scalar soit egal a .true. par defaut

M 76 mars/libf/phymars/vdifc.F
M 76 mesoscale/NOTES.txt
M 76 mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
commentaires et notes

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

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

Download in other formats:

Original Format