Context Navigation

vdifc.F @ 1648

Last change on this file since 1648 was 1647, checked in by jvatant, 9 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File size: 15.2 KB

Line
1	subroutine vdifc(ngrid,nlay,nq,ppopsk,
2	& ptimestep,pcapcal,lecrit,
3	& pplay,pplev,pzlay,pzlev,pz0,
4	& pu,pv,ph,pq,ptsrf,pemis,pqsurf,
5	& pdhfi,pdqfi,pfluxsrf,
6	& pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
7	& pdqdif,pdqsdif,lastcall)
8
9	use radcommon_h, only : sigma
10	USE surfdat_h
11	USE tracer_h
12	use comcstfi_mod, only: g, r, cpp, rcp
13	use callkeys_mod, only: tracer,nosurf
14
15	implicit none
16
17	!==================================================================
18	!
19	! Purpose
20	! -------
21	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
22	!
23	! Implicit scheme
24	! We start by adding to variables x the physical tendencies
25	! already computed. We resolve the equation:
26	!
27	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
28	!
29	! Authors
30	! -------
31	! F. Hourdin, F. Forget, R. Fournier (199X)
32	! R. Wordsworth, B. Charnay (2010)
33	!
34	!==================================================================
35
36	!-----------------------------------------------------------------------
37	! declarations
38	! ------------
39
40
41	! arguments
42	! ---------
43	INTEGER ngrid,nlay
44	REAL ptimestep
45	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
46	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
47	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
48	REAL ptsrf(ngrid),pemis(ngrid)
49	REAL pdhfi(ngrid,nlay)
50	REAL pfluxsrf(ngrid)
51	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
52	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
53	REAL pq2(ngrid,nlay+1)
54
55
56	! Arguments added for condensation
57	REAL ppopsk(ngrid,nlay)
58	logical lecrit
59	REAL pz0
60
61	! Tracers
62	! --------
63	integer nq
64	real pqsurf(ngrid,nq)
65	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
66	real pdqdif(ngrid,nlay,nq)
67	real pdqsdif(ngrid,nq)
68
69	! local
70	! -----
71	integer ilev,ig,ilay,nlev
72
73	REAL z4st,zdplanck(ngrid)
74	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
75	REAL zcdv(ngrid),zcdh(ngrid)
76	REAL zcdv_true(ngrid),zcdh_true(ngrid)
77	REAL zu(ngrid,nlay),zv(ngrid,nlay)
78	REAL zh(ngrid,nlay)
79	REAL ztsrf2(ngrid)
80	REAL z1(ngrid),z2(ngrid)
81	REAL za(ngrid,nlay),zb(ngrid,nlay)
82	REAL zb0(ngrid,nlay)
83	REAL zc(ngrid,nlay),zd(ngrid,nlay)
84	REAL zcst1
85	REAL zu2!, a
86	REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
87	REAL evap(ngrid)
88	REAL zcq0(ngrid),zdq0(ngrid)
89	REAL zx_alf1(ngrid),zx_alf2(ngrid)
90
91	LOGICAL firstcall
92	SAVE firstcall
93	!$OMP THREADPRIVATE(firstcall)
94
95	LOGICAL lastcall
96
97	! variables added for CO2 condensation
98	! ------------------------------------
99	REAL hh !, zhcond(ngrid,nlay)
100	! REAL latcond,tcond1mb
101	! REAL acond,bcond
102	! SAVE acond,bcond
103	!!$OMP THREADPRIVATE(acond,bcond)
104	! DATA latcond,tcond1mb/5.9e5,136.27/
105
106	! Tracers
107	! -------
108	INTEGER iq
109	REAL zq(ngrid,nlay,nq)
110	REAL zq1temp(ngrid)
111	REAL rho(ngrid) ! near-surface air density
112	REAL qsat(ngrid)
113	DATA firstcall/.true./
114	REAL kmixmin
115
116	real, parameter :: karman=0.4
117	real cd0, roughratio
118
119	real masse, Wtot, Wdiff
120
121	real dqsdif_total(ngrid)
122	real zq0(ngrid)
123
124
125	! Coherence test
126	! --------------
127
128	IF (firstcall) THEN
129	! To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
130	! bcond=1./tcond1mb
131	! acond=r/latcond
132	! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
133	! PRINT*,' acond,bcond',acond,bcond
134
135	firstcall=.false.
136	ENDIF
137
138	!-----------------------------------------------------------------------
139	! 1. Initialisation
140	! -----------------
141
142	nlev=nlay+1
143
144	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
145	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
146	! ---------------------------------------------
147
148	DO ilay=1,nlay
149	DO ig=1,ngrid
150	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
151	ENDDO
152	ENDDO
153
154	zcst1=4.gptimestep/(R*R)
155	DO ilev=2,nlev-1
156	DO ig=1,ngrid
157	zb0(ig,ilev)=pplev(ig,ilev)*
158	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
159	s (ph(ig,ilev-1)+ph(ig,ilev))
160	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
161	s (pplay(ig,ilev-1)-pplay(ig,ilev))
162	ENDDO
163	ENDDO
164	DO ig=1,ngrid
165	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
166	ENDDO
167
168	dqsdif_total(:)=0.0
169
170	!-----------------------------------------------------------------------
171	! 2. Add the physical tendencies computed so far
172	! ----------------------------------------------
173
174	DO ilev=1,nlay
175	DO ig=1,ngrid
176	zu(ig,ilev)=pu(ig,ilev)
177	zv(ig,ilev)=pv(ig,ilev)
178	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
179	ENDDO
180	ENDDO
181	if(tracer) then
182	DO iq =1, nq
183	DO ilev=1,nlay
184	DO ig=1,ngrid
185	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
186	& pdqfi(ig,ilev,iq)*ptimestep
187	ENDDO
188	ENDDO
189	ENDDO
190	end if
191
192	!-----------------------------------------------------------------------
193	! 3. Turbulence scheme
194	! --------------------
195	!
196	! Source of turbulent kinetic energy at the surface
197	! -------------------------------------------------
198	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
199
200	DO ig=1,ngrid
201	roughratio = 1.E+0 + pzlay(ig,1)/pz0
202	cd0 = karman/log(roughratio)
203	cd0 = cd0*cd0
204	zcdv_true(ig) = cd0
205	zcdh_true(ig) = cd0
206	if (nosurf) then
207	zcdv_true(ig) = 0. !! disable sensible momentum flux
208	zcdh_true(ig) = 0. !! disable sensible heat flux
209	endif
210	ENDDO
211
212	DO ig=1,ngrid
213	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
214	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
215	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
216	ENDDO
217
218	! Turbulent diffusion coefficients in the boundary layer
219	! ------------------------------------------------------
220
221	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
222	& ,pu,pv,ph,zcdv_true
223	& ,pq2,zkv,zkh)
224
225	! Adding eddy mixing to mimic 3D general circulation in 1D
226	! R. Wordsworth & F. Forget (2010)
227	if ((ngrid.eq.1)) then
228	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
229	do ilev=1,nlay
230	do ig=1,ngrid
231	!zkh(ig,ilev) = 1.0
232	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
233	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
234	end do
235	end do
236	end if
237
238	!-----------------------------------------------------------------------
239	! 4. Implicit inversion of u
240	! --------------------------
241
242	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
243	! avec
244	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
245	! et
246	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
247	! donc les entrees sont /zcdv/ pour la condition a la limite sol
248	! et /zkv/ = Ku
249
250	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
251	CALL multipl(ngrid,zcdv,zb0,zb)
252
253	DO ig=1,ngrid
254	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
255	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
256	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
257	ENDDO
258
259	DO ilay=nlay-1,1,-1
260	DO ig=1,ngrid
261	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
262	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
263	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
264	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
265	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
266	ENDDO
267	ENDDO
268
269	DO ig=1,ngrid
270	zu(ig,1)=zc(ig,1)
271	ENDDO
272	DO ilay=2,nlay
273	DO ig=1,ngrid
274	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
275	ENDDO
276	ENDDO
277
278	!-----------------------------------------------------------------------
279	! 5. Implicit inversion of v
280	! --------------------------
281
282	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
283	! avec
284	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
285	! et
286	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
287	! donc les entrees sont /zcdv/ pour la condition a la limite sol
288	! et /zkv/ = Kv
289
290	DO ig=1,ngrid
291	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
292	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
293	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
294	ENDDO
295
296	DO ilay=nlay-1,1,-1
297	DO ig=1,ngrid
298	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
299	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
300	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
301	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
302	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
303	ENDDO
304	ENDDO
305
306	DO ig=1,ngrid
307	zv(ig,1)=zc(ig,1)
308	ENDDO
309	DO ilay=2,nlay
310	DO ig=1,ngrid
311	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
312	ENDDO
313	ENDDO
314
315	!----------------------------------------------------------------------------
316	! 6. Implicit inversion of h, not forgetting the coupling with the ground
317
318	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
319	! avec
320	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
321	! et
322	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
323	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
324	! et /zkh/ = Kh
325
326	! Using the wind modified by friction for lifting and sublimation
327	! ---------------------------------------------------------------
328	DO ig=1,ngrid
329	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
330	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
331	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
332	ENDDO
333
334	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
335	CALL multipl(ngrid,zcdh,zb0,zb)
336
337	DO ig=1,ngrid
338	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
339	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
340	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
341	ENDDO
342
343	DO ilay=nlay-1,2,-1
344	DO ig=1,ngrid
345	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
346	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
347	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
348	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
349	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
350	ENDDO
351	ENDDO
352
353	DO ig=1,ngrid
354	z1(ig)=1./(za(ig,1)+zb(ig,1)+
355	& zb(ig,2)*(1.-zd(ig,2)))
356	zc(ig,1)=(za(ig,1)*zh(ig,1)+
357	& zb(ig,2)zc(ig,2))z1(ig)
358	zd(ig,1)=zb(ig,1)*z1(ig)
359	ENDDO
360
361	! Calculate (d Planck / dT) at the interface temperature
362	! ------------------------------------------------------
363
364	z4st=4.0sigmaptimestep
365	DO ig=1,ngrid
366	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
367	ENDDO
368
369	! Calculate temperature tendency at the interface (dry case)
370	! ----------------------------------------------------------
371	! Sum of fluxes at interface at time t + \delta t gives change in T:
372	! radiative fluxes
373	! turbulent convective (sensible) heat flux
374	! flux (if any) from subsurface
375
376
377	DO ig=1,ngrid
378
379	z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
380	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
381	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
382	& +zdplanck(ig)
383	ztsrf2(ig) = z1(ig) / z2(ig)
384	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
385	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
386	ENDDO
387
388	! Recalculate temperature to top of atmosphere, starting from ground
389	! ------------------------------------------------------------------
390
391	DO ilay=2,nlay
392	DO ig=1,ngrid
393	hh = zh(ig,ilay-1)
394	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
395	ENDDO
396	ENDDO
397
398
399	!-----------------------------------------------------------------------
400	! TRACERS (no vapour)
401	! -------
402
403	if(tracer) then
404
405	! Calculate vertical flux from the bottom to the first layer (dust)
406	! -----------------------------------------------------------------
407	do ig=1,ngrid
408	rho(ig) = zb0(ig,1) /ptimestep
409	end do
410
411	call zerophys(ngrid*nq,pdqsdif)
412
413	! Implicit inversion of q
414	! -----------------------
415	do iq=1,nq
416
417	DO ig=1,ngrid
418	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
419	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
420	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
421	ENDDO
422
423	DO ilay=nlay-1,2,-1
424	DO ig=1,ngrid
425	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
426	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
427	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
428	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
429	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
430	ENDDO
431	ENDDO
432
433
434
435	DO ig=1,ngrid
436	z1(ig)=1./(za(ig,1)+
437	& zb(ig,2)*(1.-zdq(ig,2)))
438	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
439	& zb(ig,2)*zcq(ig,2)
440	& +(-pdqsdif(ig,iq))ptimestep)z1(ig)
441	! tracer flux from surface
442	! currently pdqsdif always zero here,
443	! so last line is superfluous
444	enddo
445
446
447	! Starting upward calculations for simple tracer mixing (e.g., dust)
448	do ig=1,ngrid
449	zq(ig,1,iq)=zcq(ig,1)
450	end do
451
452	do ilay=2,nlay
453	do ig=1,ngrid
454	zq(ig,ilay,iq)=zcq(ig,ilay)+
455	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
456	end do
457	end do
458
459
460	end do ! of do iq=1,nq
461	endif ! traceur
462
463
464	!-----------------------------------------------------------------------
465	! 8. Final calculation of the vertical diffusion tendencies
466	! -----------------------------------------------------------------
467
468	do ilev = 1, nlay
469	do ig=1,ngrid
470	pdudif(ig,ilev)=(zu(ig,ilev)-
471	& (pu(ig,ilev)))/ptimestep
472	pdvdif(ig,ilev)=(zv(ig,ilev)-
473	& (pv(ig,ilev)))/ptimestep
474	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
475
476	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
477	enddo
478	enddo
479
480	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
481	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
482	ENDDO
483
484	if (tracer) then
485	do iq = 1, nq
486	do ilev = 1, nlay
487	do ig=1,ngrid
488	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
489	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
490	& ptimestep
491	enddo
492	enddo
493	enddo
494
495	endif
496
497	! if(lastcall)then
498	! if(ngrid.eq.1)then
499	! print*,'Saving k.out...'
500	! OPEN(12,file='k.out',form='formatted')
501	! DO ilay=1,nlay
502	! write(12,*) zkh(1,ilay), pplay(1,ilay)
503	! ENDDO
504	! CLOSE(12)
505	! endif
506	! endif
507
508
509	return
510	end

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

Download in other formats:

Original Format