Context Navigation

vdifc.F @ 1233

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

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

Download in other formats:

Original Format