Context Navigation

vdifc.F @ 1242

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

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

Download in other formats:

Original Format