Context Navigation

vdifc_mod.F @ 3263

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

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

Download in other formats:

Original Format