Context Navigation

vdifc.F @ 787

Last change on this file since 787 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

File size: 24.0 KB

Line
1	subroutine vdifc(ngrid,nlay,nq,rnat,ppopsk,
2	& ptimestep,pcapcal,lecrit,
3	& pplay,pplev,pzlay,pzlev,pz0,
4	& pu,pv,ph,pq,ptsrf,pemis,pqsurf,
5	& pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
6	& pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
7	& pdqdif,pdqsdif,lastcall)
8
9	use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
10	use radcommon_h, only : sigma
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	!
34	!==================================================================
35
36	!-----------------------------------------------------------------------
37	! declarations
38	! ------------
39
40	#include "dimensions.h"
41	#include "dimphys.h"
42	#include "comcstfi.h"
43	#include "callkeys.h"
44
45	! arguments
46	! ---------
47	INTEGER ngrid,nlay
48	REAL ptimestep
49	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
50	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
51	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
52	REAL ptsrf(ngrid),pemis(ngrid)
53	REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
54	REAL pfluxsrf(ngrid)
55	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
56	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
57	REAL pq2(ngrid,nlay+1)
58
59	integer rnat(ngrid)
60
61	! Arguments added for condensation
62	REAL ppopsk(ngrid,nlay)
63	logical lecrit
64	REAL pz0
65
66	! Tracers
67	! --------
68	integer nq
69	real pqsurf(ngrid,nq)
70	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
71	real pdqdif(ngrid,nlay,nq)
72	real pdqsdif(ngrid,nq)
73
74	! local
75	! -----
76	integer ilev,ig,ilay,nlev
77
78	REAL z4st,zdplanck(ngrid)
79	REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
80	REAL zcdv(ngrid),zcdh(ngrid)
81	REAL zcdv_true(ngrid),zcdh_true(ngrid)
82	REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
83	REAL zh(ngrid,nlayermx)
84	REAL ztsrf2(ngrid)
85	REAL z1(ngrid),z2(ngrid)
86	REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
87	REAL zb0(ngrid,nlayermx)
88	REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
89	REAL zcst1
90	REAL zu2!, a
91	REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)
92	REAL evap(ngrid)
93	REAL zcq0(ngrid),zdq0(ngrid)
94	REAL zx_alf1(ngrid),zx_alf2(ngrid)
95
96	LOGICAL firstcall
97	SAVE firstcall
98
99	LOGICAL lastcall
100
101	! variables added for CO2 condensation
102	! ------------------------------------
103	REAL hh !, zhcond(ngrid,nlayermx)
104	! REAL latcond,tcond1mb
105	! REAL acond,bcond
106	! SAVE acond,bcond
107	! DATA latcond,tcond1mb/5.9e5,136.27/
108
109	! Tracers
110	! -------
111	INTEGER iq
112	REAL zq(ngrid,nlayermx,nq)
113	REAL zq1temp(ngrid)
114	REAL rho(ngrid) ! near-surface air density
115	REAL qsat(ngrid)
116	DATA firstcall/.true./
117	REAL kmixmin
118
119	! Variables added for implicit latent heat inclusion
120	! --------------------------------------------------
121	real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
122	real z1_Tdry(ngrid), z2_Tdry(ngrid)
123	real z1_T(ngrid), z2_T(ngrid)
124	real zb_T(ngrid)
125	real zc_T(ngrid,nlayermx)
126	real zd_T(ngrid,nlayermx)
127	real lat1(ngrid), lat2(ngrid)
128	real surfh2otot
129	logical surffluxdiag
130	integer isub ! sub-loop for precision
131
132	integer ivap, iice ! also make liq for clarity on surface...
133	save ivap, iice
134
135	real, parameter :: karman=0.4
136	real cd0, roughratio
137
138	logical forceWC
139	real masse, Wtot, Wdiff
140
141	real dqsdif_total(ngrid)
142	real zq0(ngrid)
143
144	forceWC=.true.
145	! forceWC=.false.
146
147
148	! Coherence test
149	! --------------
150
151	IF (firstcall) THEN
152	! To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
153	! bcond=1./tcond1mb
154	! acond=r/latcond
155	! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
156	! PRINT*,' acond,bcond',acond,bcond
157
158	if(water)then
159	! iliq=igcm_h2o_vap
160	ivap=igcm_h2o_vap
161	iice=igcm_h2o_ice ! simply to make the code legible
162	! to be generalised later
163	endif
164
165	firstcall=.false.
166	ENDIF
167
168	!-----------------------------------------------------------------------
169	! 1. Initialisation
170	! -----------------
171
172	nlev=nlay+1
173
174	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
175	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
176	! ---------------------------------------------
177
178	DO ilay=1,nlay
179	DO ig=1,ngrid
180	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
181	ENDDO
182	ENDDO
183
184	zcst1=4.gptimestep/(R*R)
185	DO ilev=2,nlev-1
186	DO ig=1,ngrid
187	zb0(ig,ilev)=pplev(ig,ilev)*
188	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
189	s (ph(ig,ilev-1)+ph(ig,ilev))
190	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
191	s (pplay(ig,ilev-1)-pplay(ig,ilev))
192	ENDDO
193	ENDDO
194	DO ig=1,ngrid
195	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
196	ENDDO
197
198	dqsdif_total(:)=0.0
199
200	!-----------------------------------------------------------------------
201	! 2. Add the physical tendencies computed so far
202	! ----------------------------------------------
203
204	DO ilev=1,nlay
205	DO ig=1,ngrid
206	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
207	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
208	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
209	ENDDO
210	ENDDO
211	if(tracer) then
212	DO iq =1, nq
213	DO ilev=1,nlay
214	DO ig=1,ngrid
215	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
216	& pdqfi(ig,ilev,iq)*ptimestep
217	ENDDO
218	ENDDO
219	ENDDO
220	end if
221
222	!-----------------------------------------------------------------------
223	! 3. Turbulence scheme
224	! --------------------
225	!
226	! Source of turbulent kinetic energy at the surface
227	! -------------------------------------------------
228	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
229
230	DO ig=1,ngrid
231	roughratio = 1.E+0 + pzlay(ig,1)/pz0
232	cd0 = karman/log(roughratio)
233	cd0 = cd0*cd0
234	zcdv_true(ig) = cd0
235	zcdh_true(ig) = cd0
236	ENDDO
237
238	DO ig=1,ngrid
239	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
240	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
241	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
242	ENDDO
243
244	! Turbulent diffusion coefficients in the boundary layer
245	! ------------------------------------------------------
246
247	call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay
248	& ,pu,pv,ph,zcdv_true
249	& ,pq2,zkv,zkh)
250
251	! Adding eddy mixing to mimic 3D general circulation in 1D
252	! R. Wordsworth & F. Forget (2010)
253	if ((ngrid.eq.1)) then
254	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
255	do ilev=1,nlay
256	do ig=1,ngrid
257	!zkh(ig,ilev) = 1.0
258	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
259	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
260	end do
261	end do
262	end if
263
264	!-----------------------------------------------------------------------
265	! 4. Implicit inversion of u
266	! --------------------------
267
268	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
269	! avec
270	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
271	! et
272	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
273	! donc les entrees sont /zcdv/ pour la condition a la limite sol
274	! et /zkv/ = Ku
275
276	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
277	CALL multipl(ngrid,zcdv,zb0,zb)
278
279	DO ig=1,ngrid
280	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
281	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
282	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
283	ENDDO
284
285	DO ilay=nlay-1,1,-1
286	DO ig=1,ngrid
287	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
288	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
289	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
290	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
291	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
292	ENDDO
293	ENDDO
294
295	DO ig=1,ngrid
296	zu(ig,1)=zc(ig,1)
297	ENDDO
298	DO ilay=2,nlay
299	DO ig=1,ngrid
300	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
301	ENDDO
302	ENDDO
303
304	!-----------------------------------------------------------------------
305	! 5. Implicit inversion of v
306	! --------------------------
307
308	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
309	! avec
310	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
311	! et
312	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
313	! donc les entrees sont /zcdv/ pour la condition a la limite sol
314	! et /zkv/ = Kv
315
316	DO ig=1,ngrid
317	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
318	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
319	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
320	ENDDO
321
322	DO ilay=nlay-1,1,-1
323	DO ig=1,ngrid
324	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
325	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
326	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
327	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
328	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
329	ENDDO
330	ENDDO
331
332	DO ig=1,ngrid
333	zv(ig,1)=zc(ig,1)
334	ENDDO
335	DO ilay=2,nlay
336	DO ig=1,ngrid
337	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
338	ENDDO
339	ENDDO
340
341	!----------------------------------------------------------------------------
342	! 6. Implicit inversion of h, not forgetting the coupling with the ground
343
344	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
345	! avec
346	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
347	! et
348	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
349	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
350	! et /zkh/ = Kh
351
352	! Using the wind modified by friction for lifting and sublimation
353	! ---------------------------------------------------------------
354	DO ig=1,ngrid
355	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
356	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
357	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
358	ENDDO
359
360	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
361	CALL multipl(ngrid,zcdh,zb0,zb)
362
363	DO ig=1,ngrid
364	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
365	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
366	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
367	ENDDO
368
369	DO ilay=nlay-1,2,-1
370	DO ig=1,ngrid
371	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
372	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
373	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
374	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
375	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
376	ENDDO
377	ENDDO
378
379	DO ig=1,ngrid
380	z1(ig)=1./(za(ig,1)+zb(ig,1)+
381	& zb(ig,2)*(1.-zd(ig,2)))
382	zc(ig,1)=(za(ig,1)*zh(ig,1)+
383	& zb(ig,2)zc(ig,2))z1(ig)
384	zd(ig,1)=zb(ig,1)*z1(ig)
385	ENDDO
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
406	z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
407	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
408	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
409	& +zdplanck(ig)
410	ztsrf2(ig) = z1(ig) / z2(ig)
411	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
412	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
413	ENDDO
414
415	! Recalculate temperature to top of atmosphere, starting from ground
416	! ------------------------------------------------------------------
417
418	DO ilay=2,nlay
419	DO ig=1,ngrid
420	hh = zh(ig,ilay-1)
421	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
422	ENDDO
423	ENDDO
424
425	endif ! not water
426
427	!-----------------------------------------------------------------------
428	! TRACERS (no vapour)
429	! -------
430
431	if(tracer) then
432
433	! Calculate vertical flux from the bottom to the first layer (dust)
434	! -----------------------------------------------------------------
435	do ig=1,ngrid
436	rho(ig) = zb0(ig,1) /ptimestep
437	end do
438
439	call zerophys(ngrid*nq,pdqsdif)
440
441	! Implicit inversion of q
442	! -----------------------
443	do iq=1,nq
444
445	if (iq.ne.igcm_h2o_vap) then
446
447	DO ig=1,ngrid
448	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
449	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
450	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
451	ENDDO
452
453	DO ilay=nlay-1,2,-1
454	DO ig=1,ngrid
455	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
456	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
457	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
458	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
459	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
460	ENDDO
461	ENDDO
462
463	if ((water).and.(iq.eq.iice)) then
464	! special case for water ice tracer: do not include
465	! h2o ice tracer from surface (which is set when handling
466	! h2o vapour case (see further down).
467	! zb(ig,1)=0 if iq ne ivap
468	DO ig=1,ngrid
469	z1(ig)=1./(za(ig,1)+
470	& zb(ig,2)*(1.-zdq(ig,2)))
471	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
472	& zb(ig,2)zcq(ig,2))z1(ig)
473	ENDDO
474	else ! general case
475	DO ig=1,ngrid
476	z1(ig)=1./(za(ig,1)+
477	& zb(ig,2)*(1.-zdq(ig,2)))
478	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
479	& zb(ig,2)*zcq(ig,2)
480	& +(-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)+
496	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
497	end do
498	end do
499
500	endif ! if (iq.ne.igcm_h2o_vap)
501
502	! Calculate temperature tendency including latent heat term
503	! and assuming an infinite source of water on the ground
504	! ------------------------------------------------------------------
505
506	if (water.and.(iq.eq.igcm_h2o_vap)) then
507
508	! compute evaporation efficiency
509	do ig = 1, ngrid
510	if(rnat(ig).eq.1)then
511	dryness(ig)=pqsurf(ig,ivap)+pqsurf(ig,iice)
512	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
513	dryness(ig)=MAX(0.,dryness(ig))
514	endif
515	enddo
516
517	do ig=1,ngrid
518
519	! Calculate the value of qsat at the surface (water)
520	call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
521	call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
522	call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
523	dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
524	! calculate dQsat / dT by finite differences
525	! we cannot use the updated temperature value yet...
526
527	enddo
528
529	! coefficients for q
530
531	do ig=1,ngrid
532	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
533	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
534	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
535	enddo
536
537	do ilay=nlay-1,2,-1
538	do ig=1,ngrid
539	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
540	$ zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
541	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
542	$ zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
543	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
544	enddo
545	enddo
546
547	do ig=1,ngrid
548	z1(ig)=1./(za(ig,1)+zb(ig,1)*dryness(ig)+
549	$ zb(ig,2)*(1.-zdq(ig,2)))
550	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
551	$ zb(ig,2)zcq(ig,2))z1(ig)
552	zdq(ig,1)=dryness(ig)zb(ig,1)z1(ig)
553	enddo
554
555	! calculation of h0 and h1
556	do ig=1,ngrid
557	zdq0(ig) = dqsat(ig)
558	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
559
560	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzb(ig,1)*zc(ig,1)
561	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
562	& + zb(ig,1)dryness(ig)RLVTT*
563	& ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
564
565	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
566	& +zdplanck(ig)
567	& +zb(ig,1)dryness(ig)RLVTTzdq0(ig)
568	& (1.0-zdq(ig,1))
569
570	ztsrf2(ig) = z1(ig) / z2(ig)
571	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
572	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
573	enddo
574
575	! calculation of qs and q1
576	do ig=1,ngrid
577	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf2(ig)
578	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
579	enddo
580
581	! calculation of evaporation
582	do ig=1,ngrid
583	evap(ig)= zb(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
584	dqsdif_total(ig)=evap(ig)
585	enddo
586
587	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
588	do ilay=2,nlay
589	do ig=1,ngrid
590	zq(ig,ilay,iq)=zcq(ig,ilay)
591	& +zdq(ig,ilay)*zq(ig,ilay-1,iq)
592	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
593	end do
594	end do
595
596	do ig=1,ngrid
597
598	! --------------------------------------------------------------------------
599	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
600	! The surface vapour tracer is actually liquid. To make things difficult.
601
602	if (rnat(ig).eq.0) then ! unfrozen ocean
603
604	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
605	pdqsdif(ig,iice)=0.0
606
607
608	elseif (rnat(ig).eq.1) then ! (continent)
609
610	! --------------------------------------------------------
611	! Now check if we've taken too much water from the surface
612	! This can only occur on the continent
613
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	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
618	pdqsdif(ig,iice)=max(-pqsurf(ig,iice)/ptimestep
619	& ,pdqsdif(ig,iice))
620	endif
621	! sublimation only greater than qsurf(ice)
622	! ----------------------------------------
623	! we just convert some liquid to vapour too
624	! if latent heats are the same, no big deal
625	if (-dqsdif_total(ig).gt.pqsurf(ig,iice))then
626	pdqsdif(ig,iice) = -pqsurf(ig,iice)/ptimestep ! removes all the ice!
627	pdqsdif(ig,ivap) = dqsdif_total(ig)/ptimestep
628	& - pdqsdif(ig,iice) ! take the remainder from the liquid instead
629	pdqsdif(ig,ivap) = max(-pqsurf(ig,ivap)/ptimestep
630	& ,pdqsdif(ig,ivap))
631	endif
632
633	endif ! if (rnat.ne.1)
634
635	! If water vapour is condensing, we must decide whether it forms ice or liquid.
636	if(dqsdif_total(ig).gt.0)then ! a bug was here!
637	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
638	pdqsdif(ig,iice)=0.0
639	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
640	else
641	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
642	pdqsdif(ig,ivap)=0.0
643	endif
644	endif
645
646	end do ! of DO ig=1,ngrid
647	endif ! if (water et iq=ivap)
648	end do ! of do iq=1,nq
649	endif ! traceur
650
651
652	!-----------------------------------------------------------------------
653	! 8. Final calculation of the vertical diffusion tendencies
654	! -----------------------------------------------------------------
655
656	do ilev = 1, nlay
657	do ig=1,ngrid
658	pdudif(ig,ilev)=(zu(ig,ilev)-
659	& (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
660	pdvdif(ig,ilev)=(zv(ig,ilev)-
661	& (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
662	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
663
664	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
665	enddo
666	enddo
667
668	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
669	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
670	ENDDO
671
672	if (tracer) then
673	do iq = 1, nq
674	do ilev = 1, nlay
675	do ig=1,ngrid
676	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
677	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
678	& ptimestep
679	enddo
680	enddo
681	enddo
682
683	if(water.and.forceWC)then ! force water conservation in model
684	! we calculate the difference and add it to the ground
685	! this is ugly and should be improved in the future
686	do ig=1,ngrid
687	Wtot=0.0
688	do ilay=1,nlay
689	masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
690	! Wtot=Wtot+masse*(zq(ig,ilay,iice)-
691	! & (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep))
692	Wtot=Wtot+masse*(zq(ig,ilay,ivap)-
693	& (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep))
694	enddo
695	Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
696
697	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
698	pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
699	else
700	pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff
701	endif
702	enddo
703
704	endif
705
706	endif
707
708	if(water)then
709	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
710	endif
711
712	! if(lastcall)then
713	! if(ngrid.eq.1)then
714	! print*,'Saving k.out...'
715	! OPEN(12,file='k.out',form='formatted')
716	! DO ilay=1,nlay
717	! write(12,*) zkh(1,ilay), pplay(1,ilay)
718	! ENDDO
719	! CLOSE(12)
720	! endif
721	! endif
722
723
724	return
725	end

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

Download in other formats:

Original Format