Context Navigation

vdifc.F @ 3529

Last change on this file since 3529 was 2743, checked in by aslmd, 2 years ago

made the vertical discretization for Titan LES similar to that of Mars (simple linear law for znw and refinement close to the surface). no more need to prepare a file levels, just put ztop and number of levels in namelist. a file levels is saved at the end of the znw computation so that it can be read in update_inputs_physiq_mod in the case of Titan.

added MESOSCALE precompiling flags in call_profilgases to set constant methane abundance. therefore no need for profile.def file and the number of vertical levels is set by namelist -- no need to set the same number in both levels and profile.def

outputing ustar from vdifc

ending hardcoding of ptop in update_inputs_physiq_mod for Titan and generic LES/mesoscale. streamlining grod%ptop through module_lmd_driver as ptopwrf within variables_mod module.

adding flux outputs to check surface and atmosphere energy budgets

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

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

Download in other formats:

Original Format