Context Navigation

vdifc_mod.F @ 3026

Last change on this file since 3026 was 3008, checked in by emillour, 18 months ago
Mars PCM: Some code cleanup around microphysics. Turn microphys.h into module microphys_h.F90, and while at it also turn nuclea.F, growthrate.F90 and massflowrateco2.F90 into modules. EM
File size: 43.7 KB

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

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

Download in other formats:

Original Format