Context Navigation

vdifc_mod.F @ 2538

Last change on this file since 2538 was 2531, checked in by emillour, 4 years ago
Mars GCM: Bug (wrong usage of surface pressure) fix in vdifc (bug was introduced in r2515) EM
File size: 38.4 KB

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

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

Download in other formats:

Original Format