Context Navigation

vdifc_mod.F @ 3400

Last change on this file since 3400 was 3400, checked in by llange, 4 months ago

Mars PCM
Back to the commit
Correction of the computation of the water vapor flux from the subsurface ice when thin frost is at the surface: When the frost was too thin (thickness < tol_frost), the flux from the subsurface ice to the surface/atmosphere was not computed. It is now corrected, assuming that the frost is too thin to prevent exchange between subsurface ice and the atmosphere
cleaning testphys1d (unused comments)
LL, who should be writing his PhD dissertation

File size: 61.5 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))))
1117	& .or.
1118	& (((-zdqsdif_tot(ig)*subtimestep)
1119	& .lt.(zqsurf(ig)))
1120	& .and. ((zqsurf(ig).lt.tol_frost))
1121	& .and.lag_layer)) then
1122
1123
1124	! zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1125	if((h2o_ice_depth(ig,islope).gt.0)
1126	& .and. (zqsurf(ig) .lt. tol_frost)
1127	& .and.lag_layer) then
1128	! 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
1129	! zqsurf(ig) = 0
1130	if(old_wsublimation_scheme) then
1131	resist(ig,islope)=h2o_ice_depth(ig,islope)
1132	& *zcdv(ig,islope)/d_coef(ig,islope)
1133	else
1134	resist(ig,islope)=h2o_ice_depth(ig,islope)
1135	& *zcdh(ig,islope)/d_coef(ig,islope)
1136	endif
1137
1138	zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice
1139	! Now we add the flux from the subsurface ice : rho Cd,h U(1/1+R) (q1-qsat_ssi(Tice))
1140
1141	call compute_Tice(nsoil, ptsoil(ig,:,islope),
1142	& ztsrf(ig),
1143	& h2o_ice_depth(ig,islope),
1144	& Tice(ig,islope)) ! compute ice temperature
1145
1146	call watersat(1,Tice(ig,islope),pplev(ig,1),
1147	& qsat_ssi(ig,islope))
1148
1149	qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1150	& *qsat_ssi(ig,islope)
1151
1152	! And now we solve correctly the system
1153	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1154	& zb(ig,2)*(1.-zd(ig,2)))
1155	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
1156	& zb(ig,2)*zc(ig,2) +
1157	& (-zdqsdif_tot(ig)) *subtimestep)
1158	& * z1(ig)
1159	zd(ig,1)=zb(ig,1)*z1(ig)
1160	zq1temp(ig)=zc(ig,1)+ zd(ig,1)
1161	& *qsat_ssi(ig,islope)
1162
1163	! Flux from the subsurface
1164	if(old_wsublimation_scheme) then
1165	zdqsdif_ssi_atm(ig,islope)=rho(ig)
1166	& dryness(ig)zcdv(ig,islope)
1167	& *1/(1+resist(ig,islope))
1168	& *(zq1temp(ig)-qsat_ssi(ig,islope))
1169	else
1170	zdqsdif_ssi_atm(ig,islope)=rho(ig)*
1171	& dryness(ig) zcdh(ig,islope)
1172	& *1/(1+resist(ig,islope))
1173	& *(zq1temp(ig)-qsat_ssi(ig,islope))
1174	endif
1175
1176	else
1177	! No atm <-> subsurface exchange, we do it the usual way
1178	zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1179	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
1180	zc(ig,1)=(za(ig,1)*
1181	& zq_tmp_vap(ig,1,igcm_h2o_vap)+
1182	& zb(ig,2)*zc(ig,2) +
1183	& (-zdqsdif_tot(ig)) subtimestep) z1(ig)
1184	zq1temp(ig)=zc(ig,1)
1185	endif ! Of subsurface <-> atmosphere exchange
1186	endif ! sublimating more than available frost & surface - frost exchange
1187	endif !if .not.watercaptag(ig)
1188
1189	! --------------------------------------------------------------------------------------------------------------------------------
1190	! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet),
1191	! we do not include their effect on the mass of surface frost.
1192	! --------------------------------------------------------------------------------------------------------------------------------
1193
1194	if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer
1195	& .and.(.not.adsorption_soil)) then
1196	! First case: still frost at the surface but no watercaptag
1197	if(((watercaptag(ig))).or.
1198	& (((-zdqsdif_tot(ig)*subtimestep)
1199	& .lt.(zqsurf(ig)))
1200	& .and. (zqsurf(ig).gt.tol_frost))) then
1201	! Still frost at the surface: we consider the possibility to have subsurface <-> frost exchange
1202	! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice))
1203	call compute_Tice(nsoil, ptsoil(ig,:,islope),
1204	& ztsrf(ig),
1205	& h2o_ice_depth(ig,islope),
1206	& Tice(ig,islope)) ! compute ice temperature
1207
1208	call watersat(1,Tice(ig,islope),pplev(ig,1),
1209	& qsat_ssi(ig,islope))
1210
1211	qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1212	& *qsat_ssi(ig,islope)
1213
1214	zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope)
1215	& /h2o_ice_depth(ig,islope)
1216	& rho(ig)dryness(ig)
1217	& *(qsat(ig)-qsat_ssi(ig,islope))
1218
1219	! Line to comment for now since we don't have a mass of subsurface frost
1220	! in our computation (otherwise, we would not conserve the H2O mass in
1221	! the system)
1222	if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope))
1223	& *subtimestep.ge.(zqsurf(ig))) then
1224	zdqsdif_ssi_frost(ig,islope) =
1225	& zqsurf(ig)/subtimestep
1226	& + zdqsdif_tot(ig)
1227	endif
1228
1229	zdqsdif_tot(ig) = zdqsdif_tot(ig) -
1230	& zdqsdif_ssi_frost(ig,islope)
1231	endif ! watercaptag or frost at the surface
1232	endif ! interaction frost <-> subsurface ice
1233
1234
1235	c Starting upward calculations for water :
1236	zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
1237	c Take into account the H2O latent heat impact on the surface temperature
1238	if (latentheat_surfwater) then
1239	lh=(2834.3-0.28*(ztsrf(ig)-To)-
1240	& 0.004(ztsrf(ig)-To)(ztsrf(ig)-To))*1.e+3
1241	zdtsrf(ig,islope)= zdqsdif_tot(ig)*lh
1242	& /pcapcal(ig,islope)
1243	endif ! (latentheat_surfwater) then
1244
1245	DO ilay=2,nlay
1246	zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)
1247	& *zq_tmp_vap(ig,ilay-1,iq)
1248	ENDDO
1249	c Subtimestep water budget :
1250	ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope)
1251	& + zdtsrf(ig,islope))*subtimestep
1252	zqsurf(ig)= zqsurf(ig)+(
1253	& zdqsdif_tot(ig))*subtimestep
1254
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	use surfdat_h, only: frost_albedo_threshold
1492
1493	integer,intent(in) :: naersize
1494	real,intent(in) :: dtsurf(naersize)
1495	real,intent(in) :: qsurf(naersize)
1496	logical,intent(in) :: watercaptag(naersize)
1497	real,intent(in) :: ptimestep
1498	real,intent(in) :: dtmax
1499	real :: ztsub(naersize)
1500	integer :: i
1501	integer,intent(out) :: ntsub(naersize)
1502
1503	do i=1,naersize
1504	c If not on permanent ice or thin frost layer :
1505	c do not use a subtimestep (ntsub=1)
1506	if ((qsurf(i).lt.frost_albedo_threshold).and.
1507	& (.not.watercaptag(i))) then
1508	ntsub(i) = 1
1509	else
1510	ztsub(i) = ptimestep * dtsurf(i) / dtmax
1511	ntsub(i) = ceiling(abs(ztsub(i)))
1512	endif ! (qsurf(i).eq.0) then
1513	c
1514	c write(78,), dtsurfptimestep, dtsurf, ntsub
1515	enddo! 1=1,ngrid
1516
1517
1518
1519	END SUBROUTINE make_tsub
1520
1521
1522	c====================================
1523
1524	SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice)
1525
1526	c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
1527	use comsoil_h, only: layer, mlayer
1528
1529	implicit none
1530	integer,intent(in) :: nsoil ! Number of soil layers
1531	real,intent(in) :: ptsoil(nsoil) ! Soil temperature (K)
1532	real,intent(in) :: ptsurf ! Soil temperature (K)
1533	real,intent(in) :: ice_depth ! Ice depth (m)
1534	real,intent(out) :: Tice ! Ice temperatures (K)
1535
1536	c Local
1537	integer :: ik ! Loop variables
1538	integer :: indexice ! Index of the ice
1539
1540	c Code:
1541	indexice = -1
1542	if(ice_depth.lt.mlayer(0)) then
1543	indexice = 0.
1544	else
1545	do ik = 0,nsoil-2 ! go through all the layers to find the ice locations
1546	if((mlayer(ik).le.ice_depth).and.
1547	& (mlayer(ik+1).gt.ice_depth)) then
1548	indexice = ik+1
1549	exit
1550	endif
1551	enddo
1552	endif
1553
1554	if(indexice.lt.0) then
1555	call abort_physic("vdifc - compute Tice",
1556	& "subsurface ice is below the last soil layer",1)
1557	else
1558	if(indexice .ge. 1) then ! Linear inteprolation between soil temperature
1559	Tice = (ptsoil(indexice)-ptsoil(indexice+1))
1560	& /(mlayer(indexice-1)-mlayer(indexice))
1561	& *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1)
1562	else ! Linear inteprolation between the 1st soil temperature and the surface temperature
1563	Tice = (ptsoil(1) - ptsurf)/mlayer(0)
1564	& *(ice_depth-mlayer(0)) + ptsoil(1)
1565	endif ! index ice >=0
1566	endif !indexice <0
1567
1568
1569	END SUBROUTINE compute_Tice
1570	END MODULE vdifc_mod

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

Download in other formats:

Original Format