Context Navigation

turbdiff.F90 @ 2176

Last change on this file since 2176 was 1993, checked in by jleconte, 6 years ago

29/08/2018 == JL

-watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...) which used Psat_water from watercommon.
This is now harmonized. ALl routines use Psat_water. Watersat.F has been removed, but the routine is now in watercommon for archival purpose. It is not used anymore.
-also changed the number of chars for tname in the dyn3D/infotrac.F90 to be able to run rcm1d.

File size: 27.8 KB

Line
1	subroutine turbdiff(ngrid,nlay,nq,rnat, &
2	ptimestep,pcapcal,lecrit, &
3	pplay,pplev,pzlay,pzlev,pz0, &
4	pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf, &
5	pdtfi,pdqfi,pfluxsrf, &
6	Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, &
7	pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
8
9	use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water, Lcpdqsat_water
10	use radcommon_h, only : sigma, glat
11	use surfdat_h, only: dryness
12	use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
13	use comcstfi_mod, only: rcp, g, r, cpp
14	use callkeys_mod, only: water,tracer,nosurf
15	use turb_mod, only : ustar
16	#ifdef MESOSCALE
17	use comm_wrf, only : comm_LATENT_HF
18	#endif
19	implicit none
20
21	!==================================================================
22	!
23	! Purpose
24	! -------
25	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
26	!
27	! Implicit scheme
28	! We start by adding to variables x the physical tendencies
29	! already computed. We resolve the equation:
30	!
31	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
32	!
33	! Authors
34	! -------
35	! F. Hourdin, F. Forget, R. Fournier (199X)
36	! R. Wordsworth, B. Charnay (2010)
37	! J. Leconte (2012): To f90
38	! - Rewritten the diffusion scheme to conserve total enthalpy
39	! by accounting for dissipation of turbulent kinetic energy.
40	! - Accounting for lost mean flow kinetic energy should come soon.
41	!
42	!==================================================================
43
44	!-----------------------------------------------------------------------
45	! declarations
46	! ------------
47
48	! arguments
49	! ---------
50	INTEGER,INTENT(IN) :: ngrid
51	INTEGER,INTENT(IN) :: nlay
52	REAL,INTENT(IN) :: ptimestep
53	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
54	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
55	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
56	REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
57	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
58	REAL,INTENT(IN) :: pemis(ngrid)
59	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
60	REAL,INTENT(IN) :: pfluxsrf(ngrid)
61	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
62	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
63	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
64	REAL,INTENT(OUT) :: sensibFlux(ngrid)
65	REAL,INTENT(IN) :: pcapcal(ngrid)
66	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
67	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
68	REAL,INTENT(IN) :: rnat(ngrid)
69	LOGICAL,INTENT(IN) :: lastcall ! not used
70
71	! Arguments added for condensation
72	logical,intent(in) :: lecrit ! not used.
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
107	LOGICAL,SAVE :: firstcall=.true.
108	!$OMP THREADPRIVATE(firstcall)
109
110	! Tracers
111	! -------
112	INTEGER iq
113	REAL zq(ngrid,nlay,nq)
114	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
115	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
116	REAL zdmassevap(ngrid)
117	REAL rho(ngrid) ! near-surface air density
118	REAL kmixmin
119
120	! Variables added for implicit latent heat inclusion
121	! --------------------------------------------------
122	real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
123
124	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
125	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
126
127	real, parameter :: karman=0.4
128	real cd0, roughratio
129
130	real dqsdif_total(ngrid)
131	real zq0(ngrid)
132
133
134	! Coherence test
135	! --------------
136
137	IF (firstcall) THEN
138
139	if(water)then
140	ivap=igcm_h2o_vap
141	iliq=igcm_h2o_ice
142	iliq_surf=igcm_h2o_vap
143	iice_surf=igcm_h2o_ice ! simply to make the code legible
144	! to be generalised
145	else
146	ivap=0
147	iliq=0
148	iliq_surf=0
149	iice_surf=0 ! simply to make the code legible
150	endif
151	sensibFlux(:)=0.
152
153	firstcall=.false.
154	ENDIF
155
156	!-----------------------------------------------------------------------
157	! 1. Initialisation
158	! -----------------
159
160	nlev=nlay+1
161
162	! Calculate rhodz, (P/Ps)(R/cp) 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	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
169	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
170	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
171	ENDDO
172	ENDDO
173
174	zcst1=4.gptimestep/(R*R)
175	DO ilev=2,nlev-1
176	DO ig=1,ngrid
177	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
178	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
179	ENDDO
180	ENDDO
181	DO ig=1,ngrid
182	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
183	ENDDO
184	dqsdif_total(:)=0.0
185
186	!-----------------------------------------------------------------------
187	! 2. Add the physical tendencies computed so far
188	! ----------------------------------------------
189
190	DO ilev=1,nlay
191	DO ig=1,ngrid
192	zu(ig,ilev)=pu(ig,ilev)
193	zv(ig,ilev)=pv(ig,ilev)
194	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
195	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
196	ENDDO
197	ENDDO
198	if(tracer) then
199	DO iq =1, nq
200	DO ilev=1,nlay
201	DO ig=1,ngrid
202	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
203	ENDDO
204	ENDDO
205	ENDDO
206	if (water) then
207	DO ilev=1,nlay
208	DO ig=1,ngrid
209	zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
210	ENDDO
211	ENDDO
212	Endif
213	end if
214
215	!-----------------------------------------------------------------------
216	! 3. Turbulence scheme
217	! --------------------
218	!
219	! Source of turbulent kinetic energy at the surface
220	! -------------------------------------------------
221	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
222
223	DO ig=1,ngrid
224	roughratio = 1. + pzlay(ig,1)/pz0
225	cd0 = karman/log(roughratio)
226	cd0 = cd0*cd0
227	zcdv_true(ig) = cd0
228	zcdh_true(ig) = cd0
229	if(nosurf)then
230	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
231	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
232	endif
233	ENDDO
234
235	DO ig=1,ngrid
236	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
237	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
238	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
239	ENDDO
240
241	! Turbulent diffusion coefficients in the boundary layer
242	! ------------------------------------------------------
243
244	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
245
246	! Adding eddy mixing to mimic 3D general circulation in 1D
247	! R. Wordsworth & F. Forget (2010)
248	if ((ngrid.eq.1)) then
249	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
250	do ilev=1,nlay
251	do ig=1,ngrid
252	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
253	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
254	end do
255	end do
256	end if
257
258	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
259	DO ig=1,ngrid
260	zkv(ig,1)=zcdv(ig)
261	ENDDO
262	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
263	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
264	zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
265
266	!-----------------------------------------------------------------------
267	! 4. Implicit inversion of u
268	! --------------------------
269
270	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
271	! avec
272	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
273	! et
274	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
275	! donc les entrees sont /zcdv/ pour la condition a la limite sol
276	! et /zkv/ = Ku
277
278	DO ig=1,ngrid
279	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
280	zcv(ig,nlay)=zmass(ig,nlay)zu(ig,nlay)z1(ig)
281	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
282	ENDDO
283
284	DO ilay=nlay-1,1,-1
285	DO ig=1,ngrid
286	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
287	zcv(ig,ilay)=(zmass(ig,ilay)zu(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
288	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
289	ENDDO
290	ENDDO
291
292	DO ig=1,ngrid
293	zu(ig,1)=zcv(ig,1)
294	ENDDO
295	DO ilay=2,nlay
296	DO ig=1,ngrid
297	zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1)
298	ENDDO
299	ENDDO
300
301	!-----------------------------------------------------------------------
302	! 5. Implicit inversion of v
303	! --------------------------
304
305	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
306	! avec
307	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
308	! et
309	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
310	! donc les entrees sont /zcdv/ pour la condition a la limite sol
311	! et /zkv/ = Kv
312
313	DO ig=1,ngrid
314	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
315	zcv(ig,nlay)=zmass(ig,nlay)zv(ig,nlay)z1(ig)
316	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
317	ENDDO
318
319	DO ilay=nlay-1,1,-1
320	DO ig=1,ngrid
321	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
322	zcv(ig,ilay)=(zmass(ig,ilay)zv(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
323	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
324	ENDDO
325	ENDDO
326
327	DO ig=1,ngrid
328	zv(ig,1)=zcv(ig,1)
329	ENDDO
330	DO ilay=2,nlay
331	DO ig=1,ngrid
332	zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1)
333	ENDDO
334	ENDDO
335
336	! Calcul of wind stress
337
338	DO ig=1,ngrid
339	flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
340	flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
341	ENDDO
342
343
344	!----------------------------------------------------------------------------
345	! 6. Implicit inversion of h, not forgetting the coupling with the ground
346
347	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
348	! avec
349	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
350	! et
351	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
352	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
353	! et /zkh/ = Kh
354
355	! Using the wind modified by friction for lifting and sublimation
356	! ---------------------------------------------------------------
357	DO ig=1,ngrid
358	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
359	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
360	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
361	zkh(ig,1)= zcdh(ig)
362	ustar(ig)=sqrt(zcdv_true(ig))*sqrt(zu2)
363	ENDDO
364
365
366	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
367	! ---------------------------------------------------------------
368	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
369	zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
370
371	DO ig=1,ngrid
372	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
373	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
374	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
375	ENDDO
376
377	DO ilay=nlay-1,2,-1
378	DO ig=1,ngrid
379	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
380	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
381	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
382	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
383	ENDDO
384	ENDDO
385
386	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
387	DO ig=1,ngrid
388	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
389	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
390	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
391	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
392	ENDDO
393
394
395	! Calculate (d Planck / dT) at the interface temperature
396	! ------------------------------------------------------
397
398	z4st=4.0sigmaptimestep
399	DO ig=1,ngrid
400	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
401	ENDDO
402
403	! Calculate temperature tendency at the interface (dry case)
404	! ----------------------------------------------------------
405	! Sum of fluxes at interface at time t + \delta t gives change in T:
406	! radiative fluxes
407	! turbulent convective (sensible) heat flux
408	! flux (if any) from subsurface
409
410	if(.not.water) then
411
412	DO ig=1,ngrid
413	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
414	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
415	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
416	ztsrf(ig) = z1(ig) / z2(ig)
417	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
418	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
419	ENDDO
420	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
421
422
423	! Recalculate temperature to top of atmosphere, starting from ground
424	! ------------------------------------------------------------------
425
426	DO ilay=2,nlay
427	DO ig=1,ngrid
428	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
429	ENDDO
430	ENDDO
431
432	endif ! not water
433
434	!-----------------------------------------------------------------------
435	! TRACERS (no vapour)
436	! -------
437
438	if(tracer) then
439
440	! Calculate vertical flux from the bottom to the first layer (dust)
441	! -----------------------------------------------------------------
442	do ig=1,ngrid
443	rho(ig) = zb0(ig,1) /ptimestep
444	end do
445
446	pdqsdif(:,:)=0.
447
448	! Implicit inversion of q
449	! -----------------------
450	do iq=1,nq
451
452	if (iq.ne.ivap) then
453
454	DO ig=1,ngrid
455	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
456	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
457	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
458	ENDDO
459
460	DO ilay=nlay-1,2,-1
461	DO ig=1,ngrid
462	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
463	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
464	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
465	ENDDO
466	ENDDO
467
468	if ((water).and.(iq.eq.iliq)) then
469	! special case for condensed water tracer: do not include
470	! h2o ice tracer from surface (which is set when handling
471	! h2o vapour case (see further down).
472	! zb(ig,1)=0 if iq ne ivap
473	DO ig=1,ngrid
474	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
475	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
476	ENDDO
477	else ! general case
478	do ig=1,ngrid
479	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
480	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
481	! tracer flux from surface
482	! currently pdqsdif always zero here,
483	! so last line is superfluous
484	enddo
485	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
486
487
488	! Starting upward calculations for simple tracer mixing (e.g., dust)
489	do ig=1,ngrid
490	zq(ig,1,iq)=zcq(ig,1)
491	end do
492
493	do ilay=2,nlay
494	do ig=1,ngrid
495	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
496	end do
497	end do
498
499	endif ! if (iq.ne.ivap)
500
501	! Calculate temperature tendency including latent heat term
502	! and assuming an infinite source of water on the ground
503	! ------------------------------------------------------------------
504
505	if (water.and.(iq.eq.ivap)) then
506
507	! compute evaporation efficiency
508	do ig=1,ngrid
509	if(nint(rnat(ig)).eq.1)then
510	dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf)
511	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
512	dryness(ig)=MAX(0.,dryness(ig))
513	endif
514	enddo
515
516	do ig=1,ngrid
517	! Calculate the value of qsat at the surface (water)
518	call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
519	call Lcpdqsat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig),dqsat(ig),psat_temp)
520	dqsat(ig)=dqsat(ig)*RCPD/RLVTT
521	enddo
522
523	! coefficients for q
524
525	do ig=1,ngrid
526	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
527	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
528	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
529	enddo
530
531	do ilay=nlay-1,2,-1
532	do ig=1,ngrid
533	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
534	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
535	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
536	enddo
537	enddo
538
539	do ig=1,ngrid
540	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,1)dryness(ig)+zfluxq(ig,2)(1.-zdq(ig,2)))
541	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
542	zdq(ig,1)=dryness(ig)zfluxq(ig,1)z1(ig)
543	enddo
544
545	do ig=1,ngrid
546	!calculation of surface temperature
547	zdq0(ig) = dqsat(ig)
548	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
549
550	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
551	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
552	+ zfluxq(ig,1)dryness(ig)RLVTT((zdq(ig,1)-1.0)zcq0(ig)+zcq(ig,1))
553
554	z2(ig) = pcapcal(ig) + cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
555	+ zdplanck(ig)+zfluxq(ig,1)dryness(ig)RLVTTzdq0(ig)(1.0-zdq(ig,1))
556
557	ztsrf(ig) = z1(ig) / z2(ig)
558
559	! calculation of qs and q1
560	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf(ig)
561	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
562
563	! calculation of evaporation
564	dqsdif_total(ig)=zfluxq(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
565
566	! --------------------------------------------------------
567	! Now check if we've taken too much water from the surface
568	! This can only occur on the continent
569	! If we do, we recompute Tsurf, T1 and q1 accordingly
570	if((-dqsdif_total(ig).gt.(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))).and.rnat(ig).eq.1)then
571	!water flux * ptimestep
572	dqsdif_total(ig)=-(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))
573
574	!recompute surface temperature
575	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxq(ig,1)zct(ig,1)zovExner(ig,1) &
576	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
577	+ RLVTT*dqsdif_total(ig)
578	z2(ig) = pcapcal(ig) + cppzfluxq(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
579	+ zdplanck(ig)
580	ztsrf(ig) = z1(ig) / z2(ig)
581
582	!recompute q1 with new water flux from surface
583	zq(ig,1,iq) = (zmass(ig,1)(pq(ig,1,iq)+ptimesteppdqfi(ig,1,iq)) &
584	+zfluxq(ig,2)*zcq(ig,2)-dqsdif_total(ig)) &
585	/ (zmass(ig,1)+(1.-zdq(ig,2))*zfluxq(ig,2))
586	end if
587
588	! calculation surface T tendency and T(1)
589	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
590	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
591	enddo
592
593
594	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
595	do ilay=2,nlay
596	do ig=1,ngrid
597	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
598	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
599	end do
600	end do
601
602
603	do ig=1,ngrid
604	! --------------------------------------------------------------------------
605	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
606	! The surface vapour tracer is actually liquid. To make things difficult.
607
608	if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
609
610	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
611	pdqsdif(ig,iice_surf)=0.0
612
613	elseif (nint(rnat(ig)).eq.1) then ! (continent)
614	! If water is evaporating / subliming, we take it from ice before liquid
615	! -- is this valid??
616	if(dqsdif_total(ig).lt.0)then
617	if (-dqsdif_total(ig).gt.pqsurf(ig,iice_surf))then
618	pdqsdif(ig,iice_surf) = -pqsurf(ig,iice_surf)/ptimestep ! removes all the ice!
619	pdqsdif(ig,iliq_surf) = dqsdif_total(ig)/ptimestep- pdqsdif(ig,iice_surf) ! take the remainder from the liquid instead
620	else
621	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
622	pdqsdif(ig,iliq_surf)=0.
623	end if
624	else !dqsdif_total(ig).ge.0
625	!If water vapour is condensing, we must decide whether it forms ice or liquid.
626	if(ztsrf(ig).gt.T_h2O_ice_liq)then
627	pdqsdif(ig,iice_surf)=0.0
628	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
629	else
630	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
631	pdqsdif(ig,iliq_surf)=0.0
632	endif
633	endif
634
635	elseif (nint(rnat(ig)).eq.2) then ! (continental glaciers)
636	pdqsdif(ig,iliq_surf)=0.0
637	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
638
639	endif !rnat
640	end do ! of DO ig=1,ngrid
641
642	endif ! if (water et iq=ivap)
643	end do ! of do iq=1,nq
644
645	if (water) then ! special case where we recompute water mixing without any evaporation.
646	! The difference with the first calculation then tells us where evaporated water has gone
647
648	DO ig=1,ngrid
649	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
650	zcq(ig,nlay)=zmass(ig,nlay)zqnoevap(ig,nlay)z1(ig)
651	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
652	ENDDO
653
654	DO ilay=nlay-1,2,-1
655	DO ig=1,ngrid
656	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
657	zcq(ig,ilay)=(zmass(ig,ilay)zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
658	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
659	ENDDO
660	ENDDO
661
662	do ig=1,ngrid
663	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
664	zcq(ig,1)=(zmass(ig,1)zqnoevap(ig,1)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
665	enddo
666
667	! Starting upward calculations for simple tracer mixing (e.g., dust)
668	do ig=1,ngrid
669	zqnoevap(ig,1)=zcq(ig,1)
670	end do
671
672	do ilay=2,nlay
673	do ig=1,ngrid
674	zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
675	end do
676	end do
677
678	endif ! if water
679
680
681	endif ! tracer
682
683
684	!-----------------------------------------------------------------------
685	! 8. Final calculation of the vertical diffusion tendencies
686	! -----------------------------------------------------------------
687
688	do ilev = 1, nlay
689	do ig=1,ngrid
690	pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
691	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
692	pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
693	enddo
694	enddo
695
696	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
697	sensibFlux(ig)=cppzfluxt(ig,1)/ptimestep(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
698	ENDDO
699
700	if (tracer) then
701	do iq = 1, nq
702	do ilev = 1, nlay
703	do ig=1,ngrid
704	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
705	enddo
706	enddo
707	enddo
708	if (water) then
709	do ilev = 1, nlay
710	do ig=1,ngrid
711	pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
712	enddo
713	enddo
714	do ig=1,ngrid
715	zdmassevap(ig)=SUM(pdqevap(ig,:)zmass(ig,:))ptimestep
716	end do
717	endif
718	endif
719
720	if(water)then
721	#ifndef MESOSCALE
722	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
723	#endif
724	if (tracer) then
725	#ifndef MESOSCALE
726	call writediagfi(ngrid,'evap_surf_flux','surface latent heat flux','W.m-2',2,RLVTT*dqsdif_total/ptimestep)
727	call writediagfi(ngrid,'fluxsurf_rad','total IR and VIS surface flux','W.m-2',2,pfluxsrf)
728	call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
729	#else
730	comm_LATENT_HF(:)=0.0
731	comm_LATENT_HF(1:ngrid)=RLVTT*dqsdif_total(1:ngrid)/ptimestep
732	#endif
733	endif
734	endif
735
736	! if(lastcall)then
737	! if(ngrid.eq.1)then
738	! print*,'Saving k.out...'
739	! OPEN(12,file='k.out',form='formatted')
740	! DO ilay=1,nlay
741	! write(12,*) zkh(1,ilay), pplay(1,ilay)
742	! ENDDO
743	! CLOSE(12)
744	! endif
745	! endif
746
747	end

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

Download in other formats:

Original Format