Context Navigation

vdifc_mod.F @ 3807

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

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

Download in other formats:

Original Format