Context Navigation

vdifc.F @ 532

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

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

Download in other formats:

Original Format