Context Navigation

vdifc.F @ 1292

Last change on this file since 1292 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format