Context Navigation

vdifc.F @ 804

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

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

Download in other formats:

Original Format