Context Navigation

turbdiff.F90 @ 1626

Last change on this file since 1626 was 1477, checked in by mturbet, 9 years ago
Cleaning of physiq.F90. Condense_cloud becomes Condense_co2
File size: 27.5 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
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
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	! J. Leconte (2012): To f90
35	! - Rewritten the diffusion scheme to conserve total enthalpy
36	! by accounting for dissipation of turbulent kinetic energy.
37	! - Accounting for lost mean flow kinetic energy should come soon.
38	!
39	!==================================================================
40
41	!-----------------------------------------------------------------------
42	! declarations
43	! ------------
44
45	! arguments
46	! ---------
47	INTEGER,INTENT(IN) :: ngrid
48	INTEGER,INTENT(IN) :: nlay
49	REAL,INTENT(IN) :: ptimestep
50	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
51	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
52	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
53	REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
54	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
55	REAL,INTENT(IN) :: pemis(ngrid)
56	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
57	REAL,INTENT(IN) :: pfluxsrf(ngrid)
58	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
59	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
60	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
61	REAL,INTENT(OUT) :: sensibFlux(ngrid)
62	REAL,INTENT(IN) :: pcapcal(ngrid)
63	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
64	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
65	REAL,INTENT(IN) :: rnat(ngrid)
66	LOGICAL,INTENT(IN) :: lastcall ! not used
67
68	! Arguments added for condensation
69	logical,intent(in) :: lecrit ! not used.
70	REAL,INTENT(IN) :: pz0
71
72	! Tracers
73	! --------
74	integer,intent(in) :: nq
75	real,intent(in) :: pqsurf(ngrid,nq)
76	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
77	real,intent(out) :: pdqdif(ngrid,nlay,nq)
78	real,intent(out) :: pdqsdif(ngrid,nq)
79
80	! local
81	! -----
82	integer ilev,ig,ilay,nlev
83
84	REAL z4st,zdplanck(ngrid)
85	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
86	REAL zcdv(ngrid),zcdh(ngrid)
87	REAL zcdv_true(ngrid),zcdh_true(ngrid)
88	REAL zu(ngrid,nlay),zv(ngrid,nlay)
89	REAL zh(ngrid,nlay),zt(ngrid,nlay)
90	REAL ztsrf(ngrid)
91	REAL z1(ngrid),z2(ngrid)
92	REAL zmass(ngrid,nlay)
93	REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
94	REAL zb0(ngrid,nlay)
95	REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
96	REAL zcv(ngrid,nlay),zdv(ngrid,nlay) !inversion coefficient for winds
97	REAL zct(ngrid,nlay),zdt(ngrid,nlay) !inversion coefficient for temperature
98	REAL zcq(ngrid,nlay),zdq(ngrid,nlay) !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,SAVE :: firstcall=.true.
105	!$OMP THREADPRIVATE(firstcall)
106
107	! Tracers
108	! -------
109	INTEGER iq
110	REAL zq(ngrid,nlay,nq)
111	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
112	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
113	REAL zdmassevap(ngrid)
114	REAL rho(ngrid) ! near-surface air density
115	REAL qsat(ngrid),psat(ngrid)
116	REAL kmixmin
117
118	! Variables added for implicit latent heat inclusion
119	! --------------------------------------------------
120	real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
121
122	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
123	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
124
125	real, parameter :: karman=0.4
126	real cd0, roughratio
127
128	real dqsdif_total(ngrid)
129	real zq0(ngrid)
130
131
132	! Coherence test
133	! --------------
134
135	IF (firstcall) THEN
136
137	if(water)then
138	ivap=igcm_h2o_vap
139	iliq=igcm_h2o_ice
140	iliq_surf=igcm_h2o_vap
141	iice_surf=igcm_h2o_ice ! simply to make the code legible
142	! to be generalised
143	else
144	ivap=0
145	iliq=0
146	iliq_surf=0
147	iice_surf=0 ! simply to make the code legible
148	endif
149	sensibFlux(:)=0.
150
151	firstcall=.false.
152	ENDIF
153
154	!-----------------------------------------------------------------------
155	! 1. Initialisation
156	! -----------------
157
158	nlev=nlay+1
159
160	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
161	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
162	! ---------------------------------------------
163
164	DO ilay=1,nlay
165	DO ig=1,ngrid
166	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
167	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
168	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
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)/(pt(ig,ilev-1)+pt(ig,ilev))
176	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
177	ENDDO
178	ENDDO
179	DO ig=1,ngrid
180	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
181	ENDDO
182	dqsdif_total(:)=0.0
183
184	!-----------------------------------------------------------------------
185	! 2. Add the physical tendencies computed so far
186	! ----------------------------------------------
187
188	DO ilev=1,nlay
189	DO ig=1,ngrid
190	zu(ig,ilev)=pu(ig,ilev)
191	zv(ig,ilev)=pv(ig,ilev)
192	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
193	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
194	ENDDO
195	ENDDO
196	if(tracer) then
197	DO iq =1, nq
198	DO ilev=1,nlay
199	DO ig=1,ngrid
200	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
201	ENDDO
202	ENDDO
203	ENDDO
204	if (water) then
205	DO ilev=1,nlay
206	DO ig=1,ngrid
207	zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
208	ENDDO
209	ENDDO
210	Endif
211	end if
212
213	!-----------------------------------------------------------------------
214	! 3. Turbulence scheme
215	! --------------------
216	!
217	! Source of turbulent kinetic energy at the surface
218	! -------------------------------------------------
219	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
220
221	DO ig=1,ngrid
222	roughratio = 1. + pzlay(ig,1)/pz0
223	cd0 = karman/log(roughratio)
224	cd0 = cd0*cd0
225	zcdv_true(ig) = cd0
226	zcdh_true(ig) = cd0
227	if(nosurf)then
228	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
229	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
230	endif
231	ENDDO
232
233	DO ig=1,ngrid
234	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
235	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
236	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
237	ENDDO
238
239	! Turbulent diffusion coefficients in the boundary layer
240	! ------------------------------------------------------
241
242	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
243
244	! Adding eddy mixing to mimic 3D general circulation in 1D
245	! R. Wordsworth & F. Forget (2010)
246	if ((ngrid.eq.1)) then
247	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
248	do ilev=1,nlay
249	do ig=1,ngrid
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	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
257	DO ig=1,ngrid
258	zkv(ig,1)=zcdv(ig)
259	ENDDO
260	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
261
262	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
263	zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
264
265	!-----------------------------------------------------------------------
266	! 4. Implicit inversion of u
267	! --------------------------
268
269	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
270	! avec
271	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
272	! et
273	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
274	! donc les entrees sont /zcdv/ pour la condition a la limite sol
275	! et /zkv/ = Ku
276
277	DO ig=1,ngrid
278	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
279	zcv(ig,nlay)=zmass(ig,nlay)zu(ig,nlay)z1(ig)
280	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
281	ENDDO
282
283	DO ilay=nlay-1,1,-1
284	DO ig=1,ngrid
285	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
286	zcv(ig,ilay)=(zmass(ig,ilay)zu(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
287	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
288	ENDDO
289	ENDDO
290
291	DO ig=1,ngrid
292	zu(ig,1)=zcv(ig,1)
293	ENDDO
294	DO ilay=2,nlay
295	DO ig=1,ngrid
296	zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1)
297	ENDDO
298	ENDDO
299
300	!-----------------------------------------------------------------------
301	! 5. Implicit inversion of v
302	! --------------------------
303
304	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
305	! avec
306	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
307	! et
308	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
309	! donc les entrees sont /zcdv/ pour la condition a la limite sol
310	! et /zkv/ = Kv
311
312	DO ig=1,ngrid
313	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
314	zcv(ig,nlay)=zmass(ig,nlay)zv(ig,nlay)z1(ig)
315	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
316	ENDDO
317
318	DO ilay=nlay-1,1,-1
319	DO ig=1,ngrid
320	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
321	zcv(ig,ilay)=(zmass(ig,ilay)zv(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
322	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
323	ENDDO
324	ENDDO
325
326	DO ig=1,ngrid
327	zv(ig,1)=zcv(ig,1)
328	ENDDO
329	DO ilay=2,nlay
330	DO ig=1,ngrid
331	zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1)
332	ENDDO
333	ENDDO
334
335	! Calcul of wind stress
336
337	DO ig=1,ngrid
338	flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
339	flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
340	ENDDO
341
342
343	!----------------------------------------------------------------------------
344	! 6. Implicit inversion of h, not forgetting the coupling with the ground
345
346	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
347	! avec
348	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
349	! et
350	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
351	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
352	! et /zkh/ = Kh
353
354	! Using the wind modified by friction for lifting and sublimation
355	! ---------------------------------------------------------------
356	DO ig=1,ngrid
357	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
358	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
359	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
360	zkh(ig,1)= zcdh(ig)
361	ENDDO
362
363	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
364	! ---------------------------------------------------------------
365	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
366	zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
367
368	DO ig=1,ngrid
369	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
370	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
371	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
372	ENDDO
373
374	DO ilay=nlay-1,2,-1
375	DO ig=1,ngrid
376	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
377	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
378	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
379	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
380	ENDDO
381	ENDDO
382
383	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
384	DO ig=1,ngrid
385	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
386	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
387	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
388	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
389	ENDDO
390
391
392	! Calculate (d Planck / dT) at the interface temperature
393	! ------------------------------------------------------
394
395	z4st=4.0sigmaptimestep
396	DO ig=1,ngrid
397	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
398	ENDDO
399
400	! Calculate temperature tendency at the interface (dry case)
401	! ----------------------------------------------------------
402	! Sum of fluxes at interface at time t + \delta t gives change in T:
403	! radiative fluxes
404	! turbulent convective (sensible) heat flux
405	! flux (if any) from subsurface
406
407	if(.not.water) then
408
409	DO ig=1,ngrid
410	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
411	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
412	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
413	ztsrf(ig) = z1(ig) / z2(ig)
414	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
415	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
416	ENDDO
417	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
418
419
420	! Recalculate temperature to top of atmosphere, starting from ground
421	! ------------------------------------------------------------------
422
423	DO ilay=2,nlay
424	DO ig=1,ngrid
425	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
426	ENDDO
427	ENDDO
428
429	endif ! not water
430
431	!-----------------------------------------------------------------------
432	! TRACERS (no vapour)
433	! -------
434
435	if(tracer) then
436
437	! Calculate vertical flux from the bottom to the first layer (dust)
438	! -----------------------------------------------------------------
439	do ig=1,ngrid
440	rho(ig) = zb0(ig,1) /ptimestep
441	end do
442
443	pdqsdif(:,:)=0.
444
445	! Implicit inversion of q
446	! -----------------------
447	do iq=1,nq
448
449	if (iq.ne.ivap) then
450
451	DO ig=1,ngrid
452	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
453	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
454	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
455	ENDDO
456
457	DO ilay=nlay-1,2,-1
458	DO ig=1,ngrid
459	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
460	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
461	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
462	ENDDO
463	ENDDO
464
465	if ((water).and.(iq.eq.iliq)) then
466	! special case for condensed water tracer: do not include
467	! h2o ice tracer from surface (which is set when handling
468	! h2o vapour case (see further down).
469	! zb(ig,1)=0 if iq ne ivap
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))*z1(ig)
473	ENDDO
474	else ! general case
475	do ig=1,ngrid
476	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
477	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
478	! tracer flux from surface
479	! currently pdqsdif always zero here,
480	! so last line is superfluous
481	enddo
482	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
483
484
485	! Starting upward calculations for simple tracer mixing (e.g., dust)
486	do ig=1,ngrid
487	zq(ig,1,iq)=zcq(ig,1)
488	end do
489
490	do ilay=2,nlay
491	do ig=1,ngrid
492	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
493	end do
494	end do
495
496	endif ! if (iq.ne.ivap)
497
498	! Calculate temperature tendency including latent heat term
499	! and assuming an infinite source of water on the ground
500	! ------------------------------------------------------------------
501
502	if (water.and.(iq.eq.ivap)) then
503
504	! compute evaporation efficiency
505	do ig=1,ngrid
506	if(nint(rnat(ig)).eq.1)then
507	dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf)
508	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
509	dryness(ig)=MAX(0.,dryness(ig))
510	endif
511	enddo
512
513	do ig=1,ngrid
514	! Calculate the value of qsat at the surface (water)
515	call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
516	call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
517	call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
518	dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
519	! calculate dQsat / dT by finite differences
520	! we cannot use the updated temperature value yet...
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	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
722	if (tracer) then
723	call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
724	endif
725	endif
726
727	! if(lastcall)then
728	! if(ngrid.eq.1)then
729	! print*,'Saving k.out...'
730	! OPEN(12,file='k.out',form='formatted')
731	! DO ilay=1,nlay
732	! write(12,*) zkh(1,ilay), pplay(1,ilay)
733	! ENDDO
734	! CLOSE(12)
735	! endif
736	! endif
737
738	end

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

Download in other formats:

Original Format