Context Navigation

turbdiff.F90 @ 1243

Last change on this file since 1243 was 1194, checked in by sglmd, 11 years ago
Latitude-dependent gravity field added. Option oblate = .true. in callphys.def, and three additional variables needed: J2, Rmean and MassPlanet?.
File size: 27.0 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	pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf, &
6	Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, &
7	pdqdif,pdqevap,pdqsdif,lastcall)
8
9	use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
10	use radcommon_h, only : sigma, glat
11	USE surfdat_h
12	USE comgeomfi_h
13	USE tracer_h
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	! J. Leconte (2012): To f90
34	! - Rewritten the diffusion scheme to conserve total enthalpy
35	! by accounting for dissipation of turbulent kinetic energy.
36	! - Accounting for lost mean flow kinetic energy should come soon.
37	!
38	!==================================================================
39
40	!-----------------------------------------------------------------------
41	! declarations
42	! ------------
43
44	#include "dimensions.h"
45	#include "dimphys.h"
46	#include "comcstfi.h"
47	#include "callkeys.h"
48
49	! arguments
50	! ---------
51	INTEGER ngrid,nlay
52	REAL ptimestep
53	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
54	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
55	REAL pu(ngrid,nlay),pv(ngrid,nlay)
56	REAL pt(ngrid,nlay),ppopsk(ngrid,nlay)
57	REAL ptsrf(ngrid),pemis(ngrid)
58	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
59	REAL pdtfi(ngrid,nlay)
60	REAL pfluxsrf(ngrid)
61	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
62	REAL pdtdif(ngrid,nlay)
63	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
64	REAL pq2(ngrid,nlay+1)
65
66	integer rnat(ngrid)
67
68	! Arguments added for condensation
69	logical lecrit
70	REAL pz0
71
72	! Tracers
73	! --------
74	integer nq
75	real pqsurf(ngrid,nq)
76	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
77	real pdqdif(ngrid,nlay,nq)
78	real pdqsdif(ngrid,nq)
79
80	! local
81	! -----
82	integer ilev,ig,ilay,nlev
83
84	REAL z4st,zdplanck(ngrid)
85	REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
86	REAL zcdv(ngrid),zcdh(ngrid)
87	REAL zcdv_true(ngrid),zcdh_true(ngrid)
88	REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
89	REAL zh(ngrid,nlayermx),zt(ngrid,nlayermx)
90	REAL ztsrf(ngrid)
91	REAL z1(ngrid),z2(ngrid)
92	REAL zmass(ngrid,nlayermx)
93	REAL zfluxv(ngrid,nlayermx),zfluxt(ngrid,nlayermx),zfluxq(ngrid,nlayermx)
94	REAL zb0(ngrid,nlayermx)
95	REAL zExner(ngrid,nlayermx),zovExner(ngrid,nlayermx)
96	REAL zcv(ngrid,nlayermx),zdv(ngrid,nlayermx) !inversion coefficient for winds
97	REAL zct(ngrid,nlayermx),zdt(ngrid,nlayermx) !inversion coefficient for temperature
98	REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx) !inversion coefficient for tracers
99	REAL zcst1
100	REAL zu2!, a
101	REAL zcq0(ngrid),zdq0(ngrid)
102	REAL zx_alf1(ngrid),zx_alf2(ngrid)
103
104	LOGICAL firstcall
105	SAVE firstcall
106
107	LOGICAL lastcall
108	! Tracers
109	! -------
110	INTEGER iq
111	REAL zq(ngrid,nlayermx,nq)
112	REAL zqnoevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
113	REAL pdqevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
114	REAL zdmassevap(ngrid)
115	REAL rho(ngrid) ! near-surface air density
116	REAL qsat(ngrid),psat(ngrid)
117	DATA firstcall/.true./
118	REAL kmixmin
119
120	! Variables added for implicit latent heat inclusion
121	! --------------------------------------------------
122	real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
123
124	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
125
126	real, parameter :: karman=0.4
127	real cd0, roughratio
128
129	real dqsdif_total(ngrid)
130	real zq0(ngrid)
131
132
133	! Coherence test
134	! --------------
135
136	IF (firstcall) THEN
137
138	if(water)then
139	ivap=igcm_h2o_vap
140	iliq=igcm_h2o_ice
141	iliq_surf=igcm_h2o_vap
142	iice_surf=igcm_h2o_ice ! simply to make the code legible
143	! to be generalised
144	else
145	ivap=0
146	iliq=0
147	iliq_surf=0
148	iice_surf=0 ! simply to make the code legible
149	endif
150	sensibFlux(:)=0.
151
152	firstcall=.false.
153	ENDIF
154
155	!-----------------------------------------------------------------------
156	! 1. Initialisation
157	! -----------------
158
159	nlev=nlay+1
160
161	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
162	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
163	! ---------------------------------------------
164
165	DO ilay=1,nlay
166	DO ig=1,ngrid
167	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
168	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
169	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
170	ENDDO
171	ENDDO
172
173	zcst1=4.gptimestep/(R*R)
174	DO ilev=2,nlev-1
175	DO ig=1,ngrid
176	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
177	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
178	ENDDO
179	ENDDO
180	DO ig=1,ngrid
181	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
182	ENDDO
183	dqsdif_total(:)=0.0
184
185	!-----------------------------------------------------------------------
186	! 2. Add the physical tendencies computed so far
187	! ----------------------------------------------
188
189	DO ilev=1,nlay
190	DO ig=1,ngrid
191	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
192	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
193	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
194	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
195	ENDDO
196	ENDDO
197	if(tracer) then
198	DO iq =1, nq
199	DO ilev=1,nlay
200	DO ig=1,ngrid
201	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
202	ENDDO
203	ENDDO
204	ENDDO
205	if (water) then
206	DO ilev=1,nlay
207	DO ig=1,ngrid
208	zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
209	ENDDO
210	ENDDO
211	Endif
212	end if
213
214	!-----------------------------------------------------------------------
215	! 3. Turbulence scheme
216	! --------------------
217	!
218	! Source of turbulent kinetic energy at the surface
219	! -------------------------------------------------
220	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
221
222	DO ig=1,ngrid
223	roughratio = 1. + pzlay(ig,1)/pz0
224	cd0 = karman/log(roughratio)
225	cd0 = cd0*cd0
226	zcdv_true(ig) = cd0
227	zcdh_true(ig) = cd0
228	if(nosurf)then
229	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
230	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
231	endif
232	ENDDO
233
234	DO ig=1,ngrid
235	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
236	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
237	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
238	ENDDO
239
240	! Turbulent diffusion coefficients in the boundary layer
241	! ------------------------------------------------------
242
243	call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
244
245	! Adding eddy mixing to mimic 3D general circulation in 1D
246	! R. Wordsworth & F. Forget (2010)
247	if ((ngrid.eq.1)) then
248	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
249	do ilev=1,nlay
250	do ig=1,ngrid
251	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
252	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
253	end do
254	end do
255	end if
256
257	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
258	DO ig=1,ngrid
259	zkv(ig,1)=zcdv(ig)
260	ENDDO
261	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
262
263	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
264	zfluxv(1:ngrid,1:nlayermx)=zkv(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx)
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
337
338	!----------------------------------------------------------------------------
339	! 6. Implicit inversion of h, not forgetting the coupling with the ground
340
341	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
342	! avec
343	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
344	! et
345	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
346	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
347	! et /zkh/ = Kh
348
349	! Using the wind modified by friction for lifting and sublimation
350	! ---------------------------------------------------------------
351	DO ig=1,ngrid
352	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
353	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
354	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
355	zkh(ig,1)= zcdh(ig)
356	ENDDO
357
358	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
359	! ---------------------------------------------------------------
360	zfluxq(1:ngrid,1:nlayermx)=zkh(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
361	zfluxt(1:ngrid,1:nlayermx)=zfluxq(1:ngrid,1:nlayermx)*zExner(1:ngrid,1:nlayermx)
362
363	DO ig=1,ngrid
364	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
365	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
366	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
367	ENDDO
368
369	DO ilay=nlay-1,2,-1
370	DO ig=1,ngrid
371	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
372	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
373	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
374	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
375	ENDDO
376	ENDDO
377
378	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
379	DO ig=1,ngrid
380	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
381	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
382	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
383	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
384	ENDDO
385
386
387	! Calculate (d Planck / dT) at the interface temperature
388	! ------------------------------------------------------
389
390	z4st=4.0sigmaptimestep
391	DO ig=1,ngrid
392	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
393	ENDDO
394
395	! Calculate temperature tendency at the interface (dry case)
396	! ----------------------------------------------------------
397	! Sum of fluxes at interface at time t + \delta t gives change in T:
398	! radiative fluxes
399	! turbulent convective (sensible) heat flux
400	! flux (if any) from subsurface
401
402	if(.not.water) then
403
404	DO ig=1,ngrid
405	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
406	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
407	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
408	ztsrf(ig) = z1(ig) / z2(ig)
409	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
410	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
411	ENDDO
412	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
413
414
415	! Recalculate temperature to top of atmosphere, starting from ground
416	! ------------------------------------------------------------------
417
418	DO ilay=2,nlay
419	DO ig=1,ngrid
420	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
421	ENDDO
422	ENDDO
423
424	endif ! not water
425
426	!-----------------------------------------------------------------------
427	! TRACERS (no vapour)
428	! -------
429
430	if(tracer) then
431
432	! Calculate vertical flux from the bottom to the first layer (dust)
433	! -----------------------------------------------------------------
434	do ig=1,ngrid
435	rho(ig) = zb0(ig,1) /ptimestep
436	end do
437
438	pdqsdif(:,:)=0.
439
440	! Implicit inversion of q
441	! -----------------------
442	do iq=1,nq
443
444	if (iq.ne.ivap) then
445
446	DO ig=1,ngrid
447	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
448	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
449	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
450	ENDDO
451
452	DO ilay=nlay-1,2,-1
453	DO ig=1,ngrid
454	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
455	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
456	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
457	ENDDO
458	ENDDO
459
460	if ((water).and.(iq.eq.iliq)) then
461	! special case for condensed water tracer: do not include
462	! h2o ice tracer from surface (which is set when handling
463	! h2o vapour case (see further down).
464	! zb(ig,1)=0 if iq ne ivap
465	DO ig=1,ngrid
466	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
467	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
468	ENDDO
469	else ! general case
470	do ig=1,ngrid
471	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
472	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
473	! tracer flux from surface
474	! currently pdqsdif always zero here,
475	! so last line is superfluous
476	enddo
477	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
478
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
493	! Calculate temperature tendency including latent heat term
494	! and assuming an infinite source of water on the ground
495	! ------------------------------------------------------------------
496
497	if (water.and.(iq.eq.ivap)) then
498
499	! compute evaporation efficiency
500	do ig=1,ngrid
501	if(rnat(ig).eq.1)then
502	dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf)
503	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
504	dryness(ig)=MAX(0.,dryness(ig))
505	endif
506	enddo
507
508	do ig=1,ngrid
509	! Calculate the value of qsat at the surface (water)
510	call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
511	call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
512	call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
513	dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
514	! calculate dQsat / dT by finite differences
515	! we cannot use the updated temperature value yet...
516	enddo
517
518	! coefficients for q
519
520	do ig=1,ngrid
521	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
522	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
523	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
524	enddo
525
526	do ilay=nlay-1,2,-1
527	do ig=1,ngrid
528	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
529	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
530	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
531	enddo
532	enddo
533
534	do ig=1,ngrid
535	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,1)dryness(ig)+zfluxq(ig,2)(1.-zdq(ig,2)))
536	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
537	zdq(ig,1)=dryness(ig)zfluxq(ig,1)z1(ig)
538	enddo
539
540	do ig=1,ngrid
541	!calculation of surface temperature
542	zdq0(ig) = dqsat(ig)
543	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
544
545	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
546	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
547	+ zfluxq(ig,1)dryness(ig)RLVTT((zdq(ig,1)-1.0)zcq0(ig)+zcq(ig,1))
548
549	z2(ig) = pcapcal(ig) + cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
550	+ zdplanck(ig)+zfluxq(ig,1)dryness(ig)RLVTTzdq0(ig)(1.0-zdq(ig,1))
551
552	ztsrf(ig) = z1(ig) / z2(ig)
553
554	! calculation of qs and q1
555	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf(ig)
556	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
557
558	! calculation of evaporation
559	dqsdif_total(ig)=zfluxq(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
560
561	! --------------------------------------------------------
562	! Now check if we've taken too much water from the surface
563	! This can only occur on the continent
564	! If we do, we recompute Tsurf, T1 and q1 accordingly
565	if((-dqsdif_total(ig).gt.(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))).and.rnat(ig).eq.1)then
566	!water flux * ptimestep
567	dqsdif_total(ig)=-(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))
568
569	!recompute surface temperature
570	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzfluxq(ig,1)zct(ig,1)zovExner(ig,1) &
571	+ zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep &
572	+ RLVTT*dqsdif_total(ig)
573	z2(ig) = pcapcal(ig) + cppzfluxq(ig,1)(1.-zovExner(ig,1)*zdt(ig,1)) &
574	+ zdplanck(ig)
575	ztsrf(ig) = z1(ig) / z2(ig)
576
577	!recompute q1 with new water flux from surface
578	zq(ig,1,iq) = (zmass(ig,1)(pq(ig,1,iq)+ptimesteppdqfi(ig,1,iq)) &
579	+zfluxq(ig,2)*zcq(ig,2)-dqsdif_total(ig)) &
580	/ (zmass(ig,1)+(1.-zdq(ig,2))*zfluxq(ig,2))
581	end if
582
583	! calculation surface T tendency and T(1)
584	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
585	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
586	enddo
587
588
589	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
590	do ilay=2,nlay
591	do ig=1,ngrid
592	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
593	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
594	end do
595	end do
596
597
598	do ig=1,ngrid
599	! --------------------------------------------------------------------------
600	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
601	! The surface vapour tracer is actually liquid. To make things difficult.
602
603	if (rnat(ig).eq.0) then ! unfrozen ocean
604
605	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
606	pdqsdif(ig,iice_surf)=0.0
607
608	elseif (rnat(ig).eq.1) then ! (continent)
609	! If water is evaporating / subliming, we take it from ice before liquid
610	! -- is this valid??
611	if(dqsdif_total(ig).lt.0)then
612	if (-dqsdif_total(ig).gt.pqsurf(ig,iice_surf))then
613	pdqsdif(ig,iice_surf) = -pqsurf(ig,iice_surf)/ptimestep ! removes all the ice!
614	pdqsdif(ig,iliq_surf) = dqsdif_total(ig)/ptimestep- pdqsdif(ig,iice_surf) ! take the remainder from the liquid instead
615	else
616	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
617	pdqsdif(ig,iliq_surf)=0.
618	end if
619	else !dqsdif_total(ig).ge.0
620	!If water vapour is condensing, we must decide whether it forms ice or liquid.
621	if(ztsrf(ig).gt.T_h2O_ice_liq)then
622	pdqsdif(ig,iice_surf)=0.0
623	pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
624	else
625	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
626	pdqsdif(ig,iliq_surf)=0.0
627	endif
628	endif
629
630	elseif (rnat(ig).eq.2) then ! (continental glaciers)
631	pdqsdif(ig,iliq_surf)=0.0
632	pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
633
634	endif !rnat
635	end do ! of DO ig=1,ngrid
636
637	endif ! if (water et iq=ivap)
638	end do ! of do iq=1,nq
639
640	if (water) then ! special case where we recompute water mixing without any evaporation.
641	! The difference with the first calculation then tells us where evaporated water has gone
642
643	DO ig=1,ngrid
644	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
645	zcq(ig,nlay)=zmass(ig,nlay)zqnoevap(ig,nlay)z1(ig)
646	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
647	ENDDO
648
649	DO ilay=nlay-1,2,-1
650	DO ig=1,ngrid
651	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
652	zcq(ig,ilay)=(zmass(ig,ilay)zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
653	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
654	ENDDO
655	ENDDO
656
657	do ig=1,ngrid
658	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
659	zcq(ig,1)=(zmass(ig,1)zqnoevap(ig,1)+zfluxq(ig,2)zcq(ig,2))*z1(ig)
660	enddo
661
662	! Starting upward calculations for simple tracer mixing (e.g., dust)
663	do ig=1,ngrid
664	zqnoevap(ig,1)=zcq(ig,1)
665	end do
666
667	do ilay=2,nlay
668	do ig=1,ngrid
669	zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
670	end do
671	end do
672
673	endif ! if water
674
675
676	endif ! tracer
677
678
679	!-----------------------------------------------------------------------
680	! 8. Final calculation of the vertical diffusion tendencies
681	! -----------------------------------------------------------------
682
683	do ilev = 1, nlay
684	do ig=1,ngrid
685	pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
686	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
687	pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
688	enddo
689	enddo
690
691	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
692	sensibFlux(ig)=cppzfluxt(ig,1)/ptimestep(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
693	ENDDO
694
695	if (tracer) then
696	do iq = 1, nq
697	do ilev = 1, nlay
698	do ig=1,ngrid
699	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
700	enddo
701	enddo
702	enddo
703	if (water) then
704	do ilev = 1, nlay
705	do ig=1,ngrid
706	pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
707	enddo
708	enddo
709	do ig=1,ngrid
710	zdmassevap(ig)=SUM(pdqevap(ig,:)zmass(ig,:))ptimestep
711	end do
712	endif
713	endif
714
715	if(water)then
716	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
717	if (tracer) then
718	call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
719	endif
720	endif
721
722	! if(lastcall)then
723	! if(ngrid.eq.1)then
724	! print*,'Saving k.out...'
725	! OPEN(12,file='k.out',form='formatted')
726	! DO ilay=1,nlay
727	! write(12,*) zkh(1,ilay), pplay(1,ilay)
728	! ENDDO
729	! CLOSE(12)
730	! endif
731	! endif
732
733
734	return
735	end

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

Download in other formats:

Original Format