Context Navigation

vdifc.F @ 745

Last change on this file since 745 was 742, checked in by tnavarro, 13 years ago
latent heat for sublimating ice turned off: at high obliquity it gives horrible, frightful, nightmarish ground temperatures
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	DO ilev=1,nlay
282	DO ig=1,ngrid
283	ztcond(ig,ilev)=
284	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
285	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
286	ENDDO
287	ENDDO
288
289	ztcond(:,nlay+1)=ztcond(:,nlay)
290
291	if (callcond) then
292	DO ilev=1,nlay
293	DO ig=1,ngrid
294	! zhcond(ig,ilev) =
295	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
296	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
297	END DO
298	END DO
299	else
300	call zerophys(ngrid*nlay,zhcond)
301	end if
302
303
304	c-----------------------------------------------------------------------
305	c 2. ajout des tendances physiques
306	c -----------------------------
307
308	DO ilev=1,nlay
309	DO ig=1,ngrid
310	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
311	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
312	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
313	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
314	ENDDO
315	ENDDO
316	if(tracer) then
317	DO iq =1, nq
318	DO ilev=1,nlay
319	DO ig=1,ngrid
320	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
321	ENDDO
322	ENDDO
323	ENDDO
324	end if
325
326	c-----------------------------------------------------------------------
327	c 3. schema de turbulence
328	c --------------------
329
330	c ** source d'energie cinetique turbulente a la surface
331	c (condition aux limites du schema de diffusion turbulente
332	c dans la couche limite
333	c ---------------------
334
335	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
336	& ,zcdv_true,zcdh_true)
337
338	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
339
340	IF (callrichsl) THEN
341	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
342	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
343	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
344	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
345
346	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
347	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
348
349	tst(:)=0.
350	DO ig=1,ngrid
351	IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u=t=0
352	tst(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ust(ig)
353	ENDIF
354	ENDDO
355
356	ELSE
357	zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
358	zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
359	ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
360	tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
361	ENDIF
362
363	! Some usefull diagnostics for the new surface layer parametrization :
364
365	! call WRITEDIAGFI(ngridmx,'pcdv',
366	! & 'momentum drag','no units',
367	! & 2,zcdv_true)
368	! call WRITEDIAGFI(ngridmx,'pcdh',
369	! & 'heat drag','no units',
370	! & 2,zcdh_true)
371	! call WRITEDIAGFI(ngridmx,'ust',
372	! & 'friction velocity','m/s',2,ust)
373	! call WRITEDIAGFI(ngridmx,'tst',
374	! & 'friction temperature','K',2,tst)
375	! call WRITEDIAGFI(ngridmx,'rm-1',
376	! & 'aerodyn momentum conductance','m/s',
377	! & 2,zcdv)
378	! call WRITEDIAGFI(ngridmx,'rh-1',
379	! & 'aerodyn heat conductance','m/s',
380	! & 2,zcdh)
381
382	c ** schema de diffusion turbulente dans la couche limite
383	c ----------------------------------------------------
384	IF (.not. calltherm) THEN
385
386	CALL vdif_kc(ptimestep,g,pzlev,pzlay
387	& ,pu,pv,ph,zcdv_true
388	& ,pq2,zkv,zkh,zq)
389
390	ELSE
391
392	pt(:,:)=ph(:,:)*ppopsk(:,:)
393	CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
394	s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ust
395	s ,9)
396
397	ENDIF
398
399	if ((doubleq).and.(ngrid.eq.1)) then
400	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
401	do ilev=1,nlay
402	do ig=1,ngrid
403	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
404	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
405	end do
406	end do
407	end if
408
409	c ** diagnostique pour le schema de turbulence
410	c -----------------------------------------
411
412	IF(lecrit) THEN
413	PRINT*
414	PRINT*,'Diagnostic for the vertical turbulent mixing'
415	PRINT*,'Cd for momentum and potential temperature'
416
417	PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
418	PRINT*,'Mixing coefficient for momentum and pot.temp.'
419	DO ilev=1,nlay
420	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
421	ENDDO
422	ENDIF
423
424
425
426
427	c-----------------------------------------------------------------------
428	c 4. inversion pour l'implicite sur u
429	c --------------------------------
430
431	c ** l'equation est
432	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
433	c avec
434	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
435	c et
436	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
437	c donc les entrees sont /zcdv/ pour la condition a la limite sol
438	c et /zkv/ = Ku
439
440	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
441	CALL multipl(ngrid,zcdv,zb0,zb)
442
443	DO ig=1,ngrid
444	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
445	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
446	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
447	ENDDO
448
449	DO ilay=nlay-1,1,-1
450	DO ig=1,ngrid
451	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
452	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
453	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
454	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
455	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
456	ENDDO
457	ENDDO
458
459	DO ig=1,ngrid
460	zu(ig,1)=zc(ig,1)
461	ENDDO
462	DO ilay=2,nlay
463	DO ig=1,ngrid
464	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
465	ENDDO
466	ENDDO
467
468
469
470
471
472	c-----------------------------------------------------------------------
473	c 5. inversion pour l'implicite sur v
474	c --------------------------------
475
476	c ** l'equation est
477	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
478	c avec
479	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
480	c et
481	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
482	c donc les entrees sont /zcdv/ pour la condition a la limite sol
483	c et /zkv/ = Kv
484
485	DO ig=1,ngrid
486	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
487	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
488	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
489	ENDDO
490
491	DO ilay=nlay-1,1,-1
492	DO ig=1,ngrid
493	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
494	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
495	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
496	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
497	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
498	ENDDO
499	ENDDO
500
501	DO ig=1,ngrid
502	zv(ig,1)=zc(ig,1)
503	ENDDO
504	DO ilay=2,nlay
505	DO ig=1,ngrid
506	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
507	ENDDO
508	ENDDO
509
510
511
512
513
514	c-----------------------------------------------------------------------
515	c 6. inversion pour l'implicite sur h sans oublier le couplage
516	c avec le sol (conduction)
517	c ------------------------
518
519	c ** l'equation est
520	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
521	c avec
522	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
523	c et
524	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
525	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
526	c et /zkh/ = Kh
527	c -------------
528
529	c Mass variation scheme:
530	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
531	CALL multipl(ngrid,zcdh,zb0,zb)
532
533	c on initialise dm c
534
535	pdtc(:,:)=0.
536	zt(:,:)=0.
537	dmice(:,:)=0.
538
539	c ** calcul de (d Planck / dT) a la temperature d'interface
540	c ------------------------------------------------------
541
542	z4st=4.5.67e-8ptimestep
543	IF (tke_heat_flux .eq. 0.) THEN
544	DO ig=1,ngrid
545	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
546	ENDDO
547	ELSE
548	zdplanck(:)=0.
549	ENDIF
550
551	! calcul de zc et zd pour la couche top en prenant en compte le terme
552	! de variation de masse (on fait une boucle pour que ça converge)
553
554	! Identification des points de grilles qui ont besoin de la correction
555
556	llnt(:)=1
557	#ifdef MESOSCALE
558	IF (.not.flag_LES) THEN
559	#endif
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
571	#ifdef MESOSCALE
572	ENDIF
573	#endif
574
575	DO ig=1,ngrid
576
577	! Initialization of z1 and zd, which do not depend on dmice
578
579	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
580	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
581
582	DO ilay=nlay-1,1,-1
583	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
584	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
585	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
586	ENDDO
587
588	! Convergence loop
589
590	DO j=1,llnt(ig)
591
592	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
593	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
594	& -betam(ig,nlay)*dmice(ig,nlay)
595	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
596	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
597
598	! calcul de zc et zd pour les couches du haut vers le bas
599
600	DO ilay=nlay-1,1,-1
601	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
602	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
603	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
604	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
605	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
606	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
607	ENDDO
608
609	c ** calcul de la temperature_d'interface et de sa tendance.
610	c on ecrit que la somme des flux est nulle a l'interface
611	c a t + \delta t,
612	c c'est a dire le flux radiatif a {t + \delta t}
613	c + le flux turbulent a {t + \delta t}
614	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
615	c (notation K dt = /cpp*b/)
616	c + le flux dans le sol a t
617	c + l'evolution du flux dans le sol lorsque la temperature d'interface
618	c passe de sa valeur a t a sa valeur a {t + \delta t}.
619	c ----------------------------------------------------
620
621	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
622	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
623	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
624	ztsrf2(ig)=z1(ig)/z2(ig)
625	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
626	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
627
628	c ** et a partir de la temperature au sol on remonte
629	c -----------------------------------------------
630
631	DO ilay=2,nlay
632	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
633	ENDDO
634	DO ilay=1,nlay
635	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
636	ENDDO
637
638	c Condensation/sublimation in the atmosphere
639	c ------------------------------------------
640	c (computation of zcondicea and dmice)
641
642	zcondicea(ig,:)=0.
643	pdtc(ig,:)=0.
644
645	DO l=nlay , 1, -1
646	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
647	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
648	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
649	& ccondpdtc(ig,l)
650	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
651	END IF
652	ENDDO
653
654	ENDDO !of Do j=1,XXX
655
656	ENDDO !of Do ig=1,nlay
657
658	pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
659
660	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
661	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zhs(ig,1)-ztsrf2(ig))
662	ENDDO
663
664	c-----------------------------------------------------------------------
665	c TRACERS
666	c -------
667
668	if(tracer) then
669
670	c Using the wind modified by friction for lifting and sublimation
671	c ----------------------------------------------------------------
672
673	! This is computed above and takes into account surface-atmosphere flux
674	! enhancement by subgrid gustiness and atmospheric-stability related
675	! variations of transfer coefficients.
676
677	! DO ig=1,ngrid
678	! zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
679	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
680	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
681	! ENDDO
682
683	c Calcul du flux vertical au bas de la premiere couche (dust) :
684	c -----------------------------------------------------------
685	do ig=1,ngridmx
686	rho(ig) = zb0(ig,1) /ptimestep
687	c zb(ig,1) = 0.
688	end do
689	c Dust lifting:
690	if (lifting) then
691	#ifndef MESOSCALE
692	if (doubleq.AND.submicron) then
693	do ig=1,ngrid
694	c if(co2ice(ig).lt.1) then
695	pdqsdif(ig,igcm_dust_mass) =
696	& -alpha_lift(igcm_dust_mass)
697	pdqsdif(ig,igcm_dust_number) =
698	& -alpha_lift(igcm_dust_number)
699	pdqsdif(ig,igcm_dust_submicron) =
700	& -alpha_lift(igcm_dust_submicron)
701	c end if
702	end do
703	else if (doubleq) then
704	do ig=1,ngrid
705	if(co2ice(ig).lt.1) then ! soulevement pas constant
706	pdqsdif(ig,igcm_dust_mass) =
707	& -alpha_lift(igcm_dust_mass)
708	pdqsdif(ig,igcm_dust_number) =
709	& -alpha_lift(igcm_dust_number)
710	end if
711	end do
712	else if (submicron) then
713	do ig=1,ngrid
714	pdqsdif(ig,igcm_dust_submicron) =
715	& -alpha_lift(igcm_dust_submicron)
716	end do
717	else
718	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
719	& pdqsdif)
720	endif !doubleq.AND.submicron
721	#else
722	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
723	& pdqsdif)
724	#endif
725	else
726	pdqsdif(1:ngrid,1:nq) = 0.
727	end if
728
729	c OU calcul de la valeur de q a la surface (water) :
730	c ----------------------------------------
731	if (water) then
732	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
733	end if
734
735	c Inversion pour l'implicite sur q
736	c --------------------------------
737	do iq=1,nq
738	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
739
740	if ((water).and.(iq.eq.igcm_h2o_vap)) then
741	c This line is required to account for turbulent transport
742	c from surface (e.g. ice) to mid-layer of atmosphere:
743	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
744	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
745	else ! (re)-initialize zb(:,1)
746	zb(1:ngrid,1)=0
747	end if
748
749	DO ig=1,ngrid
750	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
751	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
752	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
753	ENDDO
754
755	DO ilay=nlay-1,2,-1
756	DO ig=1,ngrid
757	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
758	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
759	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
760	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
761	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
762	ENDDO
763	ENDDO
764
765	if (water.and.(iq.eq.igcm_h2o_ice)) then
766	! special case for water ice tracer: do not include
767	! h2o ice tracer from surface (which is set when handling
768	! h2o vapour case (see further down).
769	DO ig=1,ngrid
770	z1(ig)=1./(za(ig,1)+zb(ig,1)+
771	$ zb(ig,2)*(1.-zd(ig,2)))
772	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
773	$ zb(ig,2)zc(ig,2)) z1(ig)
774	ENDDO
775	else ! general case
776	DO ig=1,ngrid
777	z1(ig)=1./(za(ig,1)+zb(ig,1)+
778	$ zb(ig,2)*(1.-zd(ig,2)))
779	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
780	$ zb(ig,2)*zc(ig,2) +
781	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
782	ENDDO
783	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
784
785	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
786	c Calculation for turbulent exchange with the surface (for ice)
787	DO ig=1,ngrid
788	zd(ig,1)=zb(ig,1)*z1(ig)
789	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
790
791	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
792	& *(zq1temp(ig)-qsat(ig))
793	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
794	END DO
795
796	DO ig=1,ngrid
797	if(.not.watercaptag(ig)) then
798	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
799	& .gt.pqsurf(ig,igcm_h2o_ice)) then
800	c write(,)'on sublime plus que qsurf!'
801	pdqsdif(ig,igcm_h2o_ice)=
802	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
803	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
804	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
805	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
806	$ zb(ig,2)*zc(ig,2) +
807	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
808	zq1temp(ig)=zc(ig,1)
809	endif
810	endif ! if (.not.watercaptag(ig))
811	c Starting upward calculations for water :
812	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
813
814	!c Take into account H2O latent heat in surface energy budget
815	!c We solve dT/dt = (2834.3-0.28(T-To)-0.004(T-To)^2)1e3iceflux/cpp
816	! tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
817	!
818	! alpha = exp(-4abs(T1-T2)pdqsdif(ig,igcm_h2o_ice)
819	! & *ptimestep/pcapcal(ig))
820	!
821	! tsrf_lw(ig) = (tsrf_lw(ig)(T2-alphaT1)+T1T2(alpha-1))
822	! & /(tsrf_lw(ig)(1-alpha)+alphaT2-T1) ! surface temperature at t+1
823	!
824	! pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
825
826	if(pqsurf(ig,igcm_h2o_ice)
827	& +pdqsdif(ig,igcm_h2o_ice)*ptimestep
828	& .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
829	& pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
830
831	ENDDO ! of DO ig=1,ngrid
832	ELSE
833	c Starting upward calculations for simple mixing of tracer (dust)
834	DO ig=1,ngrid
835	zq(ig,1,iq)=zc(ig,1)
836	ENDDO
837	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
838
839	DO ilay=2,nlay
840	DO ig=1,ngrid
841	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
842	ENDDO
843	ENDDO
844	enddo ! of do iq=1,nq
845	end if ! of if(tracer)
846
847
848	c-----------------------------------------------------------------------
849	c 8. calcul final des tendances de la diffusion verticale
850	c ----------------------------------------------------
851
852	DO ilev = 1, nlay
853	DO ig=1,ngrid
854	pdudif(ig,ilev)=( zu(ig,ilev)-
855	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
856	pdvdif(ig,ilev)=( zv(ig,ilev)-
857	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
858	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
859	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
860	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
861	ENDDO
862	ENDDO
863
864	if (tracer) then
865	DO iq = 1, nq
866	DO ilev = 1, nlay
867	DO ig=1,ngrid
868	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
869	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
870	ENDDO
871	ENDDO
872	ENDDO
873	end if
874
875	c ** diagnostique final
876	c ------------------
877
878	IF(lecrit) THEN
879	PRINT*,'In vdif'
880	PRINT*,'Ts (t) and Ts (t+st)'
881	WRITE(*,'(a10,3a15)')
882	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
883	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
884	DO ilev=1,nlay
885	WRITE(*,'(4f15.7)')
886	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
887	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
888
889	ENDDO
890	ENDIF
891
892	RETURN
893	END

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

Download in other formats:

Original Format