Context Navigation

vdifc_mod.F @ 3400

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

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

File size: 61.5 KB

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

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

Download in other formats:

Original Format