Context Navigation

vdifc.F @ 513

Last change on this file since 513 was 501, checked in by acolaitis, 14 years ago

More robust gustiness speed computation. Sorry for changing it so often, this time it should be pretty final... This formulation make the 1D model match LES result very well on high vertical resolution. The offset in mixed layer temperature observed on normal resolution is definitely due to vertical resolution (mixing does not occur high enough in the PBL) and not to surface layer underestimating sensible heat flux.

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

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

Download in other formats:

Original Format