Context Navigation

vdifc_mod.F @ 3404

Last change on this file since 3404 was 3404, checked in by llange, 3 months ago
Mars PCM Following previous commit, fixing a bug: more frost than what is a the surface could sublimes (when working with lag layer only). It is corrected, not in a elegant way... LL
File size: 61.6 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	if ((-zdqsdif_tot(ig)*subtimestep)
1124	& .ge.(zqsurf(ig))) then
1125	zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1126	endif
1127
1128	if((h2o_ice_depth(ig,islope).gt.0)
1129	& .and. (zqsurf(ig) .lt. tol_frost)
1130	& .and.lag_layer) then
1131	! 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
1132	! zqsurf(ig) = 0
1133	if(old_wsublimation_scheme) then
1134	resist(ig,islope)=h2o_ice_depth(ig,islope)
1135	& *zcdv(ig,islope)/d_coef(ig,islope)
1136	else
1137	resist(ig,islope)=h2o_ice_depth(ig,islope)
1138	& *zcdh(ig,islope)/d_coef(ig,islope)
1139	endif
1140
1141	zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice
1142	! Now we add the flux from the subsurface ice : rho Cd,h U(1/1+R) (q1-qsat_ssi(Tice))
1143
1144	call compute_Tice(nsoil, ptsoil(ig,:,islope),
1145	& ztsrf(ig),
1146	& h2o_ice_depth(ig,islope),
1147	& Tice(ig,islope)) ! compute ice temperature
1148
1149	call watersat(1,Tice(ig,islope),pplev(ig,1),
1150	& qsat_ssi(ig,islope))
1151
1152	qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1153	& *qsat_ssi(ig,islope)
1154
1155	! And now we solve correctly the system
1156	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1157	& zb(ig,2)*(1.-zd(ig,2)))
1158	zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
1159	& zb(ig,2)*zc(ig,2) +
1160	& (-zdqsdif_tot(ig)) *subtimestep)
1161	& * z1(ig)
1162	zd(ig,1)=zb(ig,1)*z1(ig)
1163	zq1temp(ig)=zc(ig,1)+ zd(ig,1)
1164	& *qsat_ssi(ig,islope)
1165
1166	! Flux from the subsurface
1167	if(old_wsublimation_scheme) then
1168	zdqsdif_ssi_atm(ig,islope)=rho(ig)
1169	& dryness(ig)zcdv(ig,islope)
1170	& *1/(1+resist(ig,islope))
1171	& *(zq1temp(ig)-qsat_ssi(ig,islope))
1172	else
1173	zdqsdif_ssi_atm(ig,islope)=rho(ig)*
1174	& dryness(ig) zcdh(ig,islope)
1175	& *1/(1+resist(ig,islope))
1176	& *(zq1temp(ig)-qsat_ssi(ig,islope))
1177	endif
1178
1179	else
1180	! No atm <-> subsurface exchange, we do it the usual way
1181	zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1182	z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
1183	zc(ig,1)=(za(ig,1)*
1184	& zq_tmp_vap(ig,1,igcm_h2o_vap)+
1185	& zb(ig,2)*zc(ig,2) +
1186	& (-zdqsdif_tot(ig)) subtimestep) z1(ig)
1187	zq1temp(ig)=zc(ig,1)
1188	endif ! Of subsurface <-> atmosphere exchange
1189	endif ! sublimating more than available frost & surface - frost exchange
1190	endif !if .not.watercaptag(ig)
1191
1192	! --------------------------------------------------------------------------------------------------------------------------------
1193	! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet),
1194	! we do not include their effect on the mass of surface frost.
1195	! --------------------------------------------------------------------------------------------------------------------------------
1196
1197	if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer
1198	& .and.(.not.adsorption_soil)) then
1199	! First case: still frost at the surface but no watercaptag
1200	if(((watercaptag(ig))).or.
1201	& (((-zdqsdif_tot(ig)*subtimestep)
1202	& .lt.(zqsurf(ig)))
1203	& .and. (zqsurf(ig).gt.tol_frost))) then
1204	! Still frost at the surface: we consider the possibility to have subsurface <-> frost exchange
1205	! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice))
1206	call compute_Tice(nsoil, ptsoil(ig,:,islope),
1207	& ztsrf(ig),
1208	& h2o_ice_depth(ig,islope),
1209	& Tice(ig,islope)) ! compute ice temperature
1210
1211	call watersat(1,Tice(ig,islope),pplev(ig,1),
1212	& qsat_ssi(ig,islope))
1213
1214	qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1215	& *qsat_ssi(ig,islope)
1216
1217	zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope)
1218	& /h2o_ice_depth(ig,islope)
1219	& rho(ig)dryness(ig)
1220	& *(qsat(ig)-qsat_ssi(ig,islope))
1221
1222	! Line to comment for now since we don't have a mass of subsurface frost
1223	! in our computation (otherwise, we would not conserve the H2O mass in
1224	! the system)
1225	if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope))
1226	& *subtimestep.ge.(zqsurf(ig))) then
1227	zdqsdif_ssi_frost(ig,islope) =
1228	& zqsurf(ig)/subtimestep
1229	& + zdqsdif_tot(ig)
1230	endif
1231
1232	zdqsdif_tot(ig) = zdqsdif_tot(ig) -
1233	& zdqsdif_ssi_frost(ig,islope)
1234	endif ! watercaptag or frost at the surface
1235	endif ! interaction frost <-> subsurface ice
1236
1237
1238	c Starting upward calculations for water :
1239	zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
1240	c Take into account the H2O latent heat impact on the surface temperature
1241	if (latentheat_surfwater) then
1242	lh=(2834.3-0.28*(ztsrf(ig)-To)-
1243	& 0.004(ztsrf(ig)-To)(ztsrf(ig)-To))*1.e+3
1244	zdtsrf(ig,islope)= zdqsdif_tot(ig)*lh
1245	& /pcapcal(ig,islope)
1246	endif ! (latentheat_surfwater) then
1247
1248	DO ilay=2,nlay
1249	zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)
1250	& *zq_tmp_vap(ig,ilay-1,iq)
1251	ENDDO
1252	c Subtimestep water budget :
1253	ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope)
1254	& + zdtsrf(ig,islope))*subtimestep
1255	zqsurf(ig)= zqsurf(ig)+(
1256	& zdqsdif_tot(ig))*subtimestep
1257
1258	if (zqsurf(ig)<0 .and.(.not.watercaptag(ig))) then
1259	zqsurf(ig)=0
1260	endif
1261	zdqsdif_ssi_atm_tot(ig,islope) =
1262	& zdqsdif_ssi_atm_tot(ig,islope)
1263	& + zdqsdif_ssi_atm(ig,islope)*subtimestep
1264	zdqsdif_ssi_frost_tot(ig,islope) =
1265	& zdqsdif_ssi_frost_tot(ig,islope)
1266	& +zdqsdif_ssi_frost(ig,islope)*subtimestep
1267	c Monitoring instantaneous latent heat flux in W.m-2 :
1268	zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
1269	& (zdtsrf(ig,islope)*pcapcal(ig,islope))
1270	& *subtimestep
1271
1272	c We ensure that surface temperature can't rise above the solid domain if there
1273	c is still ice on the surface (oldschool)
1274	if(zqsurf(ig)
1275	& +zdqsdif_tot(ig)*subtimestep
1276	& .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
1277	zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
1278	& (To-ztsrf(ig))/subtimestep) ! ice melt case
1279	endif
1280
1281	c End of the subtimestep
1282	ENDDO ! tsub=1,nsubtimestep
1283
1284	c Integration of subtimestep temp and water budget :
1285	c (btw could also compute the post timestep temp and ice
1286	c by simply adding the subtimestep trend instead of this)
1287	surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep
1288	pdtsrf(ig,islope)= (ztsrf(ig) -
1289	& ptsrf(ig,islope))/ptimestep
1290	pdqsdif(ig,igcm_h2o_ice,islope)=
1291	& (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/
1292	& cos(pi*def_slope_mean(islope)/180.))
1293	& /ptimestep
1294
1295	zdqsdif_ssi_atm_tot(ig,islope)=
1296	& zdqsdif_ssi_atm_tot(ig,islope)/ptimestep
1297	zdqsdif_ssi_frost_tot(ig,islope)=
1298	& zdqsdif_ssi_frost_tot(ig,islope)/ptimestep
1299	zdqsdif_ssi_tot(ig,islope) =
1300	& zdqsdif_ssi_atm_tot(ig,islope)
1301	& + zdqsdif_ssi_frost_tot(ig,islope)
1302	c if subliming more than qsurf(ice) and on watercaptag, water
1303	c sublimates from watercap reservoir (dwatercap_dif is <0)
1304	if(watercaptag(ig)) then
1305	if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep)
1306	& .gt.(pqsurf(ig,igcm_h2o_ice,islope)
1307	& /cos(pi*def_slope_mean(islope)/180.))) then
1308	dwatercap_dif(ig,islope)=
1309	& pdqsdif(ig,igcm_h2o_ice,islope)+
1310	& (pqsurf(ig,igcm_h2o_ice,islope) /
1311	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1312	pdqsdif(ig,igcm_h2o_ice,islope)=
1313	& - (pqsurf(ig,igcm_h2o_ice,islope)/
1314	& cos(pi*def_slope_mean(islope)/180.))/ptimestep
1315	endif! ((-pdqsdif(ig)*ptimestep)
1316	endif !(watercaptag(ig)) then
1317	zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:)
1318	ENDDO ! of DO ig=1,ngrid
1319	ENDDO ! islope
1320	c Some grid box averages: interface with the atmosphere
1321	DO ig = 1,ngrid
1322	DO ilay = 1,nlay
1323	zq(ig,ilay,iq) = 0.
1324	DO islope = 1,nslope
1325	zq(ig,ilay,iq) = zq(ig,ilay,iq) +
1326	$ zq_slope_vap(ig,ilay,iq,islope) *
1327	$ subslope_dist(ig,islope)
1328	ENDDO
1329	ENDDO
1330	ENDDO
1331	! Recompute values in kg/m^2 slopped
1332	DO ig = 1,ngrid
1333	DO islope = 1,nslope
1334	pdqsdif(ig,igcm_h2o_ice,islope) =
1335	& pdqsdif(ig,igcm_h2o_ice,islope)
1336	& * cos(pi*def_slope_mean(islope)/180.)
1337
1338	dwatercap_dif(ig,islope) =
1339	& dwatercap_dif(ig,islope)
1340	& * cos(pi*def_slope_mean(islope)/180.)
1341	ENDDO
1342	ENDDO
1343
1344	END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
1345
1346	c --------- end of h2o_vap ----------------------------
1347
1348	c --------- hdo_vap -----------------------------------
1349
1350	c hdo_ice has already been with along h2o_ice
1351	c amongst "normal" tracers (ie not h2o_vap)
1352
1353	if (hdo.and.(iq.eq.igcm_hdo_vap)) then
1354	zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1355	zb(1:ngrid,1)=0
1356
1357	DO ig=1,ngrid
1358	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1359	zc(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
1360	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1361	ENDDO
1362
1363	DO ilay=nlay-1,2,-1
1364	DO ig=1,ngrid
1365	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1366	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1367	zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
1368	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
1369	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1370	ENDDO
1371	ENDDO
1372	hdoflux_meshavg(:) = 0.
1373	DO islope = 1,nslope
1374
1375	pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope)
1376	& /cos(pi*def_slope_mean(islope)/180.)
1377
1378	call watersat(ngrid,pdtsrf(:,islope)*ptimestep +
1379	& ptsrf(:,islope),pplev(:,1),qsat_tmp)
1380
1381	CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
1382	& zt,pplay,zq,pqsurf(:,:,islope),
1383	& saved_h2o_vap,qsat_tmp,
1384	& pdqsdif_tmphdo,
1385	& dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.),
1386	& hdoflux(:,islope))
1387
1388	pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) *
1389	& cos(pi*def_slope_mean(islope)/180.)
1390	DO ig = 1,ngrid
1391	hdoflux_meshavg(ig) = hdoflux_meshavg(ig) +
1392	& hdoflux(ig,islope)*subslope_dist(ig,islope)
1393
1394	ENDDO !ig = 1,ngrid
1395	ENDDO !islope = 1,nslope
1396
1397	DO ig=1,ngrid
1398	z1(ig)=1./(za(ig,1)+zb(ig,1)+
1399	$ zb(ig,2)*(1.-zd(ig,2)))
1400	zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
1401	$ zb(ig,2)*zc(ig,2) +
1402	$ (-hdoflux_meshavg(ig)) ptimestep) z1(ig) !tracer flux from surface
1403	ENDDO
1404
1405	DO ig=1,ngrid
1406	zq(ig,1,iq)=zc(ig,1)
1407	DO ilay=2,nlay
1408	zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
1409	ENDDO
1410	ENDDO
1411	endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
1412
1413	c --------- end of hdo ----------------------------
1414
1415	enddo ! of do iq=1,nq
1416
1417	c --------- end of tracers ----------------------------
1418
1419	call write_output("surf_h2o_lh",
1420	& "Ground ice latent heat flux",
1421	& "W.m-2",surf_h2o_lh(:,iflat))
1422
1423	call write_output('zdqsdif_ssi_frost_tot',
1424	& 'Flux between frost and subsurface ice','kg.m-2.s-1',
1425	& zdqsdif_ssi_frost_tot(:,iflat))
1426
1427	call write_output('zdqsdif_ssi_atm_tot',
1428	& 'Flux between atmosphere and subsurface ice','kg.m-2.s-1',
1429	& zdqsdif_ssi_atm_tot(:,iflat))
1430
1431	call write_output('zdqsdif_ssi_tot',
1432	& 'Total flux echange with subsurface ice','kg.m-2.s-1',
1433	& zdqsdif_ssi_tot(:,iflat))
1434
1435
1436	C Diagnostic output for HDO
1437	! if (hdo) then
1438	! CALL write_output('hdoflux',
1439	! & 'hdoflux',
1440	! & ' ',hdoflux_meshavg(:))
1441	! CALL write_output('h2oflux',
1442	! & 'h2oflux',
1443	! & ' ',h2oflux(:))
1444	! endif
1445
1446	c-----------------------------------------------------------------------
1447	c 8. calcul final des tendances de la diffusion verticale
1448	c ----------------------------------------------------
1449
1450	DO ilev = 1, nlay
1451	DO ig=1,ngrid
1452	pdudif(ig,ilev)=( zu(ig,ilev)-
1453	$ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
1454	pdvdif(ig,ilev)=( zv(ig,ilev)-
1455	$ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
1456	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
1457	$ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
1458	pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
1459	ENDDO
1460	ENDDO
1461
1462	pdqdif(1:ngrid,1:nlay,1:nq)=(zq(1:ngrid,1:nlay,1:nq)-
1463	& (pq(1:ngrid,1:nlay,1:nq)
1464	& +pdqfi(1:ngrid,1:nlay,1:nq)
1465	& *ptimestep))/ptimestep
1466
1467	c ** diagnostique final
1468	c ------------------
1469
1470	IF(lecrit) THEN
1471	PRINT*,'In vdif'
1472	PRINT*,'Ts (t) and Ts (t+st)'
1473	WRITE(*,'(a10,3a15)')
1474	s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
1475	PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1)
1476	DO ilev=1,nlay
1477	WRITE(*,'(4f15.7)')
1478	s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
1479	s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
1480
1481	ENDDO
1482	ENDIF
1483
1484	END SUBROUTINE vdifc
1485
1486	c====================================
1487
1488	SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
1489	$ dtmax,watercaptag,ntsub)
1490
1491	c Pas de temps adaptatif en estimant le taux de sublimation
1492	c et en adaptant avec un critere "dtmax" du chauffage a accomoder
1493	c dtmax est regle empiriquement (pour l'instant) a 0.5 K
1494	use surfdat_h, only: frost_albedo_threshold
1495
1496	integer,intent(in) :: naersize
1497	real,intent(in) :: dtsurf(naersize)
1498	real,intent(in) :: qsurf(naersize)
1499	logical,intent(in) :: watercaptag(naersize)
1500	real,intent(in) :: ptimestep
1501	real,intent(in) :: dtmax
1502	real :: ztsub(naersize)
1503	integer :: i
1504	integer,intent(out) :: ntsub(naersize)
1505
1506	do i=1,naersize
1507	c If not on permanent ice or thin frost layer :
1508	c do not use a subtimestep (ntsub=1)
1509	if ((qsurf(i).lt.frost_albedo_threshold).and.
1510	& (.not.watercaptag(i))) then
1511	ntsub(i) = 1
1512	else
1513	ztsub(i) = ptimestep * dtsurf(i) / dtmax
1514	ntsub(i) = ceiling(abs(ztsub(i)))
1515	endif ! (qsurf(i).eq.0) then
1516	c
1517	c write(78,), dtsurfptimestep, dtsurf, ntsub
1518	enddo! 1=1,ngrid
1519
1520
1521
1522	END SUBROUTINE make_tsub
1523
1524
1525	c====================================
1526
1527	SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice)
1528
1529	c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
1530	use comsoil_h, only: layer, mlayer
1531
1532	implicit none
1533	integer,intent(in) :: nsoil ! Number of soil layers
1534	real,intent(in) :: ptsoil(nsoil) ! Soil temperature (K)
1535	real,intent(in) :: ptsurf ! Soil temperature (K)
1536	real,intent(in) :: ice_depth ! Ice depth (m)
1537	real,intent(out) :: Tice ! Ice temperatures (K)
1538
1539	c Local
1540	integer :: ik ! Loop variables
1541	integer :: indexice ! Index of the ice
1542
1543	c Code:
1544	indexice = -1
1545	if(ice_depth.lt.mlayer(0)) then
1546	indexice = 0.
1547	else
1548	do ik = 0,nsoil-2 ! go through all the layers to find the ice locations
1549	if((mlayer(ik).le.ice_depth).and.
1550	& (mlayer(ik+1).gt.ice_depth)) then
1551	indexice = ik+1
1552	exit
1553	endif
1554	enddo
1555	endif
1556
1557	if(indexice.lt.0) then
1558	call abort_physic("vdifc - compute Tice",
1559	& "subsurface ice is below the last soil layer",1)
1560	else
1561	if(indexice .ge. 1) then ! Linear inteprolation between soil temperature
1562	Tice = (ptsoil(indexice)-ptsoil(indexice+1))
1563	& /(mlayer(indexice-1)-mlayer(indexice))
1564	& *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1)
1565	else ! Linear inteprolation between the 1st soil temperature and the surface temperature
1566	Tice = (ptsoil(1) - ptsurf)/mlayer(0)
1567	& *(ice_depth-mlayer(0)) + ptsoil(1)
1568	endif ! index ice >=0
1569	endif !indexice <0
1570
1571
1572	END SUBROUTINE compute_Tice
1573	END MODULE vdifc_mod

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

Download in other formats:

Original Format