Context Navigation

vdifc.F @ 499

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

M 498 libf/phymars/calltherm_interface.F90
^{------------------> Added computation of wstar. Wmax is no longer an output.}

M 498 libf/phymars/physiq.F
^{------------------> Minor changes.}

M 498 libf/phymars/vdif_cd.F
^{------------------> Changed computation of bulk Richardson number to allow more unstable configurations.}

M 498 libf/phymars/vdifc.F
^{------------------> Changed computation of surface conductance coefficients following a new approach. Comparisons with LES are much better.}

File size: 26.1 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	& (14.75sqrt(hfmax(:)) - 10.hfmax(:) + 1.7hfmax(:)2)*2) ! subgrid gustiness added by quadratic interpolation with a coeff beta
323	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
324	& (20.6hfmax(:) - 22.8hfmax(:)*2 + 9.hfmax(:)4)2) ! LES comparisons. This parameter is linked to the thermals settings)
325
326	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
327	& (14.75sqrt(hfmax(:)) - 10.hfmax(:) + 1.7hfmax(:)2)*2)
328	tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
329	! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
330	! tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
331	ELSE
332	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
333	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
334	ENDIF
335
336
337	! Some usefull diagnostics for the new surface layer parametrization :
338
339	! call WRITEDIAGFI(ngridmx,'pcdv',
340	! & 'momentum drag','no units',
341	! & 2,zcdv_true)
342	! call WRITEDIAGFI(ngridmx,'pcdh',
343	! & 'heat drag','no units',
344	! & 2,zcdh_true)
345	! call WRITEDIAGFI(ngridmx,'ust',
346	! & 'friction velocity','m/s',2,ust)
347	! call WRITEDIAGFI(ngridmx,'tst',
348	! & 'friction temperature','K',2,tst)
349	! call WRITEDIAGFI(ngridmx,'rm-1',
350	! & 'aerodyn momentum conductance','m/s',
351	! & 2,zcdv)
352	! call WRITEDIAGFI(ngridmx,'rh-1',
353	! & 'aerodyn heat conductance','m/s',
354	! & 2,zcdh)
355
356	c ** schema de diffusion turbulente dans la couche limite
357	c ----------------------------------------------------
358	!
359	CALL vdif_kc(ptimestep,g,pzlev,pzlay
360	& ,pu,pv,ph,zcdv_true
361	& ,pq2,zkv,zkh,zq)
362	!
363	! pt(:,:)=ph(:,:)*ppopsk(:,:)
364	! CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
365	! s ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
366	! s ,8)
367
368	if ((doubleq).and.(ngrid.eq.1)) then
369	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
370	do ilev=1,nlay
371	do ig=1,ngrid
372	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
373	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
374	end do
375	end do
376	end if
377
378	c ** diagnostique pour le schema de turbulence
379	c -----------------------------------------
380
381	IF(lecrit) THEN
382	PRINT*
383	PRINT*,'Diagnostic for the vertical turbulent mixing'
384	PRINT*,'Cd for momentum and potential temperature'
385
386	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
387	PRINT*,'Mixing coefficient for momentum and pot.temp.'
388	DO ilev=1,nlay
389	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
390	ENDDO
391	ENDIF
392
393
394
395
396	c-----------------------------------------------------------------------
397	c 4. inversion pour l'implicite sur u
398	c --------------------------------
399
400	c ** l'equation est
401	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
402	c avec
403	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
404	c et
405	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
406	c donc les entrees sont /zcdv/ pour la condition a la limite sol
407	c et /zkv/ = Ku
408
409	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
410	CALL multipl(ngrid,zcdv,zb0,zb)
411
412	DO ig=1,ngrid
413	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
414	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
415	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
416	ENDDO
417
418	DO ilay=nlay-1,1,-1
419	DO ig=1,ngrid
420	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
421	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
422	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
423	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
424	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
425	ENDDO
426	ENDDO
427
428	DO ig=1,ngrid
429	zu(ig,1)=zc(ig,1)
430	ENDDO
431	DO ilay=2,nlay
432	DO ig=1,ngrid
433	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
434	ENDDO
435	ENDDO
436
437
438
439
440
441	c-----------------------------------------------------------------------
442	c 5. inversion pour l'implicite sur v
443	c --------------------------------
444
445	c ** l'equation est
446	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
447	c avec
448	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
449	c et
450	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
451	c donc les entrees sont /zcdv/ pour la condition a la limite sol
452	c et /zkv/ = Kv
453
454	DO ig=1,ngrid
455	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
456	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
457	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
458	ENDDO
459
460	DO ilay=nlay-1,1,-1
461	DO ig=1,ngrid
462	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
463	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
464	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
465	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
466	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
467	ENDDO
468	ENDDO
469
470	DO ig=1,ngrid
471	zv(ig,1)=zc(ig,1)
472	ENDDO
473	DO ilay=2,nlay
474	DO ig=1,ngrid
475	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
476	ENDDO
477	ENDDO
478
479
480
481
482
483	c-----------------------------------------------------------------------
484	c 6. inversion pour l'implicite sur h sans oublier le couplage
485	c avec le sol (conduction)
486	c ------------------------
487
488	c ** l'equation est
489	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
490	c avec
491	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
492	c et
493	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
494	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
495	c et /zkh/ = Kh
496	c -------------
497
498	c Mass variation scheme:
499	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
500	CALL multipl(ngrid,zcdh,zb0,zb)
501
502	c on initialise dm c
503
504	pdtc(:,:)=0.
505	zt(:,:)=0.
506	dmice(:,:)=0.
507
508	c ** calcul de (d Planck / dT) a la temperature d'interface
509	c ------------------------------------------------------
510
511	z4st=4.5.67e-8ptimestep
512	DO ig=1,ngrid
513	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
514	ENDDO
515
516
517	! calcul de zc et zd pour la couche top en prenant en compte le terme
518	! de variation de masse (on fait une boucle pour que ça converge)
519
520	! Identification des points de grilles qui ont besoin de la correction
521
522	llnt(:)=1
523
524	DO ig=1,ngrid
525	DO l=1,nlay
526	if(zh(ig,l) .lt. zhcond(ig,l)) then
527	llnt(ig)=300
528	! 200 and 100 do not go beyond month 9 with normal dissipation
529	goto 5
530	endif
531	ENDDO
532	5 continue
533	ENDDO
534
535	DO ig=1,ngrid
536
537	! Initialization of z1 and zd, which do not depend on dmice
538
539	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
540	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
541
542	DO ilay=nlay-1,1,-1
543	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
544	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
545	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
546	ENDDO
547
548	! Convergence loop
549
550	DO j=1,llnt(ig)
551
552	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
553	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
554	& -betam(ig,nlay)*dmice(ig,nlay)
555	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
556	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
557
558	! calcul de zc et zd pour les couches du haut vers le bas
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	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
564	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
565	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
566	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
567	ENDDO
568
569	c ** calcul de la temperature_d'interface et de sa tendance.
570	c on ecrit que la somme des flux est nulle a l'interface
571	c a t + \delta t,
572	c c'est a dire le flux radiatif a {t + \delta t}
573	c + le flux turbulent a {t + \delta t}
574	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
575	c (notation K dt = /cpp*b/)
576	c + le flux dans le sol a t
577	c + l'evolution du flux dans le sol lorsque la temperature d'interface
578	c passe de sa valeur a t a sa valeur a {t + \delta t}.
579	c ----------------------------------------------------
580
581	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
582	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
583	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
584	ztsrf2(ig)=z1(ig)/z2(ig)
585	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
586
587	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
588
589	c ** et a partir de la temperature au sol on remonte
590	c -----------------------------------------------
591
592	DO ilay=2,nlay
593	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
594	ENDDO
595	DO ilay=1,nlay
596	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
597	ENDDO
598
599	c Condensation/sublimation in the atmosphere
600	c ------------------------------------------
601	c (computation of zcondicea and dmice)
602
603	zcondicea(ig,:)=0.
604	pdtc(ig,:)=0.
605
606	DO l=nlay , 1, -1
607	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
608	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
609	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
610	& ccondpdtc(ig,l)
611	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
612	END IF
613	ENDDO
614
615	ENDDO !of Do j=1,XXX
616
617	ENDDO !of Do ig=1,nlay
618
619	pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
620
621
622	c-----------------------------------------------------------------------
623	c TRACERS
624	c -------
625
626	if(tracer) then
627
628	c Using the wind modified by friction for lifting and sublimation
629	c ----------------------------------------------------------------
630	DO ig=1,ngrid
631	zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
632	zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
633	zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
634	ENDDO
635
636	c Calcul du flux vertical au bas de la premiere couche (dust) :
637	c -----------------------------------------------------------
638	do ig=1,ngridmx
639	rho(ig) = zb0(ig,1) /ptimestep
640	c zb(ig,1) = 0.
641	end do
642	c Dust lifting:
643	if (lifting) then
644	#ifndef MESOSCALE
645	if (doubleq.AND.submicron) then
646	do ig=1,ngrid
647	c if(co2ice(ig).lt.1) then
648	pdqsdif(ig,igcm_dust_mass) =
649	& -alpha_lift(igcm_dust_mass)
650	pdqsdif(ig,igcm_dust_number) =
651	& -alpha_lift(igcm_dust_number)
652	pdqsdif(ig,igcm_dust_submicron) =
653	& -alpha_lift(igcm_dust_submicron)
654	c end if
655	end do
656	else if (doubleq) then
657	do ig=1,ngrid
658	!!! soulevement constant
659	pdqsdif(ig,igcm_dust_mass) =
660	& -alpha_lift(igcm_dust_mass)
661	pdqsdif(ig,igcm_dust_number) =
662	& -alpha_lift(igcm_dust_number)
663	end do
664	else if (submicron) then
665	do ig=1,ngrid
666	pdqsdif(ig,igcm_dust_submicron) =
667	& -alpha_lift(igcm_dust_submicron)
668	end do
669	else
670	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
671	& pdqsdif)
672	endif !doubleq.AND.submicron
673	#else
674	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
675	& pdqsdif)
676	#endif
677	else
678	pdqsdif(1:ngrid,1:nq) = 0.
679	end if
680
681	c OU calcul de la valeur de q a la surface (water) :
682	c ----------------------------------------
683	if (water) then
684	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
685	end if
686
687	c Inversion pour l'implicite sur q
688	c --------------------------------
689	do iq=1,nq
690	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
691
692	if ((water).and.(iq.eq.igcm_h2o_vap)) then
693	c This line is required to account for turbulent transport
694	c from surface (e.g. ice) to mid-layer of atmosphere:
695	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
696	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
697	else ! (re)-initialize zb(:,1)
698	zb(1:ngrid,1)=0
699	end if
700
701	DO ig=1,ngrid
702	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
703	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
704	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
705	ENDDO
706
707	DO ilay=nlay-1,2,-1
708	DO ig=1,ngrid
709	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
710	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
711	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
712	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
713	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
714	ENDDO
715	ENDDO
716
717	if (water.and.(iq.eq.igcm_h2o_ice)) then
718	! special case for water ice tracer: do not include
719	! h2o ice tracer from surface (which is set when handling
720	! h2o vapour case (see further down).
721	DO ig=1,ngrid
722	z1(ig)=1./(za(ig,1)+zb(ig,1)+
723	$ zb(ig,2)*(1.-zd(ig,2)))
724	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
725	$ zb(ig,2)zc(ig,2)) z1(ig)
726	ENDDO
727	else ! general case
728	DO ig=1,ngrid
729	z1(ig)=1./(za(ig,1)+zb(ig,1)+
730	$ zb(ig,2)*(1.-zd(ig,2)))
731	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
732	$ zb(ig,2)*zc(ig,2) +
733	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
734	ENDDO
735	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
736
737	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
738	c Calculation for turbulent exchange with the surface (for ice)
739	DO ig=1,ngrid
740	zd(ig,1)=zb(ig,1)*z1(ig)
741	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
742
743	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
744	& *(zq1temp(ig)-qsat(ig))
745	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
746	END DO
747
748	DO ig=1,ngrid
749	if(.not.watercaptag(ig)) then
750	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
751	& .gt.pqsurf(ig,igcm_h2o_ice)) then
752	c write(,)'on sublime plus que qsurf!'
753	pdqsdif(ig,igcm_h2o_ice)=
754	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
755	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
756	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
757	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
758	$ zb(ig,2)*zc(ig,2) +
759	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
760	zq1temp(ig)=zc(ig,1)
761	endif
762	endif ! if (.not.watercaptag(ig))
763	c Starting upward calculations for water :
764	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
765	ENDDO ! of DO ig=1,ngrid
766	ELSE
767	c Starting upward calculations for simple mixing of tracer (dust)
768	DO ig=1,ngrid
769	zq(ig,1,iq)=zc(ig,1)
770	ENDDO
771	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
772
773	DO ilay=2,nlay
774	DO ig=1,ngrid
775	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
776	ENDDO
777	ENDDO
778	enddo ! of do iq=1,nq
779	end if ! of if(tracer)
780
781	c-----------------------------------------------------------------------
782	c 8. calcul final des tendances de la diffusion verticale
783	c ----------------------------------------------------
784
785	DO ilev = 1, nlay
786	DO ig=1,ngrid
787	pdudif(ig,ilev)=( zu(ig,ilev)-
788	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
789	pdvdif(ig,ilev)=( zv(ig,ilev)-
790	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
791	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
792	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
793	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
794	ENDDO
795	ENDDO
796
797	if (tracer) then
798	DO iq = 1, nq
799	DO ilev = 1, nlay
800	DO ig=1,ngrid
801	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
802	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
803	ENDDO
804	ENDDO
805	ENDDO
806	end if
807
808	c ** diagnostique final
809	c ------------------
810
811	IF(lecrit) THEN
812	PRINT*,'In vdif'
813	PRINT*,'Ts (t) and Ts (t+st)'
814	WRITE(*,'(a10,3a15)')
815	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
816	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
817	DO ilev=1,nlay
818	WRITE(*,'(4f15.7)')
819	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
820	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
821
822	ENDDO
823	ENDIF
824
825	RETURN
826	END

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

Download in other formats:

Original Format