Context Navigation

vdifc_mod.F @ 3307

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

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

Download in other formats:

Original Format