Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

vdifc_mod.F

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

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

Download in other formats:

Original Format