Context Navigation

turbdiff_mod.F90 @ 3608

Last change on this file since 3608 was 3228, checked in by tbertrand, 12 months ago
Pluto GCM: Cleaning code and adding extra options for newstart.F
File size: 18.3 KB

Rev	Line
[3184]	1	module turbdiff_mod
[3197]	2
[3184]	3	implicit none
[3197]	4
[3184]	5	contains
[3197]	6
[3184]	7	subroutine turbdiff(ngrid,nlay,nq, &
[3197]	8	ptimestep,pcapcal, &
[3184]	9	pplay,pplev,pzlay,pzlev,pz0, &
	10	pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf, &
	11	pdtfi,pdqfi,pfluxsrf, &
	12	Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, &
	13	pdqdif,pdqevap,pdqsdif,flux_u,flux_v)
	14
	15	use radcommon_h, only : sigma, glat
	16	use surfdat_h, only: dryness
	17	use comcstfi_mod, only: rcp, g, r, cpp
	18	use callkeys_mod, only: tracer,nosurf,kmixmin
	19	use turb_mod, only : ustar
	20	#ifdef MESOSCALE
	21	use comm_wrf, only : comm_LATENT_HF
	22	#endif
	23	implicit none
	24
	25	!==================================================================
[3197]	26	!
[3184]	27	! Purpose
	28	! -------
	29	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
[3197]	30	!
[3184]	31	! Implicit scheme
	32	! We start by adding to variables x the physical tendencies
	33	! already computed. We resolve the equation:
	34	!
	35	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
[3197]	36	!
[3184]	37	! Authors
[3197]	38	! -------
[3184]	39	! F. Hourdin, F. Forget, R. Fournier (199X)
	40	! R. Wordsworth, B. Charnay (2010)
	41	! J. Leconte (2012): To f90
	42	! - Rewritten the diffusion scheme to conserve total enthalpy
	43	! by accounting for dissipation of turbulent kinetic energy.
	44	! - Accounting for lost mean flow kinetic energy should come soon.
[3197]	45	!
[3184]	46	!==================================================================
	47
	48	!-----------------------------------------------------------------------
	49	! declarations
	50	! ------------
	51
	52	! arguments
	53	! ---------
	54	INTEGER,INTENT(IN) :: ngrid
	55	INTEGER,INTENT(IN) :: nlay
	56	REAL,INTENT(IN) :: ptimestep
	57	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	58	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	59	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
	60	REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
	61	REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
	62	REAL,INTENT(IN) :: pemis(ngrid)
	63	REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
	64	REAL,INTENT(IN) :: pfluxsrf(ngrid)
	65	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
	66	REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
	67	REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
	68	REAL,INTENT(OUT) :: sensibFlux(ngrid)
	69	REAL,INTENT(IN) :: pcapcal(ngrid)
	70	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
	71	REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
	72
	73	REAL,INTENT(IN) :: pz0
	74
	75	! Tracers
	76	! --------
[3197]	77	integer,intent(in) :: nq
[3184]	78	real,intent(in) :: pqsurf(ngrid,nq)
[3197]	79	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	80	real,intent(out) :: pdqdif(ngrid,nlay,nq)
	81	real,intent(out) :: pdqsdif(ngrid,nq)
	82
[3184]	83	! local
	84	! -----
	85	integer ilev,ig,ilay,nlev
	86
	87	REAL z4st,zdplanck(ngrid)
	88	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
	89	REAL zcdv(ngrid),zcdh(ngrid)
	90	REAL zcdv_true(ngrid),zcdh_true(ngrid)
	91	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	92	REAL zh(ngrid,nlay),zt(ngrid,nlay)
	93	REAL ztsrf(ngrid)
	94	REAL z1(ngrid),z2(ngrid)
	95	REAL zmass(ngrid,nlay)
	96	REAL zfluxv(ngrid,nlay),zfluxt(ngrid,nlay),zfluxq(ngrid,nlay)
	97	REAL zb0(ngrid,nlay)
	98	REAL zExner(ngrid,nlay),zovExner(ngrid,nlay)
	99	REAL zcv(ngrid,nlay),zdv(ngrid,nlay) !inversion coefficient for winds
	100	REAL zct(ngrid,nlay),zdt(ngrid,nlay) !inversion coefficient for temperature
	101	REAL zcq(ngrid,nlay),zdq(ngrid,nlay) !inversion coefficient for tracers
	102	REAL zcst1
	103	REAL zu2!, a
	104	REAL zcq0(ngrid),zdq0(ngrid)
	105	REAL zx_alf1(ngrid),zx_alf2(ngrid)
	106	! 1D eddy diffusion coefficient
	107	REAL kzz_eddy(nlay)
	108	REAL pmin_kzz
	109
	110	LOGICAL,SAVE :: firstcall=.true.
	111	!$OMP THREADPRIVATE(firstcall)
[3197]	112
[3184]	113	! Tracers
	114	! -------
	115	INTEGER iq
	116	REAL zq(ngrid,nlay,nq)
	117	REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
	118	REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
	119	REAL zdmassevap(ngrid)
	120	REAL rho(ngrid) ! near-surface air density
	121
	122	! Variables added for implicit latent heat inclusion
	123	! --------------------------------------------------
	124	real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
	125
	126	integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
	127	!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
	128
	129	real, parameter :: karman=0.4
	130	real cd0, roughratio
	131
[3197]	132	real dqsdif_total(ngrid)
	133	real zq0(ngrid)
[3184]	134
	135
	136	! Coherence test
	137	! --------------
	138
	139	IF (firstcall) THEN
[3228]	140	ivap=1
[3198]	141	iliq=0
	142	iliq_surf=0
	143	iice_surf=0 ! simply to make the code legible
[3184]	144	sensibFlux(:)=0.
	145	firstcall=.false.
	146	ENDIF
	147
	148	!-----------------------------------------------------------------------
	149	! 1. Initialisation
	150	! -----------------
	151
	152	nlev=nlay+1
	153
	154	! Calculate rhodz, (P/Ps)(R/cp) and dtrho/dz=dtrho*2 g/dp
	155	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
	156	! ---------------------------------------------
	157
	158	DO ilay=1,nlay
	159	DO ig=1,ngrid
	160	zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
	161	zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
	162	zovExner(ig,ilay)=1./ppopsk(ig,ilay)
	163	ENDDO
	164	ENDDO
	165
	166	zcst1=4.gptimestep/(R*R)
	167	DO ilev=2,nlev-1
	168	DO ig=1,ngrid
	169	zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev))
	170	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev))
	171	ENDDO
	172	ENDDO
	173	DO ig=1,ngrid
	174	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
	175	ENDDO
	176	dqsdif_total(:)=0.0
	177
	178	!-----------------------------------------------------------------------
	179	! 2. Add the physical tendencies computed so far
	180	! ----------------------------------------------
	181
	182	DO ilev=1,nlay
	183	DO ig=1,ngrid
	184	zu(ig,ilev)=pu(ig,ilev)
	185	zv(ig,ilev)=pv(ig,ilev)
	186	zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
	187	zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
	188	ENDDO
	189	ENDDO
	190	if(tracer) then
	191	DO iq =1, nq
	192	DO ilev=1,nlay
	193	DO ig=1,ngrid
	194	zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep
	195	ENDDO
	196	ENDDO
	197	ENDDO
	198	end if
	199
	200	!-----------------------------------------------------------------------
	201	! 3. Turbulence scheme
	202	! --------------------
	203	!
	204	! Source of turbulent kinetic energy at the surface
[3197]	205	! -------------------------------------------------
[3184]	206	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
	207
	208	DO ig=1,ngrid
	209	roughratio = 1. + pzlay(ig,1)/pz0
	210	cd0 = karman/log(roughratio)
	211	cd0 = cd0*cd0
	212	zcdv_true(ig) = cd0
	213	zcdh_true(ig) = cd0
	214	if(nosurf)then
	215	zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
[3197]	216	zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
[3184]	217	endif
	218	ENDDO
	219
	220	DO ig=1,ngrid
	221	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	222	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	223	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
	224	ENDDO
	225
	226	! Turbulent diffusion coefficients in the boundary layer
[3197]	227	! ------------------------------------------------------
[3184]	228
[3197]	229	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
	230
[3184]	231	! Adding eddy mixing to mimic 3D general circulation in 1D
	232	! R. Wordsworth & F. Forget (2010)
	233	if ((ngrid.eq.1)) then
	234	! kmixmin minimum eddy mix coeff in 1D
	235	! set up in inifis_mod.F90 - default value 1.0e-2
	236	do ilev=1,nlay
	237
	238	! Here to code your specific eddy mix coeff in 1D
	239	! Earth example that can be uncommented below
	240	! -------------------------------------------------
	241	! ====== Earth kzz from Zahnle et al. 2006 ======
	242	! -------------------------------------------------
	243	! if(pzlev(1,ilev).le.11.0e3) then
	244	! kzz_eddy(ilev)=10.0
	245	! pmin_kzz=pplev(1,ilev)exp((pzlev(1,ilev)-11.0e3)g/(r*zt(1,ilev)))
	246	! else
	247	! kzz_eddy(ilev)=0.1(pplev(1,ilev)/pmin_kzz)*(-0.5)
	248	! kzz_eddy(ilev)=min(kzz_eddy(ilev),100.0)
	249	! endif
	250	! do ig=1,ngrid
	251	! zkh(ig,ilev) = max(kzz_eddy(ilev),zkh(ig,ilev))
[3197]	252	! zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev))
[3184]	253	! end do
	254
	255	do ig=1,ngrid
	256	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
[3197]	257	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
[3184]	258	end do
	259	end do
	260	end if
	261
	262	!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
	263	DO ig=1,ngrid
	264	zkv(ig,1)=zcdv(ig)
	265	ENDDO
	266	!we treat only winds, energy and tracers coefficients will be computed with upadted winds
	267	!JL12 calculate the flux coefficients (tables multiplied elements by elements)
	268	zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
[3197]	269
[3184]	270	!-----------------------------------------------------------------------
	271	! 4. Implicit inversion of u
	272	! --------------------------
	273
	274	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	275	! avec
	276	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	277	! et
	278	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	279	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	280	! et /zkv/ = Ku
	281
	282	DO ig=1,ngrid
	283	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
	284	zcv(ig,nlay)=zmass(ig,nlay)zu(ig,nlay)z1(ig)
	285	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
	286	ENDDO
	287
	288	DO ilay=nlay-1,1,-1
	289	DO ig=1,ngrid
	290	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
	291	zcv(ig,ilay)=(zmass(ig,ilay)zu(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
	292	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
	293	ENDDO
	294	ENDDO
	295
	296	DO ig=1,ngrid
	297	zu(ig,1)=zcv(ig,1)
	298	ENDDO
	299	DO ilay=2,nlay
	300	DO ig=1,ngrid
	301	zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1)
	302	ENDDO
	303	ENDDO
	304
	305	!-----------------------------------------------------------------------
	306	! 5. Implicit inversion of v
	307	! --------------------------
	308
	309	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	310	! avec
	311	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	312	! et
	313	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	314	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	315	! et /zkv/ = Kv
	316
	317	DO ig=1,ngrid
[3197]	318	z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
[3184]	319	zcv(ig,nlay)=zmass(ig,nlay)zv(ig,nlay)z1(ig)
	320	zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
	321	ENDDO
	322
	323	DO ilay=nlay-1,1,-1
	324	DO ig=1,ngrid
	325	z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1)))
	326	zcv(ig,ilay)=(zmass(ig,ilay)zv(ig,ilay)+zfluxv(ig,ilay+1)zcv(ig,ilay+1))*z1(ig)
	327	zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig)
	328	ENDDO
	329	ENDDO
	330
	331	DO ig=1,ngrid
	332	zv(ig,1)=zcv(ig,1)
	333	ENDDO
	334	DO ilay=2,nlay
	335	DO ig=1,ngrid
	336	zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1)
	337	ENDDO
	338	ENDDO
	339
	340	! Calcul of wind stress
	341
	342	DO ig=1,ngrid
	343	flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
	344	flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
	345	ENDDO
	346
	347
	348	!----------------------------------------------------------------------------
	349	! 6. Implicit inversion of h, not forgetting the coupling with the ground
	350
	351	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	352	! avec
	353	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	354	! et
	355	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	356	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
	357	! et /zkh/ = Kh
	358
	359	! Using the wind modified by friction for lifting and sublimation
	360	! ---------------------------------------------------------------
	361	DO ig=1,ngrid
	362	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	363	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
	364	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
	365	zkh(ig,1)= zcdh(ig)
	366	ustar(ig)=sqrt(zcdv_true(ig))*sqrt(zu2)
	367	ENDDO
	368
	369
	370	! JL12 calculate the flux coefficients (tables multiplied elements by elements)
	371	! ---------------------------------------------------------------
	372	zfluxq(1:ngrid,1:nlay)=zkh(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay) !JL12 we save zfluxq which doesn't need the Exner factor
	373	zfluxt(1:ngrid,1:nlay)=zfluxq(1:ngrid,1:nlay)*zExner(1:ngrid,1:nlay)
	374
	375	DO ig=1,ngrid
	376	z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay))
	377	zct(ig,nlay)=zmass(ig,nlay)zt(ig,nlay)z1(ig)
	378	zdt(ig,nlay)=zfluxt(ig,nlay)zovExner(ig,nlay-1)z1(ig)
	379	ENDDO
	380
	381	DO ilay=nlay-1,2,-1
	382	DO ig=1,ngrid
	383	z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ &
	384	zfluxt(ig,ilay+1)(zovExner(ig,ilay)-zdt(ig,ilay+1)zovExner(ig,ilay+1)))
	385	zct(ig,ilay)=(zmass(ig,ilay)zt(ig,ilay)+zfluxt(ig,ilay+1)zct(ig,ilay+1)zovExner(ig,ilay+1))z1(ig)
	386	zdt(ig,ilay)=zfluxt(ig,ilay)z1(ig)zovExner(ig,ilay-1)
	387	ENDDO
	388	ENDDO
	389
	390	!JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there
	391	DO ig=1,ngrid
	392	z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ &
	393	zfluxt(ig,2)(zovExner(ig,1)-zdt(ig,2)zovExner(ig,2)))
	394	zct(ig,1)=(zmass(ig,1)zt(ig,1)+zfluxt(ig,2)zct(ig,2)zovExner(ig,2))z1(ig)
	395	zdt(ig,1)=zfluxt(ig,1)*z1(ig)
	396	ENDDO
	397
	398
	399	! Calculate (d Planck / dT) at the interface temperature
	400	! ------------------------------------------------------
	401
	402	z4st=4.0sigmaptimestep
	403	DO ig=1,ngrid
	404	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
	405	ENDDO
	406
	407	! Calculate temperature tendency at the interface (dry case)
	408	! ----------------------------------------------------------
	409	! Sum of fluxes at interface at time t + \delta t gives change in T:
	410	! radiative fluxes
	411	! turbulent convective (sensible) heat flux
	412	! flux (if any) from subsurface
	413
	414	! if(.not.water) then
	415
	416	DO ig=1,ngrid
	417	z1(ig)=pcapcal(ig)ptsrf(ig)+cppzfluxt(ig,1)zct(ig,1)zovExner(ig,1) &
[3197]	418	+ pfluxsrf(ig)ptimestep + zdplanck(ig)ptsrf(ig)
	419	z2(ig) = pcapcal(ig)+zdplanck(ig)+cppzfluxt(ig,1)(1.-zovExner(ig,1)*zdt(ig,1))
[3184]	420	ztsrf(ig) = z1(ig) / z2(ig)
	421	pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
	422	zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
	423	ENDDO
[3197]	424	! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
[3184]	425
	426
	427	! Recalculate temperature to top of atmosphere, starting from ground
	428	! ------------------------------------------------------------------
	429
	430	DO ilay=2,nlay
	431	DO ig=1,ngrid
	432	zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
	433	ENDDO
	434	ENDDO
	435
	436	! endif ! not water
	437
	438	!-----------------------------------------------------------------------
	439	! TRACERS (no vapour)
	440	! -------
	441
	442	if(tracer) then
	443
	444	! Calculate vertical flux from the bottom to the first layer (dust)
	445	! -----------------------------------------------------------------
	446	do ig=1,ngrid
	447	rho(ig) = zb0(ig,1) /ptimestep
	448	end do
	449
	450	pdqsdif(:,:)=0.
	451
	452	! Implicit inversion of q
	453	! -----------------------
[3197]	454	do iq=1,nq
[3184]	455
	456	if (iq.ne.ivap) then
	457
	458	DO ig=1,ngrid
	459	z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
	460	zcq(ig,nlay)=zmass(ig,nlay)zq(ig,nlay,iq)z1(ig)
	461	zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
[3197]	462	ENDDO
	463
[3184]	464	DO ilay=nlay-1,2,-1
	465	DO ig=1,ngrid
	466	z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
	467	zcq(ig,ilay)=(zmass(ig,ilay)zq(ig,ilay,iq)+zfluxq(ig,ilay+1)zcq(ig,ilay+1))*z1(ig)
	468	zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
	469	ENDDO
	470	ENDDO
	471
[3197]	472	do ig=1,ngrid
	473	z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
	474	zcq(ig,1)=(zmass(ig,1)zq(ig,1,iq)+zfluxq(ig,2)zcq(ig,2)+(-pdqsdif(ig,iq))ptimestep)z1(ig)
	475	! tracer flux from surface
	476	! currently pdqsdif always zero here,
	477	! so last line is superfluous
	478	enddo
	479
[3184]	480	! Starting upward calculations for simple tracer mixing (e.g., dust)
	481	do ig=1,ngrid
	482	zq(ig,1,iq)=zcq(ig,1)
	483	end do
	484
	485	do ilay=2,nlay
	486	do ig=1,ngrid
	487	zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
	488	end do
	489	end do
	490
	491	endif ! if (iq.ne.ivap)
	492	end do ! of do iq=1,nq
	493
	494	! Revmoed (AF24): Calculate temperature tendency including latent heat term
	495	! and assuming an infinite source of water on the ground
	496	! ------------------------------------------------------------------
[3197]	497
[3184]	498	endif ! tracer
	499
	500
	501	!-----------------------------------------------------------------------
	502	! 8. Final calculation of the vertical diffusion tendencies
	503	! -----------------------------------------------------------------
	504
	505	do ilev = 1, nlay
	506	do ig=1,ngrid
	507	pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
	508	pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
	509	pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
	510	enddo
	511	enddo
[3197]	512
[3184]	513	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
	514	sensibFlux(ig)=cppzfluxt(ig,1)/ptimestep(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
	515	ENDDO
	516
	517	if (tracer) then
	518	do iq = 1, nq
	519	do ilev = 1, nlay
	520	do ig=1,ngrid
	521	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
	522	enddo
	523	enddo
	524	enddo
[3197]	525	endif
[3184]	526	end subroutine turbdiff
[3197]	527
[3184]	528	end module turbdiff_mod

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

Download in other formats:

Original Format