Context Navigation

turbdiff.F90 @ 1680

Last change on this file since 1680 was 1647, checked in by jvatant, 8 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

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

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

Download in other formats:

Original Format