Context Navigation

vdifc_mod.F @ 3586

Last change on this file since 3586 was 3275, checked in by afalco, 11 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
File size: 17.1 KB

Rev	Line
[3184]	1	module vdifc_mod
[3195]	2
[3184]	3	implicit none
[3195]	4
[3184]	5	contains
[3195]	6
	7	subroutine vdifc(ngrid,nlay,nq,ppopsk,
	8	& ptimestep,pcapcal,lecrit,
[3184]	9	& pplay,pplev,pzlay,pzlev,pz0,
	10	& pu,pv,ph,pq,ptsrf,pemis,pqsurf,
	11	& pdhfi,pdqfi,pfluxsrf,
	12	& pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
	13	& pdqdif,pdqsdif)
	14
[3195]	15	! use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
[3184]	16	! & ,Psat_water, Lcpdqsat_water
	17	use radcommon_h, only : sigma
	18	use surfdat_h, only: dryness
	19	use comcstfi_mod, only: g, r, cpp, rcp
	20	use callkeys_mod, only: tracer,nosurf
	21
	22	implicit none
	23
	24	!==================================================================
[3195]	25	!
[3184]	26	! Purpose
	27	! -------
	28	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
[3195]	29	!
[3184]	30	! Implicit scheme
	31	! We start by adding to variables x the physical tendencies
	32	! already computed. We resolve the equation:
	33	!
	34	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
[3195]	35	!
[3184]	36	! Authors
[3195]	37	! -------
[3184]	38	! F. Hourdin, F. Forget, R. Fournier (199X)
	39	! R. Wordsworth, B. Charnay (2010)
[3195]	40	!
[3184]	41	!==================================================================
	42
	43	!-----------------------------------------------------------------------
	44	! declarations
	45	! ------------
	46
	47
	48	! arguments
	49	! ---------
	50	INTEGER,INTENT(IN) :: ngrid,nlay
	51	REAL,INTENT(IN) :: ptimestep
	52	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	53	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	54	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	55	REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
	56	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
	57	REAL,INTENT(IN) :: pfluxsrf(ngrid)
	58	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
	59	REAL,INTENT(OUT) :: pdhdif(ngrid,nlay)
	60	REAL,INTENT(OUT) :: pdtsrf(ngrid),sensibFlux(ngrid)
	61	REAL,INTENT(IN) :: pcapcal(ngrid)
	62	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
[3195]	63
[3184]	64	! Arguments added for condensation
	65	REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
	66	logical,intent(in) :: lecrit
	67	REAL,INTENT(IN) :: pz0
	68
	69	! Tracers
	70	! --------
[3195]	71	integer,intent(in) :: nq
[3184]	72	real,intent(in) :: pqsurf(ngrid,nq)
[3195]	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)
	76
[3184]	77	! local
	78	! -----
	79	integer ilev,ig,ilay,nlev
	80
	81	REAL z4st,zdplanck(ngrid)
	82	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
	83	REAL zcdv(ngrid),zcdh(ngrid)
	84	REAL zcdv_true(ngrid),zcdh_true(ngrid)
	85	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	86	REAL zh(ngrid,nlay)
	87	REAL ztsrf2(ngrid)
	88	REAL z1(ngrid),z2(ngrid)
	89	REAL za(ngrid,nlay),zb(ngrid,nlay)
	90	REAL zb0(ngrid,nlay)
	91	REAL zc(ngrid,nlay),zd(ngrid,nlay)
	92	REAL zcst1
	93	REAL zu2!, a
	94	REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
	95	REAL evap(ngrid)
	96	REAL zcq0(ngrid),zdq0(ngrid)
	97	REAL zx_alf1(ngrid),zx_alf2(ngrid)
	98
	99	LOGICAL firstcall
	100	SAVE firstcall
	101	!$OMP THREADPRIVATE(firstcall)
[3195]	102
[3184]	103	! variables added for N2 condensation
	104	! ------------------------------------
	105	REAL hh !, zhcond(ngrid,nlay)
	106	! REAL latcond,tcond1mb
	107	! REAL acond,bcond
	108	! SAVE acond,bcond
	109	!!$OMP THREADPRIVATE(acond,bcond)
	110	! DATA latcond,tcond1mb/5.9e5,136.27/
	111
	112	! Tracers
	113	! -------
	114	INTEGER iq
	115	REAL zq(ngrid,nlay,nq)
	116	REAL zq1temp(ngrid)
	117	REAL rho(ngrid) ! near-surface air density
	118	DATA firstcall/.true./
	119	REAL kmixmin
	120
	121	! Variables added for implicit latent heat inclusion
	122	! --------------------------------------------------
	123	real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
	124
	125	integer ivap, iice ! also make liq for clarity on surface...
	126	save ivap, iice
	127	!$OMP THREADPRIVATE(ivap,iice)
	128
	129	real, parameter :: karman=0.4
	130	real cd0, roughratio
	131
	132	logical forceWC
	133	real masse, Wtot, Wdiff
	134
[3195]	135	real dqsdif_total(ngrid)
	136	real zq0(ngrid)
[3184]	137
	138	forceWC=.true.
	139	! forceWC=.false.
	140
	141
	142	! Coherence test
	143	! --------------
	144
	145	IF (firstcall) THEN
	146	! To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
	147	! bcond=1./tcond1mb
	148	! acond=r/latcond
	149	! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
	150	! PRINT*,' acond,bcond',acond,bcond
	151
	152
	153	firstcall=.false.
	154	ENDIF
	155
	156	!-----------------------------------------------------------------------
	157	! 1. Initialisation
	158	! -----------------
	159
	160	nlev=nlay+1
	161
	162	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
	163	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
	164	! ---------------------------------------------
	165
	166	DO ilay=1,nlay
	167	DO ig=1,ngrid
	168	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
	169	ENDDO
	170	ENDDO
	171
	172	zcst1=4.gptimestep/(R*R)
	173	DO ilev=2,nlev-1
	174	DO ig=1,ngrid
	175	zb0(ig,ilev)=pplev(ig,ilev)*
	176	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
	177	s (ph(ig,ilev-1)+ph(ig,ilev))
	178	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
	179	s (pplay(ig,ilev-1)-pplay(ig,ilev))
	180	ENDDO
	181	ENDDO
	182	DO ig=1,ngrid
	183	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
	184	ENDDO
	185
	186	dqsdif_total(:)=0.0
	187
	188	!-----------------------------------------------------------------------
	189	! 2. Add the physical tendencies computed so far
	190	! ----------------------------------------------
	191
	192	DO ilev=1,nlay
	193	DO ig=1,ngrid
	194	zu(ig,ilev)=pu(ig,ilev)
	195	zv(ig,ilev)=pv(ig,ilev)
	196	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	197	ENDDO
	198	ENDDO
	199	if(tracer) then
	200	DO iq =1, nq
	201	DO ilev=1,nlay
	202	DO ig=1,ngrid
[3195]	203	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
[3184]	204	& pdqfi(ig,ilev,iq)*ptimestep
	205	ENDDO
	206	ENDDO
	207	ENDDO
	208	end if
	209
	210	!-----------------------------------------------------------------------
	211	! 3. Turbulence scheme
	212	! --------------------
	213	!
	214	! Source of turbulent kinetic energy at the surface
[3195]	215	! -------------------------------------------------
[3184]	216	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
	217
	218	DO ig=1,ngrid
	219	roughratio = 1.E+0 + pzlay(ig,1)/pz0
	220	cd0 = karman/log(roughratio)
	221	cd0 = cd0*cd0
	222	zcdv_true(ig) = cd0
	223	zcdh_true(ig) = cd0
	224	if (nosurf) then
	225	zcdv_true(ig) = 0. !! disable sensible momentum flux
	226	zcdh_true(ig) = 0. !! disable sensible heat flux
	227	endif
	228	ENDDO
	229
	230	DO ig=1,ngrid
	231	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	232	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	233	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
	234	ENDDO
	235
	236	! Turbulent diffusion coefficients in the boundary layer
[3195]	237	! ------------------------------------------------------
[3184]	238
	239	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
	240	& ,pu,pv,ph,zcdv_true
	241	& ,pq2,zkv,zkh)
	242
	243	! Adding eddy mixing to mimic 3D general circulation in 1D
	244	! R. Wordsworth & F. Forget (2010)
	245	if ((ngrid.eq.1)) then
	246	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
	247	do ilev=1,nlay
	248	do ig=1,ngrid
	249	!zkh(ig,ilev) = 1.0
	250	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
	251	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
	252	end do
	253	end do
	254	end if
	255
	256	!-----------------------------------------------------------------------
	257	! 4. Implicit inversion of u
	258	! --------------------------
	259
	260	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	261	! avec
	262	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	263	! et
	264	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	265	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	266	! et /zkv/ = Ku
[3195]	267
[3184]	268	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	269	CALL multipl(ngrid,zcdv,zb0,zb)
	270
	271	DO ig=1,ngrid
	272	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	273	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
	274	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	275	ENDDO
	276
	277	DO ilay=nlay-1,1,-1
	278	DO ig=1,ngrid
	279	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	280	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	281	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
	282	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	283	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	284	ENDDO
	285	ENDDO
	286
	287	DO ig=1,ngrid
	288	zu(ig,1)=zc(ig,1)
	289	ENDDO
	290	DO ilay=2,nlay
	291	DO ig=1,ngrid
	292	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	293	ENDDO
	294	ENDDO
	295
	296	!-----------------------------------------------------------------------
	297	! 5. Implicit inversion of v
	298	! --------------------------
	299
	300	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	301	! avec
	302	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	303	! et
	304	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	305	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	306	! et /zkv/ = Kv
	307
	308	DO ig=1,ngrid
	309	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	310	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	311	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	312	ENDDO
	313
	314	DO ilay=nlay-1,1,-1
	315	DO ig=1,ngrid
	316	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	317	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	318	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
	319	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	320	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	321	ENDDO
	322	ENDDO
	323
	324	DO ig=1,ngrid
	325	zv(ig,1)=zc(ig,1)
	326	ENDDO
	327	DO ilay=2,nlay
	328	DO ig=1,ngrid
	329	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	330	ENDDO
	331	ENDDO
	332
	333	!----------------------------------------------------------------------------
	334	! 6. Implicit inversion of h, not forgetting the coupling with the ground
	335
	336	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	337	! avec
	338	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	339	! et
	340	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	341	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
	342	! et /zkh/ = Kh
	343
	344	! Using the wind modified by friction for lifting and sublimation
	345	! ---------------------------------------------------------------
	346	DO ig=1,ngrid
	347	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	348	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
	349	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
	350	ENDDO
	351
	352	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	353	CALL multipl(ngrid,zcdh,zb0,zb)
	354
	355	DO ig=1,ngrid
	356	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	357	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
	358	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	359	ENDDO
	360
	361	DO ilay=nlay-1,2,-1
	362	DO ig=1,ngrid
	363	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	364	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
	365	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
	366	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
	367	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	368	ENDDO
	369	ENDDO
	370
	371	DO ig=1,ngrid
	372	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	373	& zb(ig,2)*(1.-zd(ig,2)))
	374	zc(ig,1)=(za(ig,1)*zh(ig,1)+
	375	& zb(ig,2)zc(ig,2))z1(ig)
	376	zd(ig,1)=zb(ig,1)*z1(ig)
	377	ENDDO
	378
	379	! Calculate (d Planck / dT) at the interface temperature
	380	! ------------------------------------------------------
	381
	382	z4st=4.0sigmaptimestep
	383	DO ig=1,ngrid
	384	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
	385	ENDDO
	386
	387	! Calculate temperature tendency at the interface (dry case)
	388	! ----------------------------------------------------------
	389	! Sum of fluxes at interface at time t + \delta t gives change in T:
	390	! radiative fluxes
	391	! turbulent convective (sensible) heat flux
	392	! flux (if any) from subsurface
	393
	394	! if(.not.water) then
	395
	396	! DO ig=1,ngrid
	397
	398	! z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
	399	! & + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
[3195]	400	! z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
[3184]	401	! & +zdplanck(ig)
	402	! ztsrf2(ig) = z1(ig) / z2(ig)
	403	! pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
	404	! zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
	405	! ENDDO
	406
	407	! ! Recalculate temperature to top of atmosphere, starting from ground
	408	! ! ------------------------------------------------------------------
	409
	410	! DO ilay=2,nlay
	411	! DO ig=1,ngrid
	412	! hh = zh(ig,ilay-1)
	413	! zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
	414	! ENDDO
	415	! ENDDO
	416
	417	! endif ! not water
	418
	419	!-----------------------------------------------------------------------
	420	! TRACERS (no vapour)
	421	! -------
	422
	423	if(tracer) then
	424
	425	! Calculate vertical flux from the bottom to the first layer (dust)
	426	! -----------------------------------------------------------------
[3195]	427	do ig=1,ngrid
[3184]	428	rho(ig) = zb0(ig,1) /ptimestep
	429	end do
	430
	431	pdqsdif(1:ngrid,1:nq)=0
	432
	433	! Implicit inversion of q
	434	! -----------------------
[3195]	435	do iq=1,nq
[3184]	436
[3275]	437	! if (iq.ne.igcm_h2o_gas) then
[3184]	438
	439	DO ig=1,ngrid
	440	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	441	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	442	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
[3195]	443	ENDDO
	444
[3184]	445	DO ilay=nlay-1,2,-1
	446	DO ig=1,ngrid
	447	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	448	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
	449	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	450	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
	451	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
	452	ENDDO
	453	ENDDO
	454
	455	DO ig=1,ngrid
	456	z1(ig)=1./(za(ig,1)+
	457	& zb(ig,2)*(1.-zdq(ig,2)))
	458	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	459	& zb(ig,2)*zcq(ig,2)
	460	& +(-pdqsdif(ig,iq))ptimestep)z1(ig)
	461	! tracer flux from surface
	462	! currently pdqsdif always zero here,
	463	! so last line is superfluous
	464	enddo
	465	! endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
	466
	467
	468	! Starting upward calculations for simple tracer mixing (e.g., dust)
	469	do ig=1,ngrid
	470	zq(ig,1,iq)=zcq(ig,1)
	471	end do
	472
	473	do ilay=2,nlay
	474	do ig=1,ngrid
	475	zq(ig,ilay,iq)=zcq(ig,ilay)+
	476	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
	477	end do
	478	end do
	479
	480	! Removed (AF24): Calculate temperature tendency including latent heat term
	481	! and assuming an infinite source of water on the ground
	482	! ------------------------------------------------------------------
	483
	484	end do ! of do iq=1,nq
	485	endif ! traceur
	486
	487
	488	!-----------------------------------------------------------------------
	489	! 8. Final calculation of the vertical diffusion tendencies
	490	! -----------------------------------------------------------------
	491
	492	do ilev = 1, nlay
	493	do ig=1,ngrid
	494	pdudif(ig,ilev)=(zu(ig,ilev)-
	495	& (pu(ig,ilev)))/ptimestep
	496	pdvdif(ig,ilev)=(zv(ig,ilev)-
	497	& (pv(ig,ilev)))/ptimestep
[3195]	498	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
[3184]	499
	500	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
	501	enddo
	502	enddo
	503
	504	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
	505	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
[3195]	506	ENDDO
[3184]	507
	508	if (tracer) then
	509	do iq = 1, nq
	510	do ilev = 1, nlay
	511	do ig=1,ngrid
	512	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
	513	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
	514	& ptimestep
	515	enddo
	516	enddo
	517	enddo
	518
	519	! if(water.and.forceWC)then ! force water conservation in model
	520	! ! we calculate the difference and add it to the ground
	521	! ! this is ugly and should be improved in the future
	522	! do ig=1,ngrid
	523	! Wtot=0.0
	524	! do ilay=1,nlay
	525	! masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
	526	! ! Wtot=Wtot+masse*(zq(ig,ilay,iice)-
	527	! ! & (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep))
	528	! Wtot=Wtot+masse*(zq(ig,ilay,ivap)-
	529	! & (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep))
	530	! enddo
	531	! Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
	532
	533	! if(ztsrf2(ig).gt.T_h2O_ice_liq)then
	534	! pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
	535	! else
	536	! pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff
	537	! endif
	538	! enddo
	539
	540	! endif
	541
[3195]	542	endif
[3184]	543
	544	! if(water)then
	545	! call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
	546	! endif
	547
	548
	549	end subroutine vdifc
	550
	551	end module vdifc_mod

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

Download in other formats:

Original Format