Context Navigation

vdifc.F @ 623

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

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

Download in other formats:

Original Format