Context Navigation

vdifc.F @ 288

Last change on this file since 288 was 284, checked in by acolaitis, 13 years ago

--- AC 07/09/2011 ---

Added new flag for the Richardson-based surface layer :

callrichsl, can be changed in callphys.def

One should always use the thermals model when using this surface layer model.
Somes cases (weakly unstable with low winds), when not using thermals, won't be well represented by the
Richardson surface layer. This stands for Mesoscale and Gcm but not for LES model.

Correct configs :

callrichsl = .true.
calltherm = .true.

callrichsl = .false.
calltherm = .false.

callrichsl = .false.
calltherm = .true.

Previously unstable config :

callrichsl = .true.
calltherm = .false.

To be able to run without thermals and with the new surface layer, a modification has been made to

physiq.F to account for gustiness in GCM and MESOSCALE for negative Richardson, so that :

callrichsl = .true.
calltherm = .false.

can now be used without problems, but is not recommended.

Consequently, callrichsl = .false. is now the default configuration for thermals.

We recall the available options in callphys.def for thermals :

outptherm = BOOLEAN (.false. by default) : outputs thermals related quantities (lots of diagfi)
nsplit_thermals = INTEGER (50 by default in gcm, 2 in mesoscale) : subtimestep for thermals model.

It is recommended to use at least 40 in the gcm, and at least 2 in the mesoscale.
The user can lower these values but should check it's log for anomalies or errors regarding
tracer transport in the thermals, or "granulosity" in the outputs for wmax, lmax and hfmax.

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
237	IF (callrichsl .eq. .true.) then
238	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
239	& +(1.wmax(:))*2 ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
240	! LES comparisons. This parameter is linked to the thermals settings)
241	ELSE
242	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
243	ENDIF
244
245	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
246	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
247	! these are the quantities to be looked at when comparing surface layers of different models
248
249	! Some usefull diagnostics for the new surface layer parametrization :
250
251	! call WRITEDIAGFI(ngridmx,'pcdv',
252	! & 'momentum drag','no units',
253	! & 2,zcdv_true)
254	! call WRITEDIAGFI(ngridmx,'pcdh',
255	! & 'heat drag','no units',
256	! & 2,zcdh_true)
257	! call WRITEDIAGFI(ngridmx,'u*',
258	! & 'friction velocity','m/s',
259	! & 2,sqrt(zcdv_true(:))*sqrt(zu2(:)))
260	! call WRITEDIAGFI(ngridmx,'T*',
261	! & 'friction temperature','K',
262	! & 2,-zcdh_true(:)*(ptsrf(:)-ph(:,1))/sqrt(zcdv_true(:)))
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)
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	if (doubleq.AND.submicron) then
497	do ig=1,ngrid
498	c if(co2ice(ig).lt.1) then
499	pdqsdif(ig,igcm_dust_mass) =
500	& -alpha_lift(igcm_dust_mass)
501	pdqsdif(ig,igcm_dust_number) =
502	& -alpha_lift(igcm_dust_number)
503	pdqsdif(ig,igcm_dust_submicron) =
504	& -alpha_lift(igcm_dust_submicron)
505	c end if
506	end do
507	else if (doubleq) then
508	do ig=1,ngrid
509	!!! soulevement constant
510	pdqsdif(ig,igcm_dust_mass) =
511	& -alpha_lift(igcm_dust_mass)
512	pdqsdif(ig,igcm_dust_number) =
513	& -alpha_lift(igcm_dust_number)
514	end do
515	else if (submicron) then
516	do ig=1,ngrid
517	pdqsdif(ig,igcm_dust_submicron) =
518	& -alpha_lift(igcm_dust_submicron)
519	end do
520	else
521	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
522	& pdqsdif)
523	endif !doubleq.AND.submicron
524	else
525	pdqsdif(1:ngrid,1:nq) = 0.
526	end if
527
528	c OU calcul de la valeur de q a la surface (water) :
529	c ----------------------------------------
530	if (water) then
531	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
532	end if
533
534	c Inversion pour l'implicite sur q
535	c --------------------------------
536	do iq=1,nq
537	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
538
539	if ((water).and.(iq.eq.igcm_h2o_vap)) then
540	c This line is required to account for turbulent transport
541	c from surface (e.g. ice) to mid-layer of atmosphere:
542	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
543	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
544	else ! (re)-initialize zb(:,1)
545	zb(1:ngrid,1)=0
546	end if
547
548	DO ig=1,ngrid
549	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
550	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
551	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
552	ENDDO
553
554	DO ilay=nlay-1,2,-1
555	DO ig=1,ngrid
556	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
557	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
558	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
559	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
560	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
561	ENDDO
562	ENDDO
563
564	if (water.and.(iq.eq.igcm_h2o_ice)) then
565	! special case for water ice tracer: do not include
566	! h2o ice tracer from surface (which is set when handling
567	! h2o vapour case (see further down).
568	DO ig=1,ngrid
569	z1(ig)=1./(za(ig,1)+zb(ig,1)+
570	$ zb(ig,2)*(1.-zd(ig,2)))
571	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
572	$ zb(ig,2)zc(ig,2)) z1(ig)
573	ENDDO
574	else ! general case
575	DO ig=1,ngrid
576	z1(ig)=1./(za(ig,1)+zb(ig,1)+
577	$ zb(ig,2)*(1.-zd(ig,2)))
578	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
579	$ zb(ig,2)*zc(ig,2) +
580	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
581	ENDDO
582	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
583
584	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
585	c Calculation for turbulent exchange with the surface (for ice)
586	DO ig=1,ngrid
587	zd(ig,1)=zb(ig,1)*z1(ig)
588	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
589
590	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
591	& *(zq1temp(ig)-qsat(ig))
592	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
593	END DO
594
595	DO ig=1,ngrid
596	if(.not.watercaptag(ig)) then
597	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
598	& .gt.pqsurf(ig,igcm_h2o_ice)) then
599	c write(,)'on sublime plus que qsurf!'
600	pdqsdif(ig,igcm_h2o_ice)=
601	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
602	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
603	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
604	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
605	$ zb(ig,2)*zc(ig,2) +
606	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
607	zq1temp(ig)=zc(ig,1)
608	endif
609	endif ! if (.not.watercaptag(ig))
610	c Starting upward calculations for water :
611	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
612	ENDDO ! of DO ig=1,ngrid
613	ELSE
614	c Starting upward calculations for simple mixing of tracer (dust)
615	DO ig=1,ngrid
616	zq(ig,1,iq)=zc(ig,1)
617	ENDDO
618	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
619
620	DO ilay=2,nlay
621	DO ig=1,ngrid
622	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
623	ENDDO
624	ENDDO
625	enddo ! of do iq=1,nq
626	end if ! of if(tracer)
627
628	c-----------------------------------------------------------------------
629	c 8. calcul final des tendances de la diffusion verticale
630	c ----------------------------------------------------
631
632	DO ilev = 1, nlay
633	DO ig=1,ngrid
634	pdudif(ig,ilev)=( zu(ig,ilev)-
635	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
636	pdvdif(ig,ilev)=( zv(ig,ilev)-
637	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
638	hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
639	$ zhcond(ig,ilev)) ! modif co2cond
640	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
641	ENDDO
642	ENDDO
643
644
645	if (tracer) then
646	DO iq = 1, nq
647	DO ilev = 1, nlay
648	DO ig=1,ngrid
649	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
650	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
651	ENDDO
652	ENDDO
653	ENDDO
654	end if
655
656	c ** diagnostique final
657	c ------------------
658
659	IF(lecrit) THEN
660	PRINT*,'In vdif'
661	PRINT*,'Ts (t) and Ts (t+st)'
662	WRITE(*,'(a10,3a15)')
663	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
664	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
665	DO ilev=1,nlay
666	WRITE(*,'(4f15.7)')
667	s ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
668	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
669
670	ENDDO
671	ENDIF
672
673	if (calltherm .and. outptherm) then
674	if (ngrid .eq. 1) then
675	call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
676	& 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
677	call WRITEDIAGFI(ngrid,'zkh','zkh',
678	& 'SI',1,zkh)
679	endif
680	endif
681
682	RETURN
683	END

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

Download in other formats:

Original Format