Context Navigation

vdifc.F @ 668

Last change on this file since 668 was 660, checked in by acolaitis, 13 years ago

Thermals: updated computation of wstar to take into account vertical turbulent heat flux from diffusion. Without that, wstar is underestimated. We use the sensible heat flux and tendancies computed in vdifc to compute the vertical turbulent heat flux contribution from the diffusion, and then compute the maximum in the PBL. The resulting total heat flux profile and wstar correspond very well to LES results.

File size: 27.3 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,wstar,zcdv_true,zcdh_true,
8	$ hfmax,sensibFlux
9	#ifdef MESOSCALE
10	& ,flag_LES
11	#endif
12	& )
13	IMPLICIT NONE
14
15	c=======================================================================
16	c
17	c subject:
18	c --------
19	c Turbulent diffusion (mixing) for potential T, U, V and tracer
20	c
21	c Shema implicite
22	c On commence par rajouter au variables x la tendance physique
23	c et on resoult en fait:
24	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
25	c
26	c author:
27	c ------
28	c Hourdin/Forget/Fournier
29	c=======================================================================
30
31	c-----------------------------------------------------------------------
32	c declarations:
33	c -------------
34
35	#include "dimensions.h"
36	#include "dimphys.h"
37	#include "comcstfi.h"
38	#include "callkeys.h"
39	#include "surfdat.h"
40	#include "comgeomfi.h"
41	#include "tracer.h"
42	#include "microphys.h"
43
44	c
45	c arguments:
46	c ----------
47
48	INTEGER ngrid,nlay
49	REAL ptimestep
50	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
51	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
52	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay),pt(ngrid,nlay)
53	REAL ptsrf(ngrid),pemis(ngrid)
54	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
55	REAL pfluxsrf(ngrid)
56	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
57	REAL pdtsrf(ngrid),pcapcal(ngrid)
58	REAL pq2(ngrid,nlay+1)
59
60	c Argument added for condensation:
61	REAL co2ice (ngrid), ppopsk(ngrid,nlay)
62	logical lecrit
63
64	REAL pz0(ngridmx) ! surface roughness length (m)
65
66	c Argument added to account for subgrid gustiness :
67
68	REAL wstar(ngridmx), hfmax(ngridmx), zi(ngridmx)
69
70	c Traceurs :
71	integer nq
72	REAL pqsurf(ngrid,nq)
73	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
74	real pdqdif(ngrid,nlay,nq)
75	real pdqsdif(ngrid,nq)
76
77	c local:
78	c ------
79
80	REAL ust(ngridmx),tst(ngridmx)
81
82	INTEGER ilev,ig,ilay,nlev
83
84	REAL z4st,zdplanck(ngridmx)
85	REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
86	REAL zkq(ngridmx,nlayermx+1)
87	REAL zcdv(ngridmx),zcdh(ngridmx)
88	REAL zcdv_true(ngridmx),zcdh_true(ngridmx) ! drag coeff are used by the LES to recompute u* and hfx
89	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
90	REAL zh(ngridmx,nlayermx)
91	REAL ztsrf2(ngridmx)
92	REAL z1(ngridmx),z2(ngridmx)
93	REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
94	REAL zb0(ngridmx,nlayermx)
95	REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
96	REAL zcst1
97	REAL zu2(ngridmx)
98
99	EXTERNAL SSUM,SCOPY
100	REAL SSUM
101	LOGICAL firstcall
102	SAVE firstcall
103
104	REAL lw
105
106	c variable added for CO2 condensation:
107	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
108	REAL hh , zhcond(ngridmx,nlayermx)
109	REAL latcond,tcond1mb
110	REAL acond,bcond
111	SAVE acond,bcond
112	DATA latcond,tcond1mb/5.9e5,136.27/
113
114	c Tracers :
115	c ~~~~~~~
116	INTEGER iq
117	REAL zq(ngridmx,nlayermx,nqmx)
118	REAL zq1temp(ngridmx)
119	REAL rho(ngridmx) ! near surface air density
120	REAL qsat(ngridmx)
121	DATA firstcall/.true./
122
123	REAL kmixmin
124
125	c Mass-variation scheme :
126	c ~~~~~~~
127
128	INTEGER j,l
129	REAL zcondicea(ngrid,nlayermx)
130	REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1)
131	REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx)
132	REAL pdtc(ngrid,nlayermx)
133	REAL zhs(ngridmx,nlayermx)
134	REAL ccond
135	SAVE ccond
136
137	c Theta_m formulation for mass-variation scheme :
138	c ~~~~~~~
139
140	INTEGER ico2,llnt(ngridmx)
141	SAVE ico2
142	REAL m_co2, m_noco2, A , B
143	SAVE A, B, m_co2, m_noco2
144	REAL vmr_co2(ngridmx,nlayermx)
145	REAL qco2,mmean
146
147	REAL sensibFlux(ngridmx)
148
149	#ifdef MESOSCALE
150	LOGICAL flag_LES !! pour LES avec isfflx!=0
151	#endif
152	c ** un petit test de coherence
153	c --------------------------
154
155	IF (firstcall) THEN
156	IF(ngrid.NE.ngridmx) THEN
157	PRINT*,'STOP dans vdifc'
158	PRINT*,'probleme de dimensions :'
159	PRINT*,'ngrid =',ngrid
160	PRINT*,'ngridmx =',ngridmx
161	STOP
162	ENDIF
163	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
164	bcond=1./tcond1mb
165	acond=r/latcond
166	ccond=cpp/(g*latcond)
167	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
168	PRINT*,' acond,bcond,ccond',acond,bcond,ccond
169
170
171	ico2=0
172
173	if (tracer) then
174	c Prepare Special treatment if one of the tracer is CO2 gas
175	do iq=1,nqmx
176	if (noms(iq).eq."co2") then
177	ico2=iq
178	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
179	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
180	c Compute A and B coefficient use to compute
181	c mean molecular mass Mair defined by
182	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
183	c 1/Mair = A*q(ico2) + B
184	A =(1/m_co2 - 1/m_noco2)
185	B=1/m_noco2
186	endif
187	enddo
188	end if
189
190	firstcall=.false.
191	ENDIF
192
193
194
195
196
197	c-----------------------------------------------------------------------
198	c 1. initialisation
199	c -----------------
200
201	nlev=nlay+1
202
203	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
204	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
205	c ----------------------------------------
206
207	DO ilay=1,nlay
208	DO ig=1,ngrid
209	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
210	! Mass variation scheme:
211	betam(ig,ilay)=-za(ig,ilay)latcond/(cppppopsk(ig,ilay))
212	ENDDO
213	ENDDO
214
215	zcst1=4.gptimestep/(r*r)
216	DO ilev=2,nlev-1
217	DO ig=1,ngrid
218	zb0(ig,ilev)=pplev(ig,ilev)*
219	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
220	s (ph(ig,ilev-1)+ph(ig,ilev))
221	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
222	s (pplay(ig,ilev-1)-pplay(ig,ilev))
223	ENDDO
224	ENDDO
225	DO ig=1,ngrid
226	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
227	ENDDO
228
229	c ** diagnostique pour l'initialisation
230	c ----------------------------------
231
232	IF(lecrit) THEN
233	ig=ngrid/2+1
234	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
235	DO ilay=1,nlay
236	WRITE(*,'(6f11.5)')
237	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
238	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
239	ENDDO
240	PRINT*,'Pression (mbar) ,altitude (km),zb'
241	DO ilev=1,nlay
242	WRITE(*,'(3f15.7)')
243	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
244	s zb0(ig,ilev)
245	ENDDO
246	ENDIF
247
248	c -----------------------------------
249	c Potential Condensation temperature:
250	c -----------------------------------
251
252	c Compute CO2 Volume mixing ratio
253	c -------------------------------
254	if (callcond.and.(ico2.ne.0)) then
255	DO ilev=1,nlay
256	DO ig=1,ngrid
257	qco2=MAX(1.E-30
258	& ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
259	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
260	mmean=1/(A*qco2 +B)
261	vmr_co2(ig,ilev) = qco2*mmean/m_co2
262	ENDDO
263	ENDDO
264	else
265	DO ilev=1,nlay
266	DO ig=1,ngrid
267	vmr_co2(ig,ilev)=0.95
268	ENDDO
269	ENDDO
270	end if
271
272	c forecast of atmospheric temperature zt and frost temperature ztcond
273	c --------------------------------------------------------------------
274
275	DO ilev=1,nlay
276	DO ig=1,ngrid
277	ztcond(ig,ilev)=
278	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
279	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
280	ENDDO
281	ENDDO
282
283	ztcond(:,nlay+1)=ztcond(:,nlay)
284
285	if (callcond) then
286	DO ilev=1,nlay
287	DO ig=1,ngrid
288	! zhcond(ig,ilev) =
289	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
290	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
291	END DO
292	END DO
293	else
294	call zerophys(ngrid*nlay,zhcond)
295	end if
296
297
298	c-----------------------------------------------------------------------
299	c 2. ajout des tendances physiques
300	c -----------------------------
301
302	DO ilev=1,nlay
303	DO ig=1,ngrid
304	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
305	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
306	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
307	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
308	ENDDO
309	ENDDO
310	if(tracer) then
311	DO iq =1, nq
312	DO ilev=1,nlay
313	DO ig=1,ngrid
314	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
315	ENDDO
316	ENDDO
317	ENDDO
318	end if
319
320	c-----------------------------------------------------------------------
321	c 3. schema de turbulence
322	c --------------------
323
324	c ** source d'energie cinetique turbulente a la surface
325	c (condition aux limites du schema de diffusion turbulente
326	c dans la couche limite
327	c ---------------------
328
329	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
330	& ,zcdv_true,zcdh_true)
331
332	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
333
334	IF (callrichsl) THEN
335	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
336	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
337	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
338	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
339
340	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
341	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
342
343	tst(:)=0.
344	DO ig=1,ngrid
345	IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u=t=0
346	tst(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ust(ig)
347	ENDIF
348	ENDDO
349
350	ELSE
351	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
352	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
353	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
354	tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
355	ENDIF
356
357	! Some usefull diagnostics for the new surface layer parametrization :
358
359	! call WRITEDIAGFI(ngridmx,'pcdv',
360	! & 'momentum drag','no units',
361	! & 2,zcdv_true)
362	! call WRITEDIAGFI(ngridmx,'pcdh',
363	! & 'heat drag','no units',
364	! & 2,zcdh_true)
365	! call WRITEDIAGFI(ngridmx,'ust',
366	! & 'friction velocity','m/s',2,ust)
367	! call WRITEDIAGFI(ngridmx,'tst',
368	! & 'friction temperature','K',2,tst)
369	! call WRITEDIAGFI(ngridmx,'rm-1',
370	! & 'aerodyn momentum conductance','m/s',
371	! & 2,zcdv)
372	! call WRITEDIAGFI(ngridmx,'rh-1',
373	! & 'aerodyn heat conductance','m/s',
374	! & 2,zcdh)
375
376	c ** schema de diffusion turbulente dans la couche limite
377	c ----------------------------------------------------
378	IF (.not. calltherm) THEN
379
380	CALL vdif_kc(ptimestep,g,pzlev,pzlay
381	& ,pu,pv,ph,zcdv_true
382	& ,pq2,zkv,zkh,zq)
383
384	ELSE
385
386	pt(:,:)=ph(:,:)*ppopsk(:,:)
387	CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
388	s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ust
389	s ,9)
390
391	ENDIF
392
393	if ((doubleq).and.(ngrid.eq.1)) then
394	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
395	do ilev=1,nlay
396	do ig=1,ngrid
397	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
398	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
399	end do
400	end do
401	end if
402
403	c ** diagnostique pour le schema de turbulence
404	c -----------------------------------------
405
406	IF(lecrit) THEN
407	PRINT*
408	PRINT*,'Diagnostic for the vertical turbulent mixing'
409	PRINT*,'Cd for momentum and potential temperature'
410
411	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
412	PRINT*,'Mixing coefficient for momentum and pot.temp.'
413	DO ilev=1,nlay
414	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
415	ENDDO
416	ENDIF
417
418
419
420
421	c-----------------------------------------------------------------------
422	c 4. inversion pour l'implicite sur u
423	c --------------------------------
424
425	c ** l'equation est
426	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
427	c avec
428	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
429	c et
430	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
431	c donc les entrees sont /zcdv/ pour la condition a la limite sol
432	c et /zkv/ = Ku
433
434	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
435	CALL multipl(ngrid,zcdv,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)zu(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)*zu(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	DO ig=1,ngrid
454	zu(ig,1)=zc(ig,1)
455	ENDDO
456	DO ilay=2,nlay
457	DO ig=1,ngrid
458	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
459	ENDDO
460	ENDDO
461
462
463
464
465
466	c-----------------------------------------------------------------------
467	c 5. inversion pour l'implicite sur v
468	c --------------------------------
469
470	c ** l'equation est
471	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
472	c avec
473	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
474	c et
475	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
476	c donc les entrees sont /zcdv/ pour la condition a la limite sol
477	c et /zkv/ = Kv
478
479	DO ig=1,ngrid
480	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
481	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
482	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
483	ENDDO
484
485	DO ilay=nlay-1,1,-1
486	DO ig=1,ngrid
487	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
488	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
489	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
490	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
491	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
492	ENDDO
493	ENDDO
494
495	DO ig=1,ngrid
496	zv(ig,1)=zc(ig,1)
497	ENDDO
498	DO ilay=2,nlay
499	DO ig=1,ngrid
500	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
501	ENDDO
502	ENDDO
503
504
505
506
507
508	c-----------------------------------------------------------------------
509	c 6. inversion pour l'implicite sur h sans oublier le couplage
510	c avec le sol (conduction)
511	c ------------------------
512
513	c ** l'equation est
514	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
515	c avec
516	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
517	c et
518	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
519	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
520	c et /zkh/ = Kh
521	c -------------
522
523	c Mass variation scheme:
524	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
525	CALL multipl(ngrid,zcdh,zb0,zb)
526
527	c on initialise dm c
528
529	pdtc(:,:)=0.
530	zt(:,:)=0.
531	dmice(:,:)=0.
532
533	c ** calcul de (d Planck / dT) a la temperature d'interface
534	c ------------------------------------------------------
535
536	z4st=4.5.67e-8ptimestep
537	IF (tke_heat_flux .eq. 0.) THEN
538	DO ig=1,ngrid
539	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
540	ENDDO
541	ELSE
542	zdplanck(:)=0.
543	ENDIF
544
545	! calcul de zc et zd pour la couche top en prenant en compte le terme
546	! de variation de masse (on fait une boucle pour que ça converge)
547
548	! Identification des points de grilles qui ont besoin de la correction
549
550	llnt(:)=1
551	#ifdef MESOSCALE
552	IF (.not.flag_LES) THEN
553	#endif
554	DO ig=1,ngrid
555	DO l=1,nlay
556	if(zh(ig,l) .lt. zhcond(ig,l)) then
557	llnt(ig)=300
558	! 200 and 100 do not go beyond month 9 with normal dissipation
559	goto 5
560	endif
561	ENDDO
562	5 continue
563	ENDDO
564
565	#ifdef MESOSCALE
566	ENDIF
567	#endif
568
569	DO ig=1,ngrid
570
571	! Initialization of z1 and zd, which do not depend on dmice
572
573	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
574	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
575
576	DO ilay=nlay-1,1,-1
577	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
578	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
579	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
580	ENDDO
581
582	! Convergence loop
583
584	DO j=1,llnt(ig)
585
586	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
587	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
588	& -betam(ig,nlay)*dmice(ig,nlay)
589	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
590	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
591
592	! calcul de zc et zd pour les couches du haut vers le bas
593
594	DO ilay=nlay-1,1,-1
595	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
596	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
597	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
598	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
599	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
600	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
601	ENDDO
602
603	c ** calcul de la temperature_d'interface et de sa tendance.
604	c on ecrit que la somme des flux est nulle a l'interface
605	c a t + \delta t,
606	c c'est a dire le flux radiatif a {t + \delta t}
607	c + le flux turbulent a {t + \delta t}
608	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
609	c (notation K dt = /cpp*b/)
610	c + le flux dans le sol a t
611	c + l'evolution du flux dans le sol lorsque la temperature d'interface
612	c passe de sa valeur a t a sa valeur a {t + \delta t}.
613	c ----------------------------------------------------
614
615	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
616	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
617	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
618	ztsrf2(ig)=z1(ig)/z2(ig)
619	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
620	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
621
622	c ** et a partir de la temperature au sol on remonte
623	c -----------------------------------------------
624
625	DO ilay=2,nlay
626	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
627	ENDDO
628	DO ilay=1,nlay
629	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
630	ENDDO
631
632	c Condensation/sublimation in the atmosphere
633	c ------------------------------------------
634	c (computation of zcondicea and dmice)
635
636	zcondicea(ig,:)=0.
637	pdtc(ig,:)=0.
638
639	DO l=nlay , 1, -1
640	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
641	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
642	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
643	& ccondpdtc(ig,l)
644	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
645	END IF
646	ENDDO
647
648	ENDDO !of Do j=1,XXX
649
650	ENDDO !of Do ig=1,nlay
651
652	pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
653
654	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
655	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zhs(ig,1)-ztsrf2(ig))
656	ENDDO
657
658	c-----------------------------------------------------------------------
659	c TRACERS
660	c -------
661
662	if(tracer) then
663
664	c Using the wind modified by friction for lifting and sublimation
665	c ----------------------------------------------------------------
666
667	! This is computed above and takes into account surface-atmosphere flux
668	! enhancement by subgrid gustiness and atmospheric-stability related
669	! variations of transfer coefficients.
670
671	! DO ig=1,ngrid
672	! zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
673	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
674	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
675	! ENDDO
676
677	c Calcul du flux vertical au bas de la premiere couche (dust) :
678	c -----------------------------------------------------------
679	do ig=1,ngridmx
680	rho(ig) = zb0(ig,1) /ptimestep
681	c zb(ig,1) = 0.
682	end do
683	c Dust lifting:
684	if (lifting) then
685	#ifndef MESOSCALE
686	if (doubleq.AND.submicron) then
687	do ig=1,ngrid
688	c if(co2ice(ig).lt.1) then
689	pdqsdif(ig,igcm_dust_mass) =
690	& -alpha_lift(igcm_dust_mass)
691	pdqsdif(ig,igcm_dust_number) =
692	& -alpha_lift(igcm_dust_number)
693	pdqsdif(ig,igcm_dust_submicron) =
694	& -alpha_lift(igcm_dust_submicron)
695	c end if
696	end do
697	else if (doubleq) then
698	do ig=1,ngrid
699	if(co2ice(ig).lt.1) then ! soulevement pas constant
700	pdqsdif(ig,igcm_dust_mass) =
701	& -alpha_lift(igcm_dust_mass)
702	pdqsdif(ig,igcm_dust_number) =
703	& -alpha_lift(igcm_dust_number)
704	end if
705	end do
706	else if (submicron) then
707	do ig=1,ngrid
708	pdqsdif(ig,igcm_dust_submicron) =
709	& -alpha_lift(igcm_dust_submicron)
710	end do
711	else
712	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
713	& pdqsdif)
714	endif !doubleq.AND.submicron
715	#else
716	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
717	& pdqsdif)
718	#endif
719	else
720	pdqsdif(1:ngrid,1:nq) = 0.
721	end if
722
723	c OU calcul de la valeur de q a la surface (water) :
724	c ----------------------------------------
725	if (water) then
726	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
727	end if
728
729	c Inversion pour l'implicite sur q
730	c --------------------------------
731	do iq=1,nq
732	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
733
734	if ((water).and.(iq.eq.igcm_h2o_vap)) then
735	c This line is required to account for turbulent transport
736	c from surface (e.g. ice) to mid-layer of atmosphere:
737	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
738	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
739	else ! (re)-initialize zb(:,1)
740	zb(1:ngrid,1)=0
741	end if
742
743	DO ig=1,ngrid
744	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
745	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
746	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
747	ENDDO
748
749	DO ilay=nlay-1,2,-1
750	DO ig=1,ngrid
751	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
752	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
753	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
754	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
755	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
756	ENDDO
757	ENDDO
758
759	if (water.and.(iq.eq.igcm_h2o_ice)) then
760	! special case for water ice tracer: do not include
761	! h2o ice tracer from surface (which is set when handling
762	! h2o vapour case (see further down).
763	DO ig=1,ngrid
764	z1(ig)=1./(za(ig,1)+zb(ig,1)+
765	$ zb(ig,2)*(1.-zd(ig,2)))
766	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
767	$ zb(ig,2)zc(ig,2)) z1(ig)
768	ENDDO
769	else ! general case
770	DO ig=1,ngrid
771	z1(ig)=1./(za(ig,1)+zb(ig,1)+
772	$ zb(ig,2)*(1.-zd(ig,2)))
773	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
774	$ zb(ig,2)*zc(ig,2) +
775	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
776	ENDDO
777	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
778
779	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
780	c Calculation for turbulent exchange with the surface (for ice)
781	DO ig=1,ngrid
782	zd(ig,1)=zb(ig,1)*z1(ig)
783	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
784
785	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
786	& *(zq1temp(ig)-qsat(ig))
787	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
788	END DO
789
790	DO ig=1,ngrid
791	if(.not.watercaptag(ig)) then
792	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
793	& .gt.pqsurf(ig,igcm_h2o_ice)) then
794	c write(,)'on sublime plus que qsurf!'
795	pdqsdif(ig,igcm_h2o_ice)=
796	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
797	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
798	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
799	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
800	$ zb(ig,2)*zc(ig,2) +
801	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
802	zq1temp(ig)=zc(ig,1)
803	endif
804	endif ! if (.not.watercaptag(ig))
805	c Starting upward calculations for water :
806	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
807	c Take into account H2o latent heat in surface energy budget
808	lw = (2834.3-0.28*(ptsrf(ig)-To)
809	& -0.004(ptsrf(ig)-To)2)1.e+3
810	pdtsrf(ig) = pdtsrf(ig)
811	& + pdqsdif(ig,igcm_h2o_ice)*lw/pcapcal(ig)
812	ENDDO ! of DO ig=1,ngrid
813	ELSE
814	c Starting upward calculations for simple mixing of tracer (dust)
815	DO ig=1,ngrid
816	zq(ig,1,iq)=zc(ig,1)
817	ENDDO
818	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
819
820	DO ilay=2,nlay
821	DO ig=1,ngrid
822	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
823	ENDDO
824	ENDDO
825	enddo ! of do iq=1,nq
826	end if ! of if(tracer)
827
828	c-----------------------------------------------------------------------
829	c 8. calcul final des tendances de la diffusion verticale
830	c ----------------------------------------------------
831
832	DO ilev = 1, nlay
833	DO ig=1,ngrid
834	pdudif(ig,ilev)=( zu(ig,ilev)-
835	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
836	pdvdif(ig,ilev)=( zv(ig,ilev)-
837	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
838	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
839	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
840	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
841	ENDDO
842	ENDDO
843
844	if (tracer) then
845	DO iq = 1, nq
846	DO ilev = 1, nlay
847	DO ig=1,ngrid
848	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
849	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
850	ENDDO
851	ENDDO
852	ENDDO
853	end if
854
855	c ** diagnostique final
856	c ------------------
857
858	IF(lecrit) THEN
859	PRINT*,'In vdif'
860	PRINT*,'Ts (t) and Ts (t+st)'
861	WRITE(*,'(a10,3a15)')
862	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
863	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
864	DO ilev=1,nlay
865	WRITE(*,'(4f15.7)')
866	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
867	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
868
869	ENDDO
870	ENDIF
871
872	RETURN
873	END

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

Download in other formats:

Original Format