Context Navigation

vdifc_mod.F @ 3167

Last change on this file since 3167 was 3167, checked in by llange, 11 months ago
Mars PCM Introducing the scheme from ATKE workshop to solve turbulent diffusion + surface layer parameterization. Works only if callatke = .true. in the callphys.def. By default, it is false and the model runs as usual with the yamada 2.5 closure Tuning of the several parameters for the ATKE in progress LL
File size: 58.8 KB

Line
1	MODULE vdifc_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE vdifc(ngrid,nlay,nsoil,nq,nqsoil,ppopsk,
8	$ ptimestep,pcapcal,lecrit,
9	$ pplay,pplev,pzlay,pzlev,pz0,
10	$ pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
11	$ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
12	$ pdudif,pdvdif,pdhdif,pdtsrf,pq2,
13	$ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
14	$ hfmax,pcondicea_co2microp,sensibFlux,
15	$ dustliftday,local_time,watercap, dwatercap_dif)
16
17	use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
18	& igcm_dust_submicron, igcm_h2o_vap,
19	& igcm_h2o_ice, alpha_lift, igcm_co2,
20	& igcm_hdo_vap, igcm_hdo_ice,
21	& igcm_stormdust_mass, igcm_stormdust_number
22	use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness,
23	& old_wsublimation_scheme
24	USE comcstfi_h, ONLY: cpp, r, rcp, g, pi
25	use watersat_mod, only: watersat
26	use turb_mod, only: turb_resolved, ustar, tstar
27	use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
28	use hdo_surfex_mod, only: hdo_surfex
29	c use geometry_mod, only: longitude_deg,latitude_deg ! Joseph
30	use dust_param_mod, only: doubleq, submicron, lifting
31	use write_output_mod, only: write_output
32	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
33	& subslope_dist,major_slope,iflat
34	use microphys_h, only: To
35	use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
36	use comsoil_h, only: layer, mlayer,adsorption_soil
37	use vdif_cd_mod, only: vdif_cd
38	use lmdz_call_atke, only: call_atke
39	IMPLICIT NONE
40
41	c=======================================================================
42	c
43	c subject:
44	c --------
45	c Turbulent diffusion (mixing) for potential T, U, V and tracer
46	c
47	c Shema implicite
48	c On commence par rajouter au variables x la tendance physique
49	c et on resoult en fait:
50	c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
51	c
52	c author:
53	c ------
54	c Hourdin/Forget/Fournier
55	c=======================================================================
56
57	c-----------------------------------------------------------------------
58	c declarations:
59	c -------------
60
61	include "callkeys.h"
62
63	c
64	c arguments:
65	c ----------
66
67	INTEGER,INTENT(IN) :: ngrid,nlay,nsoil,nqsoil
68	REAL,INTENT(IN) :: ptimestep
69	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
70	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
71	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
72	REAL,INTENT(IN) :: ph(ngrid,nlay)
73	REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope)
74	REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
75	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
76	REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope)
77	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
78	REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay)
79	REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
80	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
81	REAL,INTENT(IN) :: ptsoil(ngrid,nsoil,nslope)
82
83	c Argument added for condensation:
84	REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
85	logical,INTENT(IN) :: lecrit
86	REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1)
87
88	REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m)
89
90	c Argument added to account for subgrid gustiness :
91
92	REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
93
94	c Traceurs :
95	integer,intent(in) :: nq
96	REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope)
97	REAL :: zqsurf(ngrid) ! temporary water tracer
98	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
99	real,intent(out) :: pdqdif(ngrid,nlay,nq)
100	real,intent(out) :: pdqsdif(ngrid,nq,nslope)
101	REAL,INTENT(in) :: dustliftday(ngrid)
102	REAL,INTENT(inout) :: qsoil(ngrid,nsoil,nqsoil,nslope) !subsurface tracers
103	REAL,INTENT(in) :: local_time(ngrid)
104
105	c local:
106	c ------
107
108	REAL :: pt(ngrid,nlay)
109
110	INTEGER ilev,ig,ilay,nlev,islope,ik,lice
111
112	REAL z4st,zdplanck(ngrid)
113	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
114	REAL zkq(ngrid,nlay+1)
115	REAL zcdv(ngrid,nslope),zcdh(ngrid,nslope)
116	REAL, INTENT(OUT) :: zcdv_true(ngrid,nslope)
117	REAL, INTENT(OUT) :: zcdh_true(ngrid,nslope) ! drag coeff are used by the LES to recompute u* and hfx
118	REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid) ! drag coeffs for the major sub-grid surface
119	REAL :: zcdv_true_tmp(ngrid),zcdh_true_tmp(ngrid) ! drag coeffs (computed with wind gustiness for the major sub-grid surface
120	REAL zu(ngrid,nlay),zv(ngrid,nlay)
121	REAL zh(ngrid,nlay)
122	REAL ztsrf2(ngrid)
123	REAL z1(ngrid),z2(ngrid)
124	REAL za(ngrid,nlay),zb(ngrid,nlay)
125	REAL zb0(ngrid,nlay)
126	REAL zc(ngrid,nlay),zd(ngrid,nlay)
127	REAL zcst1
128	REAL zu2(ngrid)
129	REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located.
130	REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface
131	REAL dist_up(ngrid,nslope) !distance from ice to layer above
132	REAL dist_down(ngrid,nslope) !distance from ice to layer down
133	REAL dist_sum(ngrid,nslope) ! sum of distance
134	REAL zdqsdif_ssi(ngrid,nslope) !SSI flux
135	REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction
136
137	EXTERNAL SSUM,SCOPY
138	REAL SSUM
139	LOGICAL,SAVE :: firstcall=.true.
140
141	!$OMP THREADPRIVATE(firstcall)
142
143	c variable added for CO2 condensation:
144	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
145	REAL hh , zhcond(ngrid,nlay)
146	REAL,PARAMETER :: latcond=5.9e5
147	REAL,PARAMETER :: tcond1mb=136.27
148	REAL,SAVE :: acond,bcond
149
150	!$OMP THREADPRIVATE(acond,bcond)
151
152	c Subtimestep & implicit treatment of water vapor
153	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
154	REAL zdqsdif_surf(ngrid) ! subtimestep pdqsdif for water ice
155	REAL ztsrf(ngrid) ! temporary surface temperature in tsub
156	REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub
157	REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux
158	REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux
159
160	c For latent heat release from ground water ice sublimation
161	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
162	REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect
163	REAL lh ! latent heat, formulation given in the Technical Document:
164	! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004
165
166	c Tracers :
167	c ~~~~~~~
168	INTEGER iq
169	REAL zq(ngrid,nlay,nq)
170	REAL zq1temp(ngrid)
171	REAL rho(ngrid) ! near surface air density
172	REAL qsat(ngrid)
173	REAL qsat2(ngrid,nslope)
174	REAL resist(ngrid,nslope) !subsurface ice flux reduction coef
175
176	REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO
177	REAL hdoflux_meshavg(ngrid) ! value of vapour flux of HDO
178	! REAL h2oflux(ngrid) ! value of vapour flux of H2O
179	REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement
180	REAL saved_h2o_vap(ngrid) ! traceur d'eau avant traitement
181
182	REAL kmixmin
183
184	c Argument added for surface water ice budget:
185	REAL,INTENT(IN) :: watercap(ngrid,nslope)
186	REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope)
187
188	c Subtimestep to compute h2o latent heat flux:
189	REAL :: dtmax = 0.5 ! subtimestep temp criterion
190	INTEGER tsub ! adaptative subtimestep (seconds)
191	REAL subtimestep !ptimestep/nsubtimestep
192	INTEGER nsubtimestep(ngrid) ! number of subtimestep (int)
193
194	c Mass-variation scheme :
195	c ~~~~~~~
196
197	INTEGER j,l
198	REAL zcondicea(ngrid,nlay)
199	REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1)
200	REAL betam(ngrid,nlay),dmice(ngrid,nlay)
201	REAL pdtc(ngrid,nlay)
202	REAL zhs(ngrid,nlay)
203	REAL,SAVE :: ccond
204
205	!$OMP THREADPRIVATE(ccond)
206
207	c Theta_m formulation for mass-variation scheme :
208	c ~~~~~~~
209
210	INTEGER,SAVE :: ico2
211	INTEGER llnt(ngrid)
212	REAL,SAVE :: m_co2, m_noco2, A , B
213	REAL vmr_co2(ngrid,nlay)
214	REAL qco2,mmean(ngrid,nlay)
215
216	!$OMP THREADPRIVATE(ico2,m_co2,m_noco2,A,B)
217
218	REAL,INTENT(OUT) :: sensibFlux(ngrid)
219
220	!!MARGAUX
221	REAL DoH_vap(ngrid,nlay)
222	!! Sub-grid scale slopes
223	REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting
224	REAL :: watercap_tmp(ngrid)
225	REAL :: zq_slope_vap(ngrid,nlay,nq,nslope)
226	REAL :: zq_tmp_vap(ngrid,nlay,nq)
227	REAL :: ptsrf_tmp(ngrid)
228	REAL :: pqsurf_tmp(ngrid,nq)
229	REAL :: pdqsdif_tmphdo(ngrid,nq)
230	REAL :: qsat_tmp(ngrid)
231	INTEGER :: indmax
232
233	character*2 str2
234
235	!! Subsurface exchanges
236	LOGICAL :: exchange ! boolean to check if exchange between the subsurface and the atmosphere can occurs
237	REAL :: zdqsdif_regolith(ngrid,nslope) ! Flux from subsurface (positive pointing outwards) (kg/m^2/s)
238	REAL zq1temp_regolith(ngrid) ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
239	REAL zdqsdif_tot(ngrid) ! subtimestep pdqsdif for water ice
240	LOGICAL :: writeoutput ! boolean to say to soilexchange.F if we are at the last iteration and thus if he can write in the diagsoil
241
242	!! Water buyoncy
243	LOGICAL :: virtual
244
245
246
247	c ** un petit test de coherence
248	c --------------------------
249
250	! AS: OK firstcall absolute
251	IF (firstcall) THEN
252	c To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
253	bcond=1./tcond1mb
254	acond=r/latcond
255	ccond=cpp/(g*latcond)
256	PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
257	PRINT*,' acond,bcond,ccond',acond,bcond,ccond
258
259
260	ico2=0
261
262	c Prepare Special treatment if one of the tracer is CO2 gas
263	do iq=1,nq
264	if (noms(iq).eq."co2") then
265	ico2=iq
266	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
267	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
268	c Compute A and B coefficient use to compute
269	c mean molecular mass Mair defined by
270	c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
271	c 1/Mair = A*q(ico2) + B
272	A =(1/m_co2 - 1/m_noco2)
273	B=1/m_noco2
274	endif
275	enddo
276
277	firstcall=.false.
278	ENDIF
279
280	DO ig = 1,ngrid
281	ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig))
282	pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig))
283	ENDDO
284
285	c-----------------------------------------------------------------------
286	c 1. initialisation
287	c -----------------
288
289	nlev=nlay+1
290
291	! initialize output tendencies to zero:
292	pdudif(1:ngrid,1:nlay)=0
293	pdvdif(1:ngrid,1:nlay)=0
294	pdhdif(1:ngrid,1:nlay)=0
295	pdtsrf(1:ngrid,1:nslope)=0
296	zdtsrf(1:ngrid,1:nslope)=0
297	surf_h2o_lh(1:ngrid,1:nslope)=0
298	zsurf_h2o_lh(1:ngrid,1:nslope)=0
299	pdqdif(1:ngrid,1:nlay,1:nq)=0
300	pdqsdif(1:ngrid,1:nq,1:nslope)=0
301	pdqsdif_tmp(1:ngrid,1:nq)=0
302	zdqsdif_surf(1:ngrid)=0
303	dwatercap_dif(1:ngrid,1:nslope)=0
304	zdqsdif_regolith(1:ngrid,1:nslope)=0
305	zq1temp_regolith(1:ngrid)=0
306	zdqsdif_tot(1:ngrid)=0
307	h2o_ice_depth(1:ngrid,1:nslope)=1
308	virtual = .false.
309
310	c ** calcul de rhodz et dtrho/dz=dtrho*2 g/dp
311	c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
312	c ----------------------------------------
313
314	DO ilay=1,nlay
315	DO ig=1,ngrid
316	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
317	! Mass variation scheme:
318	betam(ig,ilay)=-za(ig,ilay)latcond/(cppppopsk(ig,ilay))
319	ENDDO
320	ENDDO
321
322	zcst1=4.gptimestep/(r*r)
323	DO ilev=2,nlev-1
324	DO ig=1,ngrid
325	zb0(ig,ilev)=pplev(ig,ilev)*
326	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
327	s (ph(ig,ilev-1)+ph(ig,ilev))
328	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
329	s (pplay(ig,ilev-1)-pplay(ig,ilev))
330	ENDDO
331	ENDDO
332	DO ig=1,ngrid
333	zb0(ig,1)=ptimesteppplev(ig,1)/(rptsrf_tmp(ig))
334	ENDDO
335
336	c ** diagnostique pour l'initialisation
337	c ----------------------------------
338
339	IF(lecrit) THEN
340	ig=ngrid/2+1
341	PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
342	DO ilay=1,nlay
343	WRITE(*,'(6f11.5)')
344	s .01pplay(ig,ilay),.001pzlay(ig,ilay),
345	s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
346	ENDDO
347	PRINT*,'Pression (mbar) ,altitude (km),zb'
348	DO ilev=1,nlay
349	WRITE(*,'(3f15.7)')
350	s .01pplev(ig,ilev),.001pzlev(ig,ilev),
351	s zb0(ig,ilev)
352	ENDDO
353	ENDIF
354
355	c -----------------------------------
356	c Potential Condensation temperature:
357	c -----------------------------------
358
359	c Compute CO2 Volume mixing ratio
360	c -------------------------------
361	if (callcond.and.(ico2.ne.0)) then
362	DO ilev=1,nlay
363	DO ig=1,ngrid
364	qco2=MAX(1.E-30
365	& ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
366	c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
367	mmean(ig,ilev)=1/(A*qco2 +B)
368	vmr_co2(ig,ilev) = qco2*mmean(ig,ilev)/m_co2
369	ENDDO
370	ENDDO
371	else
372	DO ilev=1,nlay
373	DO ig=1,ngrid
374	vmr_co2(ig,ilev)=0.95
375	ENDDO
376	ENDDO
377	end if
378
379	c forecast of atmospheric temperature zt and frost temperature ztcond
380	c --------------------------------------------------------------------
381
382	if (callcond) then
383	DO ilev=1,nlay
384	DO ig=1,ngrid
385	ztcond(ig,ilev)=
386	& 1./(bcond-acondlog(.01vmr_co2(ig,ilev)*pplay(ig,ilev)))
387	if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
388	! zhcond(ig,ilev) =
389	! & (1./(bcond-acondlog(.0095pplay(ig,ilev))))/ppopsk(ig,ilev)
390	zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
391	END DO
392	END DO
393	ztcond(:,nlay+1)=ztcond(:,nlay)
394	else
395	zhcond(:,:) = 0
396	ztcond(:,:) = 0
397	end if
398
399
400	c-----------------------------------------------------------------------
401	c 2. ajout des tendances physiques
402	c -----------------------------
403
404	DO ilev=1,nlay
405	DO ig=1,ngrid
406	zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
407	zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
408	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
409	! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
410	ENDDO
411	ENDDO
412
413	zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+
414	& pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
415
416	c-----------------------------------------------------------------------
417	c 3. schema de turbulence
418	c --------------------
419
420	c ** source d'energie cinetique turbulente a la surface
421	c (condition aux limites du schema de diffusion turbulente
422	c dans la couche limite
423	c ---------------------
424
425	CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar,
426	& ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
427	& pqsurf(:,igcm_h2o_ice,:),
428	& zcdv_true,zcdh_true)
429
430	zu2(:)=pu(:,1)pu(:,1)+pv(:,1)pv(:,1)
431
432	DO islope = 1,nslope
433	IF (callrichsl) THEN
434	zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
435	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
436	zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
437	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
438	ELSE
439	zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
440	zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
441	ENDIF
442	ENDDO
443	ustar(:) = 0
444	tstar(:) = 0
445	DO ig = 1,ngrid
446	zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
447	zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
448	zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
449	zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
450	IF (callrichsl) THEN
451	ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
452	& sqrt(zu2(ig)+(log(1.+0.7wstar(ig) +
453	& 2.3wstar(ig)2))*2)
454	IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN ! When Cd=Ch=0, u=t=0
455	tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
456	& *zcdh_tmp(ig)/ustar(ig)
457	ENDIF
458	ELSE
459	ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
460	& *sqrt(zu2(ig))
461	tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
462	& *zcdh_true(ig,major_slope(ig))
463	& /sqrt(zcdv_true(ig,major_slope(ig)))
464	ENDIF
465	ENDDO
466
467	c ** schema de diffusion turbulente dans la couche limite
468	c ----------------------------------------------------
469	pt(:,:)=ph(:,:)*ppopsk(:,:)
470	if (callyamada4) then
471	call yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
472	s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar
473	s ,9)
474
475	elseif (callatke) then
476	call call_atke(ptimestep,ngrid,nlay,zcdv_true_tmp,
477	s zcdh_true_tmp,pu(:,1),pv(:,1),ptsrf_tmp,
478	s pu,pv,pt,zq(:,1,igcm_h2o_vap),pplay,pplev,
479	s pzlay,pzlev,pq2,zkv(:,1:nlay),zkh(:,1:nlay))
480
481	zkv(:,nlay+1) = zkv(:,nlay)
482	zkh(:,nlay+1) = zkh(:,nlay)
483	else
484	call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
485	s ,pu,pv,ph,zcdv_true_tmp
486	s ,pq2,zkv,zkh,zq)
487
488	endif
489	if ((doubleq).and.(ngrid.eq.1)) then
490	kmixmin = 80. !80.! minimum eddy mix coeff in 1D
491	do ilev=2,nlay
492	do ig=1,ngrid
493	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
494	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
495	end do
496	end do
497	end if
498
499	c ** diagnostique pour le schema de turbulence
500	c -----------------------------------------
501
502	IF(lecrit) THEN
503	PRINT*
504	PRINT*,'Diagnostic for the vertical turbulent mixing'
505	PRINT*,'Cd for momentum and potential temperature'
506
507	PRINT*,zcdv_tmp(ngrid/2+1),zcdh_tmp(ngrid/2+1)
508	PRINT*,'Mixing coefficient for momentum and pot.temp.'
509	DO ilev=1,nlay
510	PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
511	ENDDO
512	ENDIF
513
514	c-----------------------------------------------------------------------
515	c 4. inversion pour l'implicite sur u
516	c --------------------------------
517
518	c ** l'equation est
519	c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
520	c avec
521	c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
522	c et
523	c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
524	c donc les entrees sont /zcdv/ pour la condition a la limite sol
525	c et /zkv/ = Ku
526
527	zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
528	zb(1:ngrid,1)=zcdv_tmp(1:ngrid)*zb0(1:ngrid,1)
529
530	DO ig=1,ngrid
531	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
532	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
533	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
534	ENDDO
535
536	DO ilay=nlay-1,1,-1
537	DO ig=1,ngrid
538	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
539	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
540	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
541	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
542	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
543	ENDDO
544	ENDDO
545
546	DO ig=1,ngrid
547	zu(ig,1)=zc(ig,1)
548	ENDDO
549	DO ilay=2,nlay
550	DO ig=1,ngrid
551	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
552	ENDDO
553	ENDDO
554
555	c-----------------------------------------------------------------------
556	c 5. inversion pour l'implicite sur v
557	c --------------------------------
558
559	c ** l'equation est
560	c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
561	c avec
562	c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
563	c et
564	c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
565	c donc les entrees sont /zcdv/ pour la condition a la limite sol
566	c et /zkv/ = Kv
567
568	DO ig=1,ngrid
569	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
570	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
571	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
572	ENDDO
573
574	DO ilay=nlay-1,1,-1
575	DO ig=1,ngrid
576	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
577	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
578	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
579	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
580	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
581	ENDDO
582	ENDDO
583
584	DO ig=1,ngrid
585	zv(ig,1)=zc(ig,1)
586	ENDDO
587	DO ilay=2,nlay
588	DO ig=1,ngrid
589	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
590	ENDDO
591	ENDDO
592
593	c-----------------------------------------------------------------------
594	c Using the wind modified by friction for lifting and sublimation
595	c ----------------------------------------------------------------
596
597	! This is computed above and takes into account surface-atmosphere flux
598	! enhancement by subgrid gustiness and atmospheric-stability related
599	! variations of transfer coefficients.
600	! Calculate Cd again with wind slowed by friction
601	c -------------------------------------------
602
603	CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar,
604	& ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
605	& pqsurf(:,igcm_h2o_ice,:),
606	& zcdv_true,zcdh_true)
607
608	zu2(:)=zu(:,1)zu(:,1)+zv(:,1)zv(:,1)
609
610	DO islope = 1,nslope
611	IF (callrichsl) THEN
612	zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
613	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
614	zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
615	& (log(1.+0.7wstar(:) + 2.3wstar(:)2))2)
616	ELSE
617	zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance
618	zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance
619	ENDIF
620	ENDDO
621	ustar(:) = 0
622	tstar(:) = 0
623	DO ig = 1,ngrid
624	zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
625	zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
626	zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
627	zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
628	IF (callrichsl) THEN
629	ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
630	& sqrt(zu2(ig)+(log(1.+0.7wstar(ig) +
631	& 2.3wstar(ig)2))*2)
632	IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN ! When Cd=Ch=0, u=t=0
633	tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
634	& *zcdh_tmp(ig)/ustar(ig)
635	ENDIF
636	ELSE
637	ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
638	& *sqrt(zu2(ig))
639	tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
640	& *zcdh_true(ig,major_slope(ig))
641	& /sqrt(zcdv_true(ig,major_slope(ig)))
642	ENDIF
643	ENDDO
644
645
646	c-----------------------------------------------------------------------
647	c 6. inversion pour l'implicite sur h sans oublier le couplage
648	c avec le sol (conduction)
649	c ------------------------
650
651	c ** l'equation est
652	c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
653	c avec
654	c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
655	c et
656	c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
657	c donc les entrees sont /zcdh/ pour la condition de raccord au sol
658	c et /zkh/ = Kh
659	c -------------
660
661	c Mass variation scheme:
662	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
663	zb(1:ngrid,1)=zcdh_tmp(1:ngrid)*zb0(1:ngrid,1)
664
665	c on initialise dm c
666
667	pdtc(:,:)=0.
668	zt(:,:)=0.
669	dmice(:,:)=0.
670
671	c ** calcul de (d Planck / dT) a la temperature d'interface
672	c ------------------------------------------------------
673
674	z4st=4.5.67e-8ptimestep
675	IF (tke_heat_flux .eq. 0.) THEN
676	DO ig=1,ngrid
677	indmax = major_slope(ig)
678	zdplanck(ig)=z4stpemis(ig,indmax)ptsrf(ig,indmax)*
679	& ptsrf(ig,indmax)*ptsrf(ig,indmax)
680	ENDDO
681	ELSE
682	zdplanck(:)=0.
683	ENDIF
684
685	! calcul de zc et zd pour la couche top en prenant en compte le terme
686	! de variation de masse (on fait une boucle pour que \E7a converge)
687
688	! Identification des points de grilles qui ont besoin de la correction
689
690	llnt(:)=1
691	IF (.not.turb_resolved) THEN
692	IF (callcond) THEN
693	DO ig=1,ngrid
694	DO l=1,nlay
695	if(zh(ig,l) .lt. zhcond(ig,l)) then
696	llnt(ig)=300
697	! 200 and 100 do not go beyond month 9 with normal dissipation
698	goto 5
699	endif
700	ENDDO
701	5 continue
702	ENDDO
703	ENDIF
704
705	ENDIF
706
707	DO ig=1,ngrid
708	indmax = major_slope(ig)
709	! Initialization of z1 and zd, which do not depend on dmice
710
711	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
712	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
713
714	DO ilay=nlay-1,1,-1
715	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
716	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
717	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
718	ENDDO
719
720	! Convergence loop
721
722	DO j=1,llnt(ig)
723
724	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
725	zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
726	& -betam(ig,nlay)*dmice(ig,nlay)
727	zc(ig,nlay)=zc(ig,nlay)*z1(ig)
728	! zd(ig,nlay)=zb(ig,nlay)*z1(ig)
729
730	! calcul de zc et zd pour les couches du haut vers le bas
731
732	DO ilay=nlay-1,1,-1
733	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
734	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
735	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
736	$ zb(ig,ilay+1)*zc(ig,ilay+1)-
737	$ betam(ig,ilay)dmice(ig,ilay))z1(ig)
738	! zd(ig,ilay)=zb(ig,ilay)*z1(ig)
739	ENDDO
740
741	c ** calcul de la temperature_d'interface et de sa tendance.
742	c on ecrit que la somme des flux est nulle a l'interface
743	c a t + \delta t,
744	c c'est a dire le flux radiatif a {t + \delta t}
745	c + le flux turbulent a {t + \delta t}
746	c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
747	c (notation K dt = /cpp*b/)
748	c + le flux dans le sol a t
749	c + l'evolution du flux dans le sol lorsque la temperature d'interface
750	c passe de sa valeur a t a sa valeur a {t + \delta t}.
751	c ----------------------------------------------------
752
753	z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax)
754	s + cppzb(ig,1)zc(ig,1)
755	s + zdplanck(ig)*ptsrf(ig,indmax)
756	s + pfluxsrf(ig,indmax)*ptimestep
757	z2(ig)= pcapcal(ig,indmax)+cppzb(ig,1)(1.-zd(ig,1))
758	s +zdplanck(ig)
759	ztsrf2(ig)=z1(ig)/z2(ig)
760	! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop
761	zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
762
763	c ** et a partir de la temperature au sol on remonte
764	c -----------------------------------------------
765
766	DO ilay=2,nlay
767	zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
768	ENDDO
769	DO ilay=1,nlay
770	zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
771	ENDDO
772
773	c Condensation/sublimation in the atmosphere
774	c ------------------------------------------
775	c (computation of zcondicea and dmice)
776
777	IF (.NOT. co2clouds) then
778	DO l=nlay , 1, -1
779	IF(zt(ig,l).LT.ztcond(ig,l)) THEN
780	pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
781	zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
782	& ccondpdtc(ig,l)
783	dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
784	END IF
785	ENDDO
786	ELSE
787	DO l=nlay , 1, -1
788	zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)*
789	c & (pplev(ig,l) - pplev(ig,l+1))/g
790	dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep
791	pdtc(ig,l)=0.
792	ENDDO
793	ENDIF
794
795	ENDDO!of Do j=1,XXX
796	pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep
797	ENDDO !of Do ig=1,ngrid
798
799	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
800	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zhs(ig,1)-ztsrf2(ig))
801	ENDDO
802
803	c Now implicit sheme on each sub-grid subslope:
804	IF (nslope.ne.1) then
805	DO islope=1,nslope
806	DO ig=1,ngrid
807	IF(islope.ne.major_slope(ig)) then
808	IF (tke_heat_flux .eq. 0.) THEN
809	zdplanck(ig)=z4stpemis(ig,islope)ptsrf(ig,islope)**3
810	ELSE
811	zdplanck(ig) = 0.
812	ENDIF
813	zb(ig,1)=zcdh(ig,islope)*zb0(ig,1)
814	z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope)
815	s + cppzb(ig,1)zc(ig,1)
816	s + zdplanck(ig)*ptsrf(ig,islope)
817	s + pfluxsrf(ig,islope)*ptimestep
818	z2(ig)= pcapcal(ig,islope)+cppzb(ig,1)(1.-zd(ig,1))
819	s +zdplanck(ig)
820	ztsrf2(ig)=z1(ig)/z2(ig)
821	pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep
822	ENDIF ! islope != indmax
823	ENDDO ! ig
824	ENDDO !islope
825	ENDIF !nslope.ne.1
826
827	c-----------------------------------------------------------------------
828	c TRACERS
829	c -------
830	c Calcul du flux vertical au bas de la premiere couche (dust) :
831	c -----------------------------------------------------------
832	do ig=1,ngrid
833	rho(ig) = zb0(ig,1) /ptimestep
834	c zb(ig,1) = 0.
835	end do
836	c Dust lifting:
837	if (lifting) then
838	#ifndef MESOSCALE
839	if (doubleq.AND.submicron) then
840	do ig=1,ngrid
841	c if(qsurf(ig,igcm_co2).lt.1) then
842	pdqsdif_tmp(ig,igcm_dust_mass) =
843	& -alpha_lift(igcm_dust_mass)
844	pdqsdif_tmp(ig,igcm_dust_number) =
845	& -alpha_lift(igcm_dust_number)
846	pdqsdif_tmp(ig,igcm_dust_submicron) =
847	& -alpha_lift(igcm_dust_submicron)
848	c end if
849	end do
850	else if (doubleq) then
851	if (dustinjection.eq.0) then !injection scheme 0 (old)
852	!or 2 (injection in CL)
853	do ig=1,ngrid
854	if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
855	pdqsdif_tmp(ig,igcm_dust_mass) =
856	& -alpha_lift(igcm_dust_mass)
857	pdqsdif_tmp(ig,igcm_dust_number) =
858	& -alpha_lift(igcm_dust_number)
859	end if
860	end do
861	elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
862	do ig=1,ngrid
863	if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
864	IF((ti_injection_sol.LE.local_time(ig)).and.
865	& (local_time(ig).LE.tf_injection_sol)) THEN
866	if (rdstorm) then !Rocket dust storm scheme
867	pdqsdif_tmp(ig,igcm_stormdust_mass) =
868	& -alpha_lift(igcm_stormdust_mass)
869	& *dustliftday(ig)
870	pdqsdif_tmp(ig,igcm_stormdust_number) =
871	& -alpha_lift(igcm_stormdust_number)
872	& *dustliftday(ig)
873	pdqsdif_tmp(ig,igcm_dust_mass)= 0.
874	pdqsdif_tmp(ig,igcm_dust_number)= 0.
875	else
876	pdqsdif_tmp(ig,igcm_dust_mass)=
877	& -dustliftday(ig)*
878	& alpha_lift(igcm_dust_mass)
879	pdqsdif_tmp(ig,igcm_dust_number)=
880	& -dustliftday(ig)*
881	& alpha_lift(igcm_dust_number)
882	endif
883	if (submicron) then
884	pdqsdif_tmp(ig,igcm_dust_submicron) = 0.
885	endif ! if (submicron)
886	ELSE ! outside dust injection time frame
887	pdqsdif_tmp(ig,igcm_dust_mass)= 0.
888	pdqsdif_tmp(ig,igcm_dust_number)= 0.
889	if (rdstorm) then
890	pdqsdif_tmp(ig,igcm_stormdust_mass)= 0.
891	pdqsdif_tmp(ig,igcm_stormdust_number)= 0.
892	end if
893	ENDIF
894
895	end if ! of if(qsurf(ig,igcm_co2).lt.1)
896	end do
897	endif ! end if dustinjection
898	else if (submicron) then
899	do ig=1,ngrid
900	pdqsdif_tmp(ig,igcm_dust_submicron) =
901	& -alpha_lift(igcm_dust_submicron)
902	end do
903	else
904	#endif
905	call dustlift(ngrid,nlay,nq,rho,zcdh_true_tmp,zcdh_tmp,
906	& pqsurf_tmp(:,igcm_co2),pdqsdif_tmp)
907	#ifndef MESOSCALE
908	endif !doubleq.AND.submicron
909	#endif
910	else
911	pdqsdif_tmp(1:ngrid,1:nq) = 0.
912	end if
913
914	c OU calcul de la valeur de q a la surface (water) :
915	c ----------------------------------------
916
917	c Inversion pour l'implicite sur q
918	c Cas des traceurs qui ne sont pas h2o_vap
919	c h2o_vap est traite plus loin avec un sous pas de temps
920	c hdo_vap est traite ensuite car dependant de h2o_vap
921	c --------------------------------
922
923	do iq=1,nq !for all tracers except water vapor
924	if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.
925	& (.not. iq.eq.igcm_hdo_vap)) then
926
927
928	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
929	zb(1:ngrid,1)=0
930
931	DO ig=1,ngrid
932	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
933	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
934	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
935	ENDDO
936
937	DO ilay=nlay-1,2,-1
938	DO ig=1,ngrid
939	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
940	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
941	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
942	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
943	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
944	ENDDO
945	ENDDO
946
947	if ((iq.eq.igcm_h2o_ice)
948	$ .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then
949
950	DO ig=1,ngrid
951	z1(ig)=1./(za(ig,1)+zb(ig,1)+
952	$ zb(ig,2)*(1.-zd(ig,2)))
953	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
954	$ zb(ig,2)zc(ig,2)) z1(ig) !special case h2o_ice
955	ENDDO
956	else ! every other tracer
957	DO ig=1,ngrid
958	z1(ig)=1./(za(ig,1)+zb(ig,1)+
959	$ zb(ig,2)*(1.-zd(ig,2)))
960	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
961	$ zb(ig,2)*zc(ig,2) +
962	$ (-pdqsdif_tmp(ig,iq)) ptimestep) z1(ig) !tracer flux from surface
963	ENDDO
964	endif !((iq.eq.igcm_h2o_ice)
965	c Starting upward calculations for simple mixing of tracer (dust)
966	DO ig=1,ngrid
967	zq(ig,1,iq)=zc(ig,1)
968	DO ilay=2,nlay
969	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
970	ENDDO
971	ENDDO
972	DO islope = 1,nslope
973	DO ig = 1,ngrid
974	pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq)
975	& * cos(pi*def_slope_mean(islope)/180.)
976	ENDDO
977	ENDDO
978
979	endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then
980	enddo ! of do iq=1,nq
981
982	c --------- h2o_vap --------------------------------
983
984
985	c Traitement de la vapeur d'eau h2o_vap
986	c Utilisation d'un sous pas de temps afin
987	c de decrire le flux de chaleur latente
988
989	do iq=1,nq
990	if ((water).and.(iq.eq.igcm_h2o_vap)) then
991
992	DO islope = 1,nslope
993	DO ig=1,ngrid
994
995
996
997	zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
998	& cos(pi*def_slope_mean(islope)/180.)
999	watercap_tmp(ig) = watercap(ig,islope)/
1000	& cos(pi*def_slope_mean(islope)/180.)
1001	ENDDO ! ig=1,ngrid
1002
1003	c make_tsub : sous pas de temps adaptatif
1004	c la subroutine est a la fin du fichier
1005	call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
1006	& ptimestep,dtmax,watercaptag,
1007	& nsubtimestep)
1008	c Calculation for turbulent exchange with the surface (for ice)
1009	c initialization of ztsrf, which is surface temperature in
1010	c the subtimestep.
1011	saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)
1012	DO ig=1,ngrid
1013	! nsubtimestep(ig)=1 !for debug
1014	subtimestep = ptimestep/nsubtimestep(ig)
1015	call write_output('subtimestep',
1016	& 'vdifc substimestep length','s',subtimestep)
1017	ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep
1018	zq_tmp_vap(ig,:,:) =zq(ig,:,:)
1019	c Debut du sous pas de temps
1020	DO tsub=1,nsubtimestep(ig)
1021	if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
1022	c C'est parti !
1023	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1024	& /float(nsubtimestep(ig))
1025	if(old_wsublimation_scheme) then
1026	zb(1:ngrid,1)=zcdv(1:ngrid,islope)*zb0(1:ngrid,1)
1027	& /float(nsubtimestep(ig))
1028	else
1029	zb(1:ngrid,1)=zcdh(1:ngrid,islope)*zb0(1:ngrid,1)
1030	& /float(nsubtimestep(ig))
1031	endif
1032	zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
1033
1034	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1035	zc(ig,nlay)=za(ig,nlay)zq_tmp_vap(ig,nlay,iq)z1(ig)
1036	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1037	DO ilay=nlay-1,2,-1
1038	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1039	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1040	zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
1041	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
1042	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1043	ENDDO
1044	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1045	$ zb(ig,2)*(1.-zd(ig,2)))
1046	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
1047	$ zb(ig,2)zc(ig,2)) z1(ig)
1048
1049	call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig))
1050	old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap)
1051	zd(ig,1)=zb(ig,1)*z1(ig)
1052	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
1053	if(old_wsublimation_scheme) then
1054	zdqsdif_surf(ig)=rho(ig)dryness(ig)zcdv(ig,islope)
1055	& *(zq1temp(ig)-qsat(ig))
1056	else
1057	zdqsdif_surf(ig)=rho(ig)dryness(ig)zcdh(ig,islope)
1058	& *(zq1temp(ig)-qsat(ig))
1059	endif
1060
1061	zdqsdif_tot(ig) = zdqsdif_surf(ig)
1062	!!! Subsurface exchange
1063	! Check for subsurface exchanges
1064	if(.not.watercaptag(ig)) then
1065	if (((-(zdqsdif_surf(ig))*
1066	& subtimestep).gt.zqsurf(ig))
1067	& .and.(pqsurf(ig,igcm_co2,islope).eq.0.)) then
1068	exchange = .true.
1069	else
1070	exchange = .false.
1071	endif
1072	else
1073	exchange = .false.
1074	endif
1075	zdqsdif_tot(ig) = zdqsdif_surf(ig)
1076
1077
1078	if (adsorption_soil) then
1079	call soilwater(1,nlay,nq,nsoil, nqsoil,
1080	& ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
1081	& exchange,qsat(ig),zq_tmp_vap(ig,:,:),
1082	& za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:),
1083	& zdqsdif_surf(ig), zqsurf(ig),
1084	& qsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
1085	& writeoutput,zdqsdif_regolith(ig,islope),
1086	& zq1temp_regolith(ig))
1087
1088
1089
1090	if(.not.watercaptag(ig)) then
1091	if (exchange) then
1092	zq1temp(ig) = zq1temp_regolith(ig)
1093	zdqsdif_tot(ig)=
1094	& -zqsurf(ig)/subtimestep
1095	else
1096	zdqsdif_tot(ig) = zdqsdif_surf(ig) +
1097	& zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
1098	endif ! of "if exchange = true"
1099	endif ! of "if not.watercaptag"
1100	endif ! adsorption
1101
1102	if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then
1103	if ((-zdqsdif_tot(ig)*subtimestep)
1104	& .gt.(zqsurf(ig))) then
1105
1106
1107	!EV subsurface ice
1108	IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
1109	& then
1110	zdqsdif_tot(ig)=
1111	& -zqsurf(ig)/subtimestep
1112	zqsurf(ig)=0
1113
1114	DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
1115	IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
1116	& (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
1117	lice = ik+1
1118	EXIT
1119	ENDIF
1120	ENDDO !of subsurface loop
1121	IF (lice .gt. 1) then !calculate the distance from the layers
1122	dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
1123	& -mlayer(lice-1))
1124	dist_down(ig,islope)=(mlayer(lice)
1125	& -h2o_ice_depth(ig,islope))
1126	dist_sum(ig,islope)=dist_up(ig,islope)
1127	& +dist_down(ig,islope)
1128	Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
1129	& *ptsoil(ig,lice-1,islope)
1130	& /dist_sum(ig,islope))+
1131	& (dist_down(ig,islope)*ptsoil(ig,lice,islope)
1132	& /dist_sum(ig,islope))
1133	ELSE
1134	Tice(ig,islope)=ptsoil(ig,1,islope)
1135	ENDIF
1136	call watersat(1,Tice(ig,1),pplev(ig,1)
1137	& ,qsat2(ig,1))
1138	qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
1139	& *qsat2(ig,1)
1140	resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig,islope)
1141	& /d_coef(ig,islope)))
1142	!write(,)'R=',resist(ig,islope)
1143	!write(,)'zice=',h2o_ice_depth(ig,islope)
1144	!write(,)'D=',d_coef(ig,islope)
1145	!write(,)'zcdh=',zcdh(ig)
1146	zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice
1147	!vdifc algorithem !!!!!needs to change reseist io to the mean
1148	!beacuse of the slopes!!!!!
1149	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1150	zc(ig,nlay)=za(ig,nlay)zq_tmp_vap(ig,nlay,iq)z1(ig)
1151	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1152	DO ilay=nlay-1,2,-1
1153	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1154	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1155	zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
1156	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
1157	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1158	ENDDO
1159	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1160	$ zb(ig,2)*(1.-zd(ig,2)))
1161	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
1162	$ zb(ig,2)zc(ig,2)) z1(ig)
1163	zd(ig,1)=zb(ig,1)*z1(ig)
1164	zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
1165	zdqsdif_ssi(ig,1)=rho(ig)dryness(ig)zcdv(ig,islope)
1166	& *(zq1temp(ig)-qeq(ig,islope))
1167	call write_output('zdq_ssi',
1168	& '','',zdqsdif_ssi(ig,1))
1169	call write_output('qeq',
1170	& '','',qeq(ig,1))
1171	call write_output('q1',
1172	& '','',zq1temp(ig))
1173	!write(,)'qeq=',qeq(ig,1)
1174	!write(,)'q1=',zq1temp(ig)
1175	!write(,)'zdqsdif_ssi=',zdqsdif_ssi(ig,1)
1176	!I should some all the interactions with the SSI accoridng
1177	!to the statistics and the evaluate
1178
1179	!write(,) "zdq_ssit", zdqsdif_ssi(ig)subtimestep
1180	if (zdqsdif_ssi(ig,1)<0) then
1181	!zdqsdif_surf(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
1182	zdqsdif_tot(ig)=zdqsdif_tot(ig)+zdqsdif_ssi(ig,1)
1183	call write_output('zdq_zdqssi',
1184	& '','',zdqsdif_tot(ig)+zdqsdif_ssi(ig,1))
1185	endif
1186	ELSE
1187	c pdqsdif > 0 : ice condensing
1188	c pdqsdif < 0 : ice subliming
1189	c write(,) "subliming more than available frost: qsurf!"
1190	zdqsdif_tot(ig)=
1191	& -zqsurf(ig)/subtimestep
1192	c write(,)'flux vers le sol=',pdqsdif(ig,nq)
1193	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
1194	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+
1195	$ zb(ig,2)*zc(ig,2) +
1196	$ (-zdqsdif_tot(ig)) subtimestep) z1(ig)
1197	zq1temp(ig)=zc(ig,1)
1198	ENDIF !if h2o_ice_depth>0 and lag_layer
1199	endif !if .not.watercaptag(ig)
1200	endif ! if sublim more than surface
1201
1202	! subsurface ice <--> frost interaction !EV
1203	IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
1204	& .and. zqsurf(ig) .gt. 0) then
1205	DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
1206	IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
1207	& (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
1208	lice = ik+1
1209	EXIT
1210	ENDIF
1211	ENDDO !of subsurface loop
1212	IF (lice .gt. 1) then !calculate the distance from the layers
1213	dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
1214	& -mlayer(lice-1))
1215	dist_down(ig,islope)=(mlayer(lice)
1216	& -h2o_ice_depth(ig,islope))
1217	dist_sum(ig,islope)=dist_up(ig,islope)
1218	& +dist_down(ig,islope)
1219	Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
1220	& *ptsoil(ig,lice-1,islope)
1221	& /dist_sum(ig,islope))+
1222	& (dist_down(ig,islope)*ptsoil(ig,lice,islope)
1223	& /dist_sum(ig,islope))
1224	ELSE
1225	Tice(ig,islope)=ptsoil(ig,1,islope)
1226	ENDIF
1227
1228	call watersat(1,Tice(ig,1),pplev(ig,1)
1229	& ,qsat2(ig,1))
1230	qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
1231	& *qsat2(ig,1)
1232	! write(,)'icedep=',h2o_ice_depth(ig,1)
1233	! write(,)'qeq=',qeq(ig,1)
1234	! write(,)'d=',d_coef(ig,1)
1235	! write(,)'qsat=',qsat(ig)
1236	! write(,)'dry=',dryness(ig)
1237	! write(,)'rho=',rho(ig)
1238	zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
1239	& /h2o_ice_depth(ig,1))
1240	& rho(ig)dryness(ig)*(qsat(ig)-qeq(ig,1))
1241	!needs to change to the mean of eq
1242	zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
1243	!!!!! zdsqdif_ssi_frosst need to be changed to an
1244	!average
1245	ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
1246	& .and. watercaptag(ig)) then
1247	DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
1248	IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
1249	& (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
1250	lice = ik+1
1251	EXIT
1252	ENDIF
1253	ENDDO !of subsurface loop
1254	IF (lice .gt. 1) then !calculate the distance from the layers
1255	dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
1256	& -mlayer(lice-1))
1257	dist_down(ig,islope)=(mlayer(lice)
1258	& -h2o_ice_depth(ig,islope))
1259	dist_sum(ig,islope)=dist_up(ig,islope)
1260	& +dist_down(ig,islope)
1261	Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
1262	& *ptsoil(ig,lice-1,islope)
1263	& /dist_sum(ig,islope))+
1264	& (dist_down(ig,islope)*ptsoil(ig,lice,islope)
1265	& /dist_sum(ig,islope))
1266	ELSE
1267	Tice(ig,islope)=ptsoil(ig,1,islope)
1268	ENDIF
1269
1270	call watersat(1,Tice(ig,1),pplev(ig,1)
1271	& ,qsat2(ig,1))
1272	qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
1273	& *qsat2(ig,1)
1274
1275	zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
1276	& /h2o_ice_depth(ig,1))
1277	& rho(ig)dryness(ig)*(qsat(ig)-qeq(ig,1))
1278	zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
1279	!needs to change to the mean of eq
1280	ENDIF
1281	! call write_output('subtimestep',
1282	! & 'vdifc substimestep length','s',subtimestep)
1283	! ENDDO !subsurface ice subslope
1284
1285
1286
1287	c Starting upward calculations for water :
1288	c Actualisation de h2o_vap dans le premier niveau
1289	zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
1290	c Take into account the H2O latent heat impact on the surface temperature
1291	if (latentheat_surfwater) then
1292	lh=(2834.3-0.28*(ztsrf(ig)-To)-
1293	& 0.004(ztsrf(ig)-To)(ztsrf(ig)-To))*1.e+3
1294	zdtsrf(ig,islope)= zdqsdif_tot(ig)*lh
1295	& /pcapcal(ig,islope)
1296	endif ! (latentheat_surfwater) then
1297
1298	DO ilay=2,nlay
1299	zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)
1300	& *zq_tmp_vap(ig,ilay-1,iq)
1301	ENDDO
1302	c Subtimestep water budget :
1303	ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope)
1304	& + zdtsrf(ig,islope))*subtimestep
1305	zqsurf(ig)= zqsurf(ig)+(
1306	& zdqsdif_tot(ig))*subtimestep
1307	if (zqsurf(ig)<0 .and.
1308	& (.not.watercaptag(ig))) then
1309	zqsurf(ig)=0
1310	endif
1311
1312	c Monitoring instantaneous latent heat flux in W.m-2 :
1313	zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
1314	& (zdtsrf(ig,islope)*pcapcal(ig,islope))
1315	& *subtimestep
1316
1317	c We ensure that surface temperature can't rise above the solid domain if there
1318	c is still ice on the surface (oldschool)
1319	if(zqsurf(ig)
1320	& +zdqsdif_tot(ig)*subtimestep
1321	& .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
1322	zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
1323	& (To-ztsrf(ig))/subtimestep) ! ice melt case
1324	endif
1325
1326	c Fin du sous pas de temps
1327	ENDDO ! tsub=1,nsubtimestep
1328
1329	c Integration of subtimestep temp and water budget :
1330	c (btw could also compute the post timestep temp and ice
1331	c by simply adding the subtimestep trend instead of this)
1332	surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep
1333	pdtsrf(ig,islope)= (ztsrf(ig) -
1334	& ptsrf(ig,islope))/ptimestep
1335	pdqsdif(ig,igcm_h2o_ice,islope)=
1336	& (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/
1337	& cos(pi*def_slope_mean(islope)/180.))
1338	& /ptimestep
1339	c if subliming more than qsurf(ice) and on watercaptag, water
1340	c sublimates from watercap reservoir (dwatercap_dif is <0)
1341	if(watercaptag(ig)) then
1342	if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep)
1343	& .gt.(pqsurf(ig,igcm_h2o_ice,islope)
1344	& /cos(pi*def_slope_mean(islope)/180.))) then
1345	dwatercap_dif(ig,islope)=
1346	& pdqsdif(ig,igcm_h2o_ice,islope)+
1347	& (pqsurf(ig,igcm_h2o_ice,islope) /
1348	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1349	pdqsdif(ig,igcm_h2o_ice,islope)=
1350	& - (pqsurf(ig,igcm_h2o_ice,islope)/
1351	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1352	endif! ((-pdqsdif(ig)*ptimestep)
1353	endif !(watercaptag(ig)) then
1354	zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:)
1355	ENDDO ! of DO ig=1,ngrid
1356	ENDDO ! islope
1357	c Some grid box averages: interface with the atmosphere
1358	DO ig = 1,ngrid
1359	DO ilay = 1,nlay
1360	zq(ig,ilay,iq) = 0.
1361	DO islope = 1,nslope
1362	zq(ig,ilay,iq) = zq(ig,ilay,iq) +
1363	$ zq_slope_vap(ig,ilay,iq,islope) *
1364	$ subslope_dist(ig,islope)
1365	ENDDO
1366	ENDDO
1367	ENDDO
1368	! Recompute values in kg/m^2 slopped
1369	DO ig = 1,ngrid
1370	DO islope = 1,nslope
1371	pdqsdif(ig,igcm_h2o_ice,islope) =
1372	& pdqsdif(ig,igcm_h2o_ice,islope)
1373	& * cos(pi*def_slope_mean(islope)/180.)
1374
1375	dwatercap_dif(ig,islope) =
1376	& dwatercap_dif(ig,islope)
1377	& * cos(pi*def_slope_mean(islope)/180.)
1378	ENDDO
1379	ENDDO
1380
1381	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
1382
1383	c --------- end of h2o_vap ----------------------------
1384
1385	c --------- hdo_vap -----------------------------------
1386
1387	c hdo_ice has already been with along h2o_ice
1388	c amongst "normal" tracers (ie not h2o_vap)
1389
1390	if (hdo.and.(iq.eq.igcm_hdo_vap)) then
1391	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1392	zb(1:ngrid,1)=0
1393
1394	DO ig=1,ngrid
1395	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1396	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
1397	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1398	ENDDO
1399
1400	DO ilay=nlay-1,2,-1
1401	DO ig=1,ngrid
1402	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1403	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1404	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
1405	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
1406	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1407	ENDDO
1408	ENDDO
1409	hdoflux_meshavg(:) = 0.
1410	DO islope = 1,nslope
1411
1412	pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope)
1413	& /cos(pi*def_slope_mean(islope)/180.)
1414
1415	call watersat(ngrid,pdtsrf(:,islope)*ptimestep +
1416	& ptsrf(:,islope),pplev(:,1),qsat_tmp)
1417
1418	CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
1419	& zt,pplay,zq,pqsurf(:,:,islope),
1420	& saved_h2o_vap,qsat_tmp,
1421	& pdqsdif_tmphdo,
1422	& dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.),
1423	& hdoflux(:,islope))
1424
1425	pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) *
1426	& cos(pi*def_slope_mean(islope)/180.)
1427	DO ig = 1,ngrid
1428	hdoflux_meshavg(ig) = hdoflux_meshavg(ig) +
1429	& hdoflux(ig,islope)*subslope_dist(ig,islope)
1430
1431	ENDDO !ig = 1,ngrid
1432	ENDDO !islope = 1,nslope
1433
1434	DO ig=1,ngrid
1435	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1436	$ zb(ig,2)*(1.-zd(ig,2)))
1437	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
1438	$ zb(ig,2)*zc(ig,2) +
1439	$ (-hdoflux_meshavg(ig)) ptimestep) z1(ig) !tracer flux from surface
1440	ENDDO
1441
1442	DO ig=1,ngrid
1443	zq(ig,1,iq)=zc(ig,1)
1444	DO ilay=2,nlay
1445	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
1446	ENDDO
1447	ENDDO
1448	endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
1449
1450	c --------- end of hdo ----------------------------
1451
1452	enddo ! of do iq=1,nq
1453
1454	c --------- end of tracers ----------------------------
1455
1456	call write_output("surf_h2o_lh",
1457	& "Ground ice latent heat flux",
1458	& "W.m-2",surf_h2o_lh(:,iflat))
1459	! call write_output('zdqsdif_ssi_frost',
1460	! & 'Flux between frost and subsurface','kg.m-2.s-1',
1461	! & zdqsdif_ssi_frost(:,1))
1462
1463	! call write_output('zdq_subtimestep',
1464	! & 'Actual flux zdqsdif_surf*subtimestep',
1465	! & 'kg.m-2',zdqsdif_tot(:)*subtimestep)
1466	! call write_output('zdq_end',
1467	! & 'Flux after all contributions',
1468	! & 'kg.m-2.s-1',zdqsdif_tot(:))
1469	C Diagnostic output for HDO
1470	! if (hdo) then
1471	! CALL write_output('hdoflux',
1472	! & 'hdoflux',
1473	! & ' ',hdoflux_meshavg(:))
1474	! CALL write_output('h2oflux',
1475	! & 'h2oflux',
1476	! & ' ',h2oflux(:))
1477	! endif
1478
1479	c-----------------------------------------------------------------------
1480	c 8. calcul final des tendances de la diffusion verticale
1481	c ----------------------------------------------------
1482
1483	DO ilev = 1, nlay
1484	DO ig=1,ngrid
1485	pdudif(ig,ilev)=( zu(ig,ilev)-
1486	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
1487	pdvdif(ig,ilev)=( zv(ig,ilev)-
1488	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
1489	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
1490	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
1491	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
1492	ENDDO
1493	ENDDO
1494
1495	pdqdif(1:ngrid,1:nlay,1:nq)=(zq(1:ngrid,1:nlay,1:nq)-
1496	& (pq(1:ngrid,1:nlay,1:nq)
1497	& +pdqfi(1:ngrid,1:nlay,1:nq)
1498	& *ptimestep))/ptimestep
1499
1500	c ** diagnostique final
1501	c ------------------
1502
1503	IF(lecrit) THEN
1504	PRINT*,'In vdif'
1505	PRINT*,'Ts (t) and Ts (t+st)'
1506	WRITE(*,'(a10,3a15)')
1507	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
1508	PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1)
1509	DO ilev=1,nlay
1510	WRITE(*,'(4f15.7)')
1511	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
1512	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
1513
1514	ENDDO
1515	ENDIF
1516
1517	END SUBROUTINE vdifc
1518
1519	c====================================
1520
1521	SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
1522	$ dtmax,watercaptag,ntsub)
1523
1524	c Pas de temps adaptatif en estimant le taux de sublimation
1525	c et en adaptant avec un critere "dtmax" du chauffage a accomoder
1526	c dtmax est regle empiriquement (pour l'instant) a 0.5 K
1527
1528	integer,intent(in) :: naersize
1529	real,intent(in) :: dtsurf(naersize)
1530	real,intent(in) :: qsurf(naersize)
1531	logical,intent(in) :: watercaptag(naersize)
1532	real,intent(in) :: ptimestep
1533	real,intent(in) :: dtmax
1534	real :: ztsub(naersize)
1535	integer :: i
1536	integer,intent(out) :: ntsub(naersize)
1537
1538	do i=1,naersize
1539	if ((qsurf(i).eq.0).and.
1540	& (.not.watercaptag(i))) then
1541	ntsub(i) = 1
1542	else
1543	ztsub(i) = ptimestep * dtsurf(i) / dtmax
1544	ntsub(i) = ceiling(abs(ztsub(i)))
1545	endif ! (qsurf(i).eq.0) then
1546	c
1547	c write(78,), dtsurfptimestep, dtsurf, ntsub
1548	enddo! 1=1,ngrid
1549
1550
1551
1552	END SUBROUTINE make_tsub
1553	END MODULE vdifc_mod

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

Download in other formats:

Original Format