Context Navigation

vdifc_mod.F @ 3167

Last change on this file since 3167 was 3167, checked in by llange, 11 months ago
Mars PCM Introducing the scheme from ATKE workshop to solve turbulent diffusion + surface layer parameterization. Works only if callatke = .true. in the callphys.def. By default, it is false and the model runs as usual with the yamada 2.5 closure Tuning of the several parameters for the ATKE in progress LL
File size: 58.8 KB

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

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

Download in other formats:

Original Format