Context Navigation

vdifc.F @ 480

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

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

Download in other formats:

Original Format