Context Navigation

vdifc_mod.F @ 3468

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

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

Download in other formats:

Original Format