Context Navigation

vdifc_mod.F @ 3118

Last change on this file since 3118 was 3115, checked in by llange, 13 months ago

Mars PCM

Introducing the possibily to compute water adsorption / desorption -routine soilwater.F90) -. FOR NOW IT WORKS ONLY IN 1D, DON'T TEST IN 3D (by default, adsorption is not called, if not specified in the callphys.def). By default, when using the adsorption, the number of subtimestep used in the water subimation scheme is 1 (otherwise, it crashes)

New grid for the soil layers (better resolution) as it is needed to solve the diffusion equations. It does not increase the CPU time.

PYM & LL

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

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

Download in other formats:

Original Format