Context Navigation

vdifc.F @ 1009

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

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

Download in other formats:

Original Format