Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

vdifc_mod.F

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

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

Download in other formats:

Original Format