Context Navigation

turbdiff.F90 @ 3567

Last change on this file since 3567 was 2291, checked in by mlefevre, 5 years ago
MESOSCALE. Implementation of mesoscale architecture for Titan physics. See MESOSCALE flags.
File size: 16.8 KB

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

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

Download in other formats:

Original Format