Context Navigation

vdifc_mod.F @ 2909

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

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

Download in other formats:

Original Format