Context Navigation

vdifc.F @ 629

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

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

Download in other formats:

Original Format