Context Navigation

vdifc_mod.F @ 3384

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

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

Download in other formats:

Original Format