Context Navigation

vdifc_mod.F @ 2461

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

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

Download in other formats:

Original Format