Context Navigation

vdifc_mod.F @ 3558

Last change on this file since 3558 was 3275, checked in by afalco, 10 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
File size: 17.1 KB

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

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

Download in other formats:

Original Format