Context Navigation

turbdiff_mod.F90 @ 3586

Last change on this file since 3586 was 3228, checked in by tbertrand, 12 months ago
Pluto GCM: Cleaning code and adding extra options for newstart.F
File size: 18.3 KB

Line
1	module turbdiff_mod
2
3	implicit none
4
5	contains
6
7	subroutine turbdiff(ngrid,nlay,nq, &
8	ptimestep,pcapcal, &
9	pplay,pplev,pzlay,pzlev,pz0, &
10	pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf, &
11	pdtfi,pdqfi,pfluxsrf, &
12	Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, &
13	pdqdif,pdqevap,pdqsdif,flux_u,flux_v)
14
15	use radcommon_h, only : sigma, glat
16	use surfdat_h, only: dryness
17	use comcstfi_mod, only: rcp, g, r, cpp
18	use callkeys_mod, only: tracer,nosurf,kmixmin
19	use turb_mod, only : ustar
20	#ifdef MESOSCALE
21	use comm_wrf, only : comm_LATENT_HF
22	#endif
23	implicit none
24
25	!==================================================================
26	!
27	! Purpose
28	! -------
29	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
30	!
31	! Implicit scheme
32	! We start by adding to variables x the physical tendencies
33	! already computed. We resolve the equation:
34	!
35	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
36	!
37	! Authors
38	! -------
39	! F. Hourdin, F. Forget, R. Fournier (199X)
40	! R. Wordsworth, B. Charnay (2010)
41	! J. Leconte (2012): To f90
42	! - Rewritten the diffusion scheme to conserve total enthalpy
43	! by accounting for dissipation of turbulent kinetic energy.
44	! - Accounting for lost mean flow kinetic energy should come soon.
45	!
46	!==================================================================
47
48	!-----------------------------------------------------------------------
49	! declarations
50	! ------------
51
52	! arguments
53	! ---------
54	INTEGER,INTENT(IN) :: ngrid
55	INTEGER,INTENT(IN) :: nlay
56	REAL,INTENT(IN) :: ptimestep
57	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
58	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
59	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
60	REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
61	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
62	REAL,INTENT(IN) :: pemis(ngrid)
63	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
64	REAL,INTENT(IN) :: pfluxsrf(ngrid)
65	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
66	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
67	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
68	REAL,INTENT(OUT) :: sensibFlux(ngrid)
69	REAL,INTENT(IN) :: pcapcal(ngrid)
70	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
71	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
72
73	REAL,INTENT(IN) :: pz0
74
75	! Tracers
76	! --------
77	integer,intent(in) :: nq
78	real,intent(in) :: pqsurf(ngrid,nq)
79	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
80	real,intent(out) :: pdqdif(ngrid,nlay,nq)
81	real,intent(out) :: pdqsdif(ngrid,nq)
82
83	! local
84	! -----
85	integer ilev,ig,ilay,nlev
86
87	REAL z4st,zdplanck(ngrid)
88	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
89	REAL zcdv(ngrid),zcdh(ngrid)
90	REAL zcdv_true(ngrid),zcdh_true(ngrid)
91	REAL zu(ngrid,nlay),zv(ngrid,nlay)
92	REAL zh(ngrid,nlay),zt(ngrid,nlay)
93	REAL ztsrf(ngrid)
94	REAL z1(ngrid),z2(ngrid)
95	REAL zmass(ngrid,nlay)
96	REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
97	REAL zb0(ngrid,nlay)
98	REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
99	REAL zcv(ngrid,nlay),zdv(ngrid,nlay) !inversion coefficient for winds
100	REAL zct(ngrid,nlay),zdt(ngrid,nlay) !inversion coefficient for temperature
101	REAL zcq(ngrid,nlay),zdq(ngrid,nlay) !inversion coefficient for tracers
102	REAL zcst1
103	REAL zu2!, a
104	REAL zcq0(ngrid),zdq0(ngrid)
105	REAL zx_alf1(ngrid),zx_alf2(ngrid)
106	! 1D eddy diffusion coefficient
107	REAL kzz_eddy(nlay)
108	REAL pmin_kzz
109
110	LOGICAL,SAVE :: firstcall=.true.
111	!$OMP THREADPRIVATE(firstcall)
112
113	! Tracers
114	! -------
115	INTEGER iq
116	REAL zq(ngrid,nlay,nq)
117	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
118	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
119	REAL zdmassevap(ngrid)
120	REAL rho(ngrid) ! near-surface air density
121
122	! Variables added for implicit latent heat inclusion
123	! --------------------------------------------------
124	real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
125
126	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
127	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
128
129	real, parameter :: karman=0.4
130	real cd0, roughratio
131
132	real dqsdif_total(ngrid)
133	real zq0(ngrid)
134
135
136	! Coherence test
137	! --------------
138
139	IF (firstcall) THEN
140	ivap=1
141	iliq=0
142	iliq_surf=0
143	iice_surf=0 ! simply to make the code legible
144	sensibFlux(:)=0.
145	firstcall=.false.
146	ENDIF
147
148	!-----------------------------------------------------------------------
149	! 1. Initialisation
150	! -----------------
151
152	nlev=nlay+1
153
154	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
155	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
156	! ---------------------------------------------
157
158	DO ilay=1,nlay
159	DO ig=1,ngrid
160	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
161	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
162	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
163	ENDDO
164	ENDDO
165
166	zcst1=4.gptimestep/(R*R)
167	DO ilev=2,nlev-1
168	DO ig=1,ngrid
169	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
170	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
171	ENDDO
172	ENDDO
173	DO ig=1,ngrid
174	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
175	ENDDO
176	dqsdif_total(:)=0.0
177
178	!-----------------------------------------------------------------------
179	! 2. Add the physical tendencies computed so far
180	! ----------------------------------------------
181
182	DO ilev=1,nlay
183	DO ig=1,ngrid
184	zu(ig,ilev)=pu(ig,ilev)
185	zv(ig,ilev)=pv(ig,ilev)
186	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
187	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
188	ENDDO
189	ENDDO
190	if(tracer) then
191	DO iq =1, nq
192	DO ilev=1,nlay
193	DO ig=1,ngrid
194	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
195	ENDDO
196	ENDDO
197	ENDDO
198	end if
199
200	!-----------------------------------------------------------------------
201	! 3. Turbulence scheme
202	! --------------------
203	!
204	! Source of turbulent kinetic energy at the surface
205	! -------------------------------------------------
206	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
207
208	DO ig=1,ngrid
209	roughratio = 1. + pzlay(ig,1)/pz0
210	cd0 = karman/log(roughratio)
211	cd0 = cd0*cd0
212	zcdv_true(ig) = cd0
213	zcdh_true(ig) = cd0
214	if(nosurf)then
215	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
216	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
217	endif
218	ENDDO
219
220	DO ig=1,ngrid
221	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
222	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
223	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
224	ENDDO
225
226	! Turbulent diffusion coefficients in the boundary layer
227	! ------------------------------------------------------
228
229	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
230
231	! Adding eddy mixing to mimic 3D general circulation in 1D
232	! R. Wordsworth & F. Forget (2010)
233	if ((ngrid.eq.1)) then
234	! kmixmin minimum eddy mix coeff in 1D
235	! set up in inifis_mod.F90 - default value 1.0e-2
236	do ilev=1,nlay
237
238	! Here to code your specific eddy mix coeff in 1D
239	! Earth example that can be uncommented below
240	! -------------------------------------------------
241	! ====== Earth kzz from Zahnle et al. 2006 ======
242	! -------------------------------------------------
243	! if(pzlev(1,ilev).le.11.0e3) then
244	! kzz_eddy(ilev)=10.0
245	! pmin_kzz=pplev(1,ilev)exp((pzlev(1,ilev)-11.0e3)g/(r*zt(1,ilev)))
246	! else
247	! kzz_eddy(ilev)=0.1(pplev(1,ilev)/pmin_kzz)*(-0.5)
248	! kzz_eddy(ilev)=min(kzz_eddy(ilev),100.0)
249	! endif
250	! do ig=1,ngrid
251	! zkh(ig,ilev) = max(kzz_eddy(ilev),zkh(ig,ilev))
252	! zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev))
253	! end do
254
255	do ig=1,ngrid
256	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
257	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
258	end do
259	end do
260	end if
261
262	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
263	DO ig=1,ngrid
264	zkv(ig,1)=zcdv(ig)
265	ENDDO
266	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
267	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
268	zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
269
270	!-----------------------------------------------------------------------
271	! 4. Implicit inversion of u
272	! --------------------------
273
274	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
275	! avec
276	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
277	! et
278	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
279	! donc les entrees sont /zcdv/ pour la condition a la limite sol
280	! et /zkv/ = Ku
281
282	DO ig=1,ngrid
283	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
284	zcv(ig,nlay)=zmass(ig,nlay)zu(ig,nlay)z1(ig)
285	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
286	ENDDO
287
288	DO ilay=nlay-1,1,-1
289	DO ig=1,ngrid
290	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
291	zcv(ig,ilay)=(zmass(ig,ilay)zu(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
292	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
293	ENDDO
294	ENDDO
295
296	DO ig=1,ngrid
297	zu(ig,1)=zcv(ig,1)
298	ENDDO
299	DO ilay=2,nlay
300	DO ig=1,ngrid
301	zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1)
302	ENDDO
303	ENDDO
304
305	!-----------------------------------------------------------------------
306	! 5. Implicit inversion of v
307	! --------------------------
308
309	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
310	! avec
311	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
312	! et
313	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
314	! donc les entrees sont /zcdv/ pour la condition a la limite sol
315	! et /zkv/ = Kv
316
317	DO ig=1,ngrid
318	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
319	zcv(ig,nlay)=zmass(ig,nlay)zv(ig,nlay)z1(ig)
320	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
321	ENDDO
322
323	DO ilay=nlay-1,1,-1
324	DO ig=1,ngrid
325	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
326	zcv(ig,ilay)=(zmass(ig,ilay)zv(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
327	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
328	ENDDO
329	ENDDO
330
331	DO ig=1,ngrid
332	zv(ig,1)=zcv(ig,1)
333	ENDDO
334	DO ilay=2,nlay
335	DO ig=1,ngrid
336	zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1)
337	ENDDO
338	ENDDO
339
340	! Calcul of wind stress
341
342	DO ig=1,ngrid
343	flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
344	flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
345	ENDDO
346
347
348	!----------------------------------------------------------------------------
349	! 6. Implicit inversion of h, not forgetting the coupling with the ground
350
351	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
352	! avec
353	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
354	! et
355	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
356	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
357	! et /zkh/ = Kh
358
359	! Using the wind modified by friction for lifting and sublimation
360	! ---------------------------------------------------------------
361	DO ig=1,ngrid
362	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
363	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
364	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
365	zkh(ig,1)= zcdh(ig)
366	ustar(ig)=sqrt(zcdv_true(ig))*sqrt(zu2)
367	ENDDO
368
369
370	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
371	! ---------------------------------------------------------------
372	zfluxq(1:ngrid,1:nlay)=zkh(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay) !JL12 we save zfluxq which doesn't need the Exner factor
373	zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
374
375	DO ig=1,ngrid
376	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
377	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
378	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
379	ENDDO
380
381	DO ilay=nlay-1,2,-1
382	DO ig=1,ngrid
383	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
384	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
385	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
386	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
387	ENDDO
388	ENDDO
389
390	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
391	DO ig=1,ngrid
392	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
393	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
394	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
395	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
396	ENDDO
397
398
399	! Calculate (d Planck / dT) at the interface temperature
400	! ------------------------------------------------------
401
402	z4st=4.0sigmaptimestep
403	DO ig=1,ngrid
404	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
405	ENDDO
406
407	! Calculate temperature tendency at the interface (dry case)
408	! ----------------------------------------------------------
409	! Sum of fluxes at interface at time t + \delta t gives change in T:
410	! radiative fluxes
411	! turbulent convective (sensible) heat flux
412	! flux (if any) from subsurface
413
414	! if(.not.water) then
415
416	DO ig=1,ngrid
417	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
418	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
419	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
420	ztsrf(ig) = z1(ig) / z2(ig)
421	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
422	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
423	ENDDO
424	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
425
426
427	! Recalculate temperature to top of atmosphere, starting from ground
428	! ------------------------------------------------------------------
429
430	DO ilay=2,nlay
431	DO ig=1,ngrid
432	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
433	ENDDO
434	ENDDO
435
436	! endif ! not water
437
438	!-----------------------------------------------------------------------
439	! TRACERS (no vapour)
440	! -------
441
442	if(tracer) then
443
444	! Calculate vertical flux from the bottom to the first layer (dust)
445	! -----------------------------------------------------------------
446	do ig=1,ngrid
447	rho(ig) = zb0(ig,1) /ptimestep
448	end do
449
450	pdqsdif(:,:)=0.
451
452	! Implicit inversion of q
453	! -----------------------
454	do iq=1,nq
455
456	if (iq.ne.ivap) then
457
458	DO ig=1,ngrid
459	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
460	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
461	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
462	ENDDO
463
464	DO ilay=nlay-1,2,-1
465	DO ig=1,ngrid
466	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
467	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
468	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
469	ENDDO
470	ENDDO
471
472	do ig=1,ngrid
473	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
474	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
475	! tracer flux from surface
476	! currently pdqsdif always zero here,
477	! so last line is superfluous
478	enddo
479
480	! Starting upward calculations for simple tracer mixing (e.g., dust)
481	do ig=1,ngrid
482	zq(ig,1,iq)=zcq(ig,1)
483	end do
484
485	do ilay=2,nlay
486	do ig=1,ngrid
487	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
488	end do
489	end do
490
491	endif ! if (iq.ne.ivap)
492	end do ! of do iq=1,nq
493
494	! Revmoed (AF24): Calculate temperature tendency including latent heat term
495	! and assuming an infinite source of water on the ground
496	! ------------------------------------------------------------------
497
498	endif ! tracer
499
500
501	!-----------------------------------------------------------------------
502	! 8. Final calculation of the vertical diffusion tendencies
503	! -----------------------------------------------------------------
504
505	do ilev = 1, nlay
506	do ig=1,ngrid
507	pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
508	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
509	pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
510	enddo
511	enddo
512
513	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
514	sensibFlux(ig)=cppzfluxt(ig,1)/ptimestep(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
515	ENDDO
516
517	if (tracer) then
518	do iq = 1, nq
519	do ilev = 1, nlay
520	do ig=1,ngrid
521	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
522	enddo
523	enddo
524	enddo
525	endif
526	end subroutine turbdiff
527
528	end module turbdiff_mod

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

Download in other formats:

Original Format