Context Navigation

vdifc.F @ 706

Last change on this file since 706 was 669, checked in by tnavarro, 13 years ago

Ice effective radius for output is weighted by ice total surface instead of ice total mass in physiq.F

New distribution for permanent ice reservoirs in surfini.F

icelocationmode = 1 allows for realistic ice distribution, no matter what resolution.
icelocationmode = 2 is predefined 32x24 or 64x48 and should be USED BY DEFAULT. It currently overestimates ice.
icelocationmode = 3 is good for quick logical definitions (paleoclimates,stability studies,etc...)

Possibility to change perihelion date in solar longitude as well as in sol in tabfi.F for newstarts
Improved latent heat release when ground ice sublimates: now a implicit scheme for ground temperature

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 T1,T2
116	SAVE T1,T2
117	DATA T1,T2/-877.5,807.5/ ! zeros of latent heat equation for ice
118
119	c Tracers :
120	c ~~~~~~~
121	INTEGER iq
122	REAL zq(ngridmx,nlayermx,nqmx)
123	REAL zq1temp(ngridmx)
124	REAL rho(ngridmx) ! near surface air density
125	REAL qsat(ngridmx)
126	DATA firstcall/.true./
127
128	REAL kmixmin
129
130	c Mass-variation scheme :
131	c ~~~~~~~
132
133	INTEGER j,l
134	REAL zcondicea(ngrid,nlayermx)
135	REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1)
136	REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx)
137	REAL pdtc(ngrid,nlayermx)
138	REAL zhs(ngridmx,nlayermx)
139	REAL ccond
140	SAVE ccond
141
142	c Theta_m formulation for mass-variation scheme :
143	c ~~~~~~~
144
145	INTEGER ico2,llnt(ngridmx)
146	SAVE ico2
147	REAL m_co2, m_noco2, A , B
148	SAVE A, B, m_co2, m_noco2
149	REAL vmr_co2(ngridmx,nlayermx)
150	REAL qco2,mmean
151
152	REAL sensibFlux(ngridmx)
153
154	#ifdef MESOSCALE
155	LOGICAL flag_LES !! pour LES avec isfflx!=0
156	#endif
157	c ** un petit test de coherence
158	c --------------------------
159
160	IF (firstcall) THEN
161	IF(ngrid.NE.ngridmx) THEN
162	PRINT*,'STOP dans vdifc'
163	PRINT*,'probleme de dimensions :'
164	PRINT*,'ngrid =',ngrid
165	PRINT*,'ngridmx =',ngridmx
166	STOP
167	ENDIF
168	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
169	bcond=1./tcond1mb
170	acond=r/latcond
171	ccond=cpp/(g*latcond)
172	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
173	PRINT*,' acond,bcond,ccond',acond,bcond,ccond
174
175
176	ico2=0
177
178	if (tracer) then
179	c Prepare Special treatment if one of the tracer is CO2 gas
180	do iq=1,nqmx
181	if (noms(iq).eq."co2") then
182	ico2=iq
183	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
184	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
185	c Compute A and B coefficient use to compute
186	c mean molecular mass Mair defined by
187	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
188	c 1/Mair = A*q(ico2) + B
189	A =(1/m_co2 - 1/m_noco2)
190	B=1/m_noco2
191	endif
192	enddo
193	end if
194
195	firstcall=.false.
196	ENDIF
197
198
199
200
201
202	c-----------------------------------------------------------------------
203	c 1. initialisation
204	c -----------------
205
206	nlev=nlay+1
207
208	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
209	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
210	c ----------------------------------------
211
212	DO ilay=1,nlay
213	DO ig=1,ngrid
214	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
215	! Mass variation scheme:
216	betam(ig,ilay)=-za(ig,ilay)latcond/(cppppopsk(ig,ilay))
217	ENDDO
218	ENDDO
219
220	zcst1=4.gptimestep/(r*r)
221	DO ilev=2,nlev-1
222	DO ig=1,ngrid
223	zb0(ig,ilev)=pplev(ig,ilev)*
224	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
225	s (ph(ig,ilev-1)+ph(ig,ilev))
226	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
227	s (pplay(ig,ilev-1)-pplay(ig,ilev))
228	ENDDO
229	ENDDO
230	DO ig=1,ngrid
231	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf(ig))
232	ENDDO
233
234	c ** diagnostique pour l'initialisation
235	c ----------------------------------
236
237	IF(lecrit) THEN
238	ig=ngrid/2+1
239	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
240	DO ilay=1,nlay
241	WRITE(*,'(6f11.5)')
242	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
243	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
244	ENDDO
245	PRINT*,'Pression (mbar) ,altitude (km),zb'
246	DO ilev=1,nlay
247	WRITE(*,'(3f15.7)')
248	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
249	s zb0(ig,ilev)
250	ENDDO
251	ENDIF
252
253	c -----------------------------------
254	c Potential Condensation temperature:
255	c -----------------------------------
256
257	c Compute CO2 Volume mixing ratio
258	c -------------------------------
259	if (callcond.and.(ico2.ne.0)) then
260	DO ilev=1,nlay
261	DO ig=1,ngrid
262	qco2=MAX(1.E-30
263	& ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
264	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
265	mmean=1/(A*qco2 +B)
266	vmr_co2(ig,ilev) = qco2*mmean/m_co2
267	ENDDO
268	ENDDO
269	else
270	DO ilev=1,nlay
271	DO ig=1,ngrid
272	vmr_co2(ig,ilev)=0.95
273	ENDDO
274	ENDDO
275	end if
276
277	c forecast of atmospheric temperature zt and frost temperature ztcond
278	c --------------------------------------------------------------------
279
280	DO ilev=1,nlay
281	DO ig=1,ngrid
282	ztcond(ig,ilev)=
283	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
284	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
285	ENDDO
286	ENDDO
287
288	ztcond(:,nlay+1)=ztcond(:,nlay)
289
290	if (callcond) then
291	DO ilev=1,nlay
292	DO ig=1,ngrid
293	! zhcond(ig,ilev) =
294	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
295	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
296	END DO
297	END DO
298	else
299	call zerophys(ngrid*nlay,zhcond)
300	end if
301
302
303	c-----------------------------------------------------------------------
304	c 2. ajout des tendances physiques
305	c -----------------------------
306
307	DO ilev=1,nlay
308	DO ig=1,ngrid
309	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
310	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
311	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
312	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
313	ENDDO
314	ENDDO
315	if(tracer) then
316	DO iq =1, nq
317	DO ilev=1,nlay
318	DO ig=1,ngrid
319	zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
320	ENDDO
321	ENDDO
322	ENDDO
323	end if
324
325	c-----------------------------------------------------------------------
326	c 3. schema de turbulence
327	c --------------------
328
329	c ** source d'energie cinetique turbulente a la surface
330	c (condition aux limites du schema de diffusion turbulente
331	c dans la couche limite
332	c ---------------------
333
334	CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
335	& ,zcdv_true,zcdh_true)
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	DO ig=1,ngrid
560	DO l=1,nlay
561	if(zh(ig,l) .lt. zhcond(ig,l)) then
562	llnt(ig)=300
563	! 200 and 100 do not go beyond month 9 with normal dissipation
564	goto 5
565	endif
566	ENDDO
567	5 continue
568	ENDDO
569
570	#ifdef MESOSCALE
571	ENDIF
572	#endif
573
574	DO ig=1,ngrid
575
576	! Initialization of z1 and zd, which do not depend on dmice
577
578	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
579	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
580
581	DO ilay=nlay-1,1,-1
582	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
583	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
584	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
585	ENDDO
586
587	! Convergence loop
588
589	DO j=1,llnt(ig)
590
591	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
592	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
593	& -betam(ig,nlay)*dmice(ig,nlay)
594	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
595	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
596
597	! calcul de zc et zd pour les couches du haut vers le bas
598
599	DO ilay=nlay-1,1,-1
600	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
601	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
602	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
603	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
604	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
605	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
606	ENDDO
607
608	c ** calcul de la temperature_d'interface et de sa tendance.
609	c on ecrit que la somme des flux est nulle a l'interface
610	c a t + \delta t,
611	c c'est a dire le flux radiatif a {t + \delta t}
612	c + le flux turbulent a {t + \delta t}
613	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
614	c (notation K dt = /cpp*b/)
615	c + le flux dans le sol a t
616	c + l'evolution du flux dans le sol lorsque la temperature d'interface
617	c passe de sa valeur a t a sa valeur a {t + \delta t}.
618	c ----------------------------------------------------
619
620	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzb(ig,1)*zc(ig,1)
621	s +zdplanck(ig)ptsrf(ig)+ pfluxsrf(ig)ptimestep
622	z2(ig)= pcapcal(ig)+cppzb(ig,1)(1.-zd(ig,1))+zdplanck(ig)
623	ztsrf2(ig)=z1(ig)/z2(ig)
624	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
625	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
626
627	c ** et a partir de la temperature au sol on remonte
628	c -----------------------------------------------
629
630	DO ilay=2,nlay
631	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
632	ENDDO
633	DO ilay=1,nlay
634	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
635	ENDDO
636
637	c Condensation/sublimation in the atmosphere
638	c ------------------------------------------
639	c (computation of zcondicea and dmice)
640
641	zcondicea(ig,:)=0.
642	pdtc(ig,:)=0.
643
644	DO l=nlay , 1, -1
645	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
646	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
647	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
648	& ccondpdtc(ig,l)
649	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
650	END IF
651	ENDDO
652
653	ENDDO !of Do j=1,XXX
654
655	ENDDO !of Do ig=1,nlay
656
657	pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
658
659	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
660	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zhs(ig,1)-ztsrf2(ig))
661	ENDDO
662
663	c-----------------------------------------------------------------------
664	c TRACERS
665	c -------
666
667	if(tracer) then
668
669	c Using the wind modified by friction for lifting and sublimation
670	c ----------------------------------------------------------------
671
672	! This is computed above and takes into account surface-atmosphere flux
673	! enhancement by subgrid gustiness and atmospheric-stability related
674	! variations of transfer coefficients.
675
676	! DO ig=1,ngrid
677	! zu2(ig)=zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
678	! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
679	! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
680	! ENDDO
681
682	c Calcul du flux vertical au bas de la premiere couche (dust) :
683	c -----------------------------------------------------------
684	do ig=1,ngridmx
685	rho(ig) = zb0(ig,1) /ptimestep
686	c zb(ig,1) = 0.
687	end do
688	c Dust lifting:
689	if (lifting) then
690	#ifndef MESOSCALE
691	if (doubleq.AND.submicron) then
692	do ig=1,ngrid
693	c if(co2ice(ig).lt.1) then
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	pdqsdif(ig,igcm_dust_submicron) =
699	& -alpha_lift(igcm_dust_submicron)
700	c end if
701	end do
702	else if (doubleq) then
703	do ig=1,ngrid
704	if(co2ice(ig).lt.1) then ! soulevement pas constant
705	pdqsdif(ig,igcm_dust_mass) =
706	& -alpha_lift(igcm_dust_mass)
707	pdqsdif(ig,igcm_dust_number) =
708	& -alpha_lift(igcm_dust_number)
709	end if
710	end do
711	else if (submicron) then
712	do ig=1,ngrid
713	pdqsdif(ig,igcm_dust_submicron) =
714	& -alpha_lift(igcm_dust_submicron)
715	end do
716	else
717	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
718	& pdqsdif)
719	endif !doubleq.AND.submicron
720	#else
721	call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
722	& pdqsdif)
723	#endif
724	else
725	pdqsdif(1:ngrid,1:nq) = 0.
726	end if
727
728	c OU calcul de la valeur de q a la surface (water) :
729	c ----------------------------------------
730	if (water) then
731	call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
732	end if
733
734	c Inversion pour l'implicite sur q
735	c --------------------------------
736	do iq=1,nq
737	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
738
739	if ((water).and.(iq.eq.igcm_h2o_vap)) then
740	c This line is required to account for turbulent transport
741	c from surface (e.g. ice) to mid-layer of atmosphere:
742	CALL multipl(ngrid,zcdv,zb0,zb(1,1))
743	CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
744	else ! (re)-initialize zb(:,1)
745	zb(1:ngrid,1)=0
746	end if
747
748	DO ig=1,ngrid
749	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
750	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
751	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
752	ENDDO
753
754	DO ilay=nlay-1,2,-1
755	DO ig=1,ngrid
756	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
757	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
758	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
759	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
760	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
761	ENDDO
762	ENDDO
763
764	if (water.and.(iq.eq.igcm_h2o_ice)) then
765	! special case for water ice tracer: do not include
766	! h2o ice tracer from surface (which is set when handling
767	! h2o vapour case (see further down).
768	DO ig=1,ngrid
769	z1(ig)=1./(za(ig,1)+zb(ig,1)+
770	$ zb(ig,2)*(1.-zd(ig,2)))
771	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
772	$ zb(ig,2)zc(ig,2)) z1(ig)
773	ENDDO
774	else ! general case
775	DO ig=1,ngrid
776	z1(ig)=1./(za(ig,1)+zb(ig,1)+
777	$ zb(ig,2)*(1.-zd(ig,2)))
778	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
779	$ zb(ig,2)*zc(ig,2) +
780	$ (-pdqsdif(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
781	ENDDO
782	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
783
784	IF ((water).and.(iq.eq.igcm_h2o_vap)) then
785	c Calculation for turbulent exchange with the surface (for ice)
786	DO ig=1,ngrid
787	zd(ig,1)=zb(ig,1)*z1(ig)
788	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
789
790	pdqsdif(ig,igcm_h2o_ice)=rho(ig)dryness(ig)zcdv(ig)
791	& *(zq1temp(ig)-qsat(ig))
792	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
793	END DO
794
795	DO ig=1,ngrid
796	if(.not.watercaptag(ig)) then
797	if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
798	& .gt.pqsurf(ig,igcm_h2o_ice)) then
799	c write(,)'on sublime plus que qsurf!'
800	pdqsdif(ig,igcm_h2o_ice)=
801	& -pqsurf(ig,igcm_h2o_ice)/ptimestep
802	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
803	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
804	zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
805	$ zb(ig,2)*zc(ig,2) +
806	$ (-pdqsdif(ig,igcm_h2o_ice)) ptimestep) z1(ig)
807	zq1temp(ig)=zc(ig,1)
808	endif
809	endif ! if (.not.watercaptag(ig))
810	c Starting upward calculations for water :
811	zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
812
813	c Take into account H2O latent heat in surface energy budget
814	c We solve dT/dt = (2834.3-0.28(T-To)-0.004(T-To)^2)1e3iceflux/cpp
815	tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
816
817	tsrf_lw(ig) = (T1+T2)*(T1+T2)
818	& - 4(T2T1 - (tsrf_lw(ig)-T1)(tsrf_lw(ig)-T2)
819	& exp( -0.25abs(T1-T2)pdqsdif(ig,igcm_h2o_ice)
820	& *ptimestep/pcapcal(ig)) )
821	tsrf_lw(ig) = (T1+T2)/2 + sqrt(tsrf_lw(ig))/2 ! surface temperature at t+1
822
823	pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
824
825	if(pqsurf(ig,igcm_h2o_ice)
826	& +pdqsdif(ig,igcm_h2o_ice)*ptimestep
827	& .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
828	& pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
829
830	ENDDO ! of DO ig=1,ngrid
831	ELSE
832	c Starting upward calculations for simple mixing of tracer (dust)
833	DO ig=1,ngrid
834	zq(ig,1,iq)=zc(ig,1)
835	ENDDO
836	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
837
838	DO ilay=2,nlay
839	DO ig=1,ngrid
840	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
841	ENDDO
842	ENDDO
843	enddo ! of do iq=1,nq
844	end if ! of if(tracer)
845
846
847	c-----------------------------------------------------------------------
848	c 8. calcul final des tendances de la diffusion verticale
849	c ----------------------------------------------------
850
851	DO ilev = 1, nlay
852	DO ig=1,ngrid
853	pdudif(ig,ilev)=( zu(ig,ilev)-
854	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
855	pdvdif(ig,ilev)=( zv(ig,ilev)-
856	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
857	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
858	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
859	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
860	ENDDO
861	ENDDO
862
863	if (tracer) then
864	DO iq = 1, nq
865	DO ilev = 1, nlay
866	DO ig=1,ngrid
867	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
868	$ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
869	ENDDO
870	ENDDO
871	ENDDO
872	end if
873
874	c ** diagnostique final
875	c ------------------
876
877	IF(lecrit) THEN
878	PRINT*,'In vdif'
879	PRINT*,'Ts (t) and Ts (t+st)'
880	WRITE(*,'(a10,3a15)')
881	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
882	PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
883	DO ilev=1,nlay
884	WRITE(*,'(4f15.7)')
885	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
886	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
887
888	ENDDO
889	ENDIF
890
891	RETURN
892	END

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

Download in other formats:

Original Format