Context Navigation

vdifc_mod.F @ 3532

Last change on this file since 3532 was 3236, checked in by gmilcareck, 11 months ago
Add virtual correction for convective adjustment and turbulent diffusion (turbdiff and vdifc) + correction of allocated tables in physiq_mod for varspec.
File size: 23.8 KB

Rev	Line
[2427]	1	module vdifc_mod
	2
	3	implicit none
	4
	5	contains
	6
[253]	7	subroutine vdifc(ngrid,nlay,nq,rnat,ppopsk,
	8	& ptimestep,pcapcal,lecrit,
	9	& pplay,pplev,pzlay,pzlev,pz0,
	10	& pu,pv,ph,pq,ptsrf,pemis,pqsurf,
[1477]	11	& pdhfi,pdqfi,pfluxsrf,
[594]	12	& pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
[2427]	13	& pdqdif,pdqsdif)
[135]	14
[1993]	15	use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
	16	& ,Psat_water, Lcpdqsat_water
[600]	17	use radcommon_h, only : sigma
[2427]	18	use surfdat_h, only: dryness
	19	use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
[1384]	20	use comcstfi_mod, only: g, r, cpp, rcp
[1397]	21	use callkeys_mod, only: water,tracer,nosurf
[135]	22
	23	implicit none
	24
[253]	25	!==================================================================
	26	!
	27	! Purpose
	28	! -------
	29	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
	30	!
	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)
	36	!
	37	! Authors
	38	! -------
	39	! F. Hourdin, F. Forget, R. Fournier (199X)
	40	! R. Wordsworth, B. Charnay (2010)
	41	!
	42	!==================================================================
[135]	43
[253]	44	!-----------------------------------------------------------------------
	45	! declarations
	46	! ------------
[135]	47
	48
[253]	49	! arguments
	50	! ---------
[2427]	51	INTEGER,INTENT(IN) :: ngrid,nlay
	52	REAL,INTENT(IN) :: ptimestep
	53	REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	54	REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	55	REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	56	REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
	57	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
	58	REAL,INTENT(IN) :: pfluxsrf(ngrid)
	59	REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
	60	REAL,INTENT(OUT) :: pdhdif(ngrid,nlay)
	61	REAL,INTENT(OUT) :: pdtsrf(ngrid),sensibFlux(ngrid)
	62	REAL,INTENT(IN) :: pcapcal(ngrid)
	63	REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
[253]	64
[2427]	65	real,intent(in) :: rnat(ngrid)
[135]	66
[253]	67	! Arguments added for condensation
[2427]	68	REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
	69	logical,intent(in) :: lecrit
	70	REAL,INTENT(IN) :: pz0
[135]	71
[253]	72	! Tracers
	73	! --------
[2427]	74	integer,intent(in) :: nq
	75	real,intent(in) :: pqsurf(ngrid,nq)
	76	real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	77	real,intent(out) :: pdqdif(ngrid,nlay,nq)
	78	real,intent(out) :: pdqsdif(ngrid,nq)
[135]	79
[253]	80	! local
	81	! -----
	82	integer ilev,ig,ilay,nlev
[135]	83
[787]	84	REAL z4st,zdplanck(ngrid)
[1308]	85	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
[787]	86	REAL zcdv(ngrid),zcdh(ngrid)
	87	REAL zcdv_true(ngrid),zcdh_true(ngrid)
[1308]	88	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	89	REAL zh(ngrid,nlay)
[787]	90	REAL ztsrf2(ngrid)
	91	REAL z1(ngrid),z2(ngrid)
[1308]	92	REAL za(ngrid,nlay),zb(ngrid,nlay)
	93	REAL zb0(ngrid,nlay)
	94	REAL zc(ngrid,nlay),zd(ngrid,nlay)
[135]	95	REAL zcst1
[253]	96	REAL zu2!, a
[1308]	97	REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
[787]	98	REAL evap(ngrid)
	99	REAL zcq0(ngrid),zdq0(ngrid)
	100	REAL zx_alf1(ngrid),zx_alf2(ngrid)
[135]	101
	102	LOGICAL firstcall
	103	SAVE firstcall
[1315]	104	!$OMP THREADPRIVATE(firstcall)
[303]	105
[253]	106	! variables added for CO2 condensation
	107	! ------------------------------------
[1308]	108	REAL hh !, zhcond(ngrid,nlay)
[253]	109	! REAL latcond,tcond1mb
	110	! REAL acond,bcond
	111	! SAVE acond,bcond
[1315]	112	!!$OMP THREADPRIVATE(acond,bcond)
[253]	113	! DATA latcond,tcond1mb/5.9e5,136.27/
[135]	114
[253]	115	! Tracers
	116	! -------
[135]	117	INTEGER iq
[1308]	118	REAL zq(ngrid,nlay,nq)
[787]	119	REAL zq1temp(ngrid)
	120	REAL rho(ngrid) ! near-surface air density
[135]	121	DATA firstcall/.true./
	122	REAL kmixmin
	123
[253]	124	! Variables added for implicit latent heat inclusion
	125	! --------------------------------------------------
[1993]	126	real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
[135]	127
[253]	128	integer ivap, iice ! also make liq for clarity on surface...
	129	save ivap, iice
[1315]	130	!$OMP THREADPRIVATE(ivap,iice)
[135]	131
[253]	132	real, parameter :: karman=0.4
	133	real cd0, roughratio
[135]	134
[253]	135	logical forceWC
	136	real masse, Wtot, Wdiff
[135]	137
[253]	138	real dqsdif_total(ngrid)
	139	real zq0(ngrid)
[135]	140
[253]	141	forceWC=.true.
	142	! forceWC=.false.
[135]	143
	144
[253]	145	! Coherence test
	146	! --------------
[135]	147
[253]	148	IF (firstcall) THEN
	149	! To compute: Tcond= 1./(bcond-acondlog(.0095p)) (p in pascal)
	150	! bcond=1./tcond1mb
	151	! acond=r/latcond
	152	! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
	153	! PRINT*,' acond,bcond',acond,bcond
	154
	155	if(water)then
	156	! iliq=igcm_h2o_vap
	157	ivap=igcm_h2o_vap
	158	iice=igcm_h2o_ice ! simply to make the code legible
	159	! to be generalised later
	160	endif
	161
	162	firstcall=.false.
	163	ENDIF
	164
	165	!-----------------------------------------------------------------------
	166	! 1. Initialisation
	167	! -----------------
	168
[135]	169	nlev=nlay+1
	170
[253]	171	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
	172	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
	173	! ---------------------------------------------
[135]	174
	175	DO ilay=1,nlay
	176	DO ig=1,ngrid
	177	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
	178	ENDDO
	179	ENDDO
	180
[253]	181	zcst1=4.gptimestep/(R*R)
[135]	182	DO ilev=2,nlev-1
	183	DO ig=1,ngrid
	184	zb0(ig,ilev)=pplev(ig,ilev)*
[253]	185	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
	186	s (ph(ig,ilev-1)+ph(ig,ilev))
[135]	187	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
[253]	188	s (pplay(ig,ilev-1)-pplay(ig,ilev))
[135]	189	ENDDO
	190	ENDDO
	191	DO ig=1,ngrid
[253]	192	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
[135]	193	ENDDO
	194
[253]	195	dqsdif_total(:)=0.0
[135]	196
[253]	197	!-----------------------------------------------------------------------
	198	! 2. Add the physical tendencies computed so far
	199	! ----------------------------------------------
[135]	200
	201	DO ilev=1,nlay
	202	DO ig=1,ngrid
[1477]	203	zu(ig,ilev)=pu(ig,ilev)
	204	zv(ig,ilev)=pv(ig,ilev)
[135]	205	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	206	ENDDO
	207	ENDDO
	208	if(tracer) then
[253]	209	DO iq =1, nq
	210	DO ilev=1,nlay
	211	DO ig=1,ngrid
	212	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
	213	& pdqfi(ig,ilev,iq)*ptimestep
	214	ENDDO
	215	ENDDO
[135]	216	ENDDO
	217	end if
	218
[253]	219	!-----------------------------------------------------------------------
	220	! 3. Turbulence scheme
	221	! --------------------
	222	!
	223	! Source of turbulent kinetic energy at the surface
	224	! -------------------------------------------------
	225	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
[135]	226
[253]	227	DO ig=1,ngrid
	228	roughratio = 1.E+0 + pzlay(ig,1)/pz0
	229	cd0 = karman/log(roughratio)
	230	cd0 = cd0*cd0
	231	zcdv_true(ig) = cd0
	232	zcdh_true(ig) = cd0
[952]	233	if (nosurf) then
	234	zcdv_true(ig) = 0. !! disable sensible momentum flux
	235	zcdh_true(ig) = 0. !! disable sensible heat flux
	236	endif
[253]	237	ENDDO
[135]	238
	239	DO ig=1,ngrid
[253]	240	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	241	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	242	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
[135]	243	ENDDO
	244
[253]	245	! Turbulent diffusion coefficients in the boundary layer
	246	! ------------------------------------------------------
[135]	247
[3236]	248	call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
	249	& ,pu,pv,pq,ph,zcdv_true
[253]	250	& ,pq2,zkv,zkh)
[135]	251
[253]	252	! Adding eddy mixing to mimic 3D general circulation in 1D
	253	! R. Wordsworth & F. Forget (2010)
[135]	254	if ((ngrid.eq.1)) then
[253]	255	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
	256	do ilev=1,nlay
	257	do ig=1,ngrid
	258	!zkh(ig,ilev) = 1.0
	259	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
	260	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
	261	end do
	262	end do
[135]	263	end if
	264
[253]	265	!-----------------------------------------------------------------------
	266	! 4. Implicit inversion of u
	267	! --------------------------
[135]	268
[253]	269	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	270	! avec
	271	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	272	! et
	273	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	274	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	275	! et /zkv/ = Ku
	276
[135]	277	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	278	CALL multipl(ngrid,zcdv,zb0,zb)
	279
	280	DO ig=1,ngrid
	281	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	282	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
	283	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	284	ENDDO
	285
	286	DO ilay=nlay-1,1,-1
	287	DO ig=1,ngrid
	288	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	289	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	290	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
[253]	291	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	292	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	293	ENDDO
	294	ENDDO
	295
	296	DO ig=1,ngrid
	297	zu(ig,1)=zc(ig,1)
	298	ENDDO
	299	DO ilay=2,nlay
	300	DO ig=1,ngrid
	301	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	302	ENDDO
	303	ENDDO
	304
[253]	305	!-----------------------------------------------------------------------
	306	! 5. Implicit inversion of v
	307	! --------------------------
[135]	308
[253]	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
[135]	316
	317	DO ig=1,ngrid
	318	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	319	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	320	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	321	ENDDO
	322
	323	DO ilay=nlay-1,1,-1
	324	DO ig=1,ngrid
	325	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	326	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	327	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
[253]	328	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	329	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	330	ENDDO
	331	ENDDO
	332
	333	DO ig=1,ngrid
	334	zv(ig,1)=zc(ig,1)
	335	ENDDO
	336	DO ilay=2,nlay
	337	DO ig=1,ngrid
	338	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	339	ENDDO
	340	ENDDO
	341
[253]	342	!----------------------------------------------------------------------------
	343	! 6. Implicit inversion of h, not forgetting the coupling with the ground
[135]	344
[253]	345	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	346	! avec
	347	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	348	! et
	349	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	350	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
	351	! et /zkh/ = Kh
[135]	352
[253]	353	! Using the wind modified by friction for lifting and sublimation
	354	! ---------------------------------------------------------------
	355	DO ig=1,ngrid
	356	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	357	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
	358	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
	359	ENDDO
	360
[135]	361	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	362	CALL multipl(ngrid,zcdh,zb0,zb)
	363
	364	DO ig=1,ngrid
	365	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	366	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
	367	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	368	ENDDO
	369
[253]	370	DO ilay=nlay-1,2,-1
[135]	371	DO ig=1,ngrid
	372	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	373	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	374	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
[253]	375	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	376	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	377	ENDDO
	378	ENDDO
	379
	380	DO ig=1,ngrid
[253]	381	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	382	& zb(ig,2)*(1.-zd(ig,2)))
	383	zc(ig,1)=(za(ig,1)*zh(ig,1)+
	384	& zb(ig,2)zc(ig,2))z1(ig)
	385	zd(ig,1)=zb(ig,1)*z1(ig)
[135]	386	ENDDO
	387
[253]	388	! Calculate (d Planck / dT) at the interface temperature
	389	! ------------------------------------------------------
[135]	390
[253]	391	z4st=4.0sigmaptimestep
[135]	392	DO ig=1,ngrid
[253]	393	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
[135]	394	ENDDO
	395
[253]	396	! Calculate temperature tendency at the interface (dry case)
	397	! ----------------------------------------------------------
	398	! Sum of fluxes at interface at time t + \delta t gives change in T:
	399	! radiative fluxes
	400	! turbulent convective (sensible) heat flux
	401	! flux (if any) from subsurface
[135]	402
[253]	403	if(.not.water) then
	404
[135]	405	DO ig=1,ngrid
[253]	406
	407	z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
	408	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
	409	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
	410	& +zdplanck(ig)
	411	ztsrf2(ig) = z1(ig) / z2(ig)
	412	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
	413	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
[135]	414	ENDDO
	415
[253]	416	! Recalculate temperature to top of atmosphere, starting from ground
	417	! ------------------------------------------------------------------
[135]	418
[253]	419	DO ilay=2,nlay
	420	DO ig=1,ngrid
	421	hh = zh(ig,ilay-1)
	422	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
	423	ENDDO
	424	ENDDO
[135]	425
[253]	426	endif ! not water
[135]	427
[253]	428	!-----------------------------------------------------------------------
	429	! TRACERS (no vapour)
	430	! -------
[135]	431
[253]	432	if(tracer) then
[135]	433
[253]	434	! Calculate vertical flux from the bottom to the first layer (dust)
	435	! -----------------------------------------------------------------
[787]	436	do ig=1,ngrid
[253]	437	rho(ig) = zb0(ig,1) /ptimestep
	438	end do
[135]	439
[2278]	440	pdqsdif(1:ngrid,1:nq)=0
[135]	441
[253]	442	! Implicit inversion of q
	443	! -----------------------
	444	do iq=1,nq
[135]	445
[253]	446	if (iq.ne.igcm_h2o_vap) then
[135]	447
	448	DO ig=1,ngrid
[253]	449	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	450	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	451	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
	452	ENDDO
	453
	454	DO ilay=nlay-1,2,-1
	455	DO ig=1,ngrid
	456	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	457	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
	458	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	459	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
	460	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
	461	ENDDO
[135]	462	ENDDO
	463
[253]	464	if ((water).and.(iq.eq.iice)) then
	465	! special case for water ice tracer: do not include
	466	! h2o ice tracer from surface (which is set when handling
	467	! h2o vapour case (see further down).
	468	! zb(ig,1)=0 if iq ne ivap
	469	DO ig=1,ngrid
	470	z1(ig)=1./(za(ig,1)+
	471	& zb(ig,2)*(1.-zdq(ig,2)))
	472	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	473	& zb(ig,2)zcq(ig,2))z1(ig)
	474	ENDDO
	475	else ! general case
	476	DO ig=1,ngrid
	477	z1(ig)=1./(za(ig,1)+
	478	& zb(ig,2)*(1.-zdq(ig,2)))
	479	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	480	& zb(ig,2)*zcq(ig,2)
	481	& +(-pdqsdif(ig,iq))ptimestep)z1(ig)
	482	! tracer flux from surface
	483	! currently pdqsdif always zero here,
	484	! so last line is superfluous
	485	enddo
	486	endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
[135]	487
	488
[253]	489	! Starting upward calculations for simple tracer mixing (e.g., dust)
	490	do ig=1,ngrid
	491	zq(ig,1,iq)=zcq(ig,1)
	492	end do
[135]	493
[253]	494	do ilay=2,nlay
	495	do ig=1,ngrid
	496	zq(ig,ilay,iq)=zcq(ig,ilay)+
	497	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
	498	end do
	499	end do
[135]	500
[253]	501	endif ! if (iq.ne.igcm_h2o_vap)
[135]	502
[253]	503	! Calculate temperature tendency including latent heat term
	504	! and assuming an infinite source of water on the ground
	505	! ------------------------------------------------------------------
[135]	506
[253]	507	if (water.and.(iq.eq.igcm_h2o_vap)) then
	508
	509	! compute evaporation efficiency
	510	do ig = 1, ngrid
[1297]	511	if(nint(rnat(ig)).eq.1)then
[253]	512	dryness(ig)=pqsurf(ig,ivap)+pqsurf(ig,iice)
	513	dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
	514	dryness(ig)=MAX(0.,dryness(ig))
	515	endif
	516	enddo
[135]	517
[253]	518	do ig=1,ngrid
	519	! Calculate the value of qsat at the surface (water)
[1993]	520	call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
	521	call Lcpdqsat_water(ptsrf(ig),pplev(ig,1),psat(ig),
	522	& qsat(ig),dqsat(ig),psat_temp)
	523	dqsat(ig)=dqsat(ig)*RCPD/RLVTT
[253]	524	enddo
	525
	526	! coefficients for q
	527
	528	do ig=1,ngrid
	529	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	530	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	531	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
	532	enddo
	533
	534	do ilay=nlay-1,2,-1
	535	do ig=1,ngrid
	536	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	537	$ zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
	538	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	539	$ zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
	540	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
	541	enddo
	542	enddo
	543
	544	do ig=1,ngrid
	545	z1(ig)=1./(za(ig,1)+zb(ig,1)*dryness(ig)+
	546	$ zb(ig,2)*(1.-zdq(ig,2)))
	547	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	548	$ zb(ig,2)zcq(ig,2))z1(ig)
	549	zdq(ig,1)=dryness(ig)zb(ig,1)z1(ig)
	550	enddo
	551
	552	! calculation of h0 and h1
	553	do ig=1,ngrid
	554	zdq0(ig) = dqsat(ig)
	555	zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
	556
	557	z1(ig) = pcapcal(ig)ptsrf(ig) +cppzb(ig,1)*zc(ig,1)
	558	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
	559	& + zb(ig,1)dryness(ig)RLVTT*
	560	& ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
	561
	562	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
	563	& +zdplanck(ig)
	564	& +zb(ig,1)dryness(ig)RLVTTzdq0(ig)
	565	& (1.0-zdq(ig,1))
	566
	567	ztsrf2(ig) = z1(ig) / z2(ig)
	568	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
	569	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
	570	enddo
	571
	572	! calculation of qs and q1
	573	do ig=1,ngrid
	574	zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf2(ig)
	575	zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
	576	enddo
	577
	578	! calculation of evaporation
	579	do ig=1,ngrid
	580	evap(ig)= zb(ig,1)dryness(ig)(zq(ig,1,ivap)-zq0(ig))
	581	dqsdif_total(ig)=evap(ig)
	582	enddo
	583
	584	! recalculate temperature and q(vap) to top of atmosphere, starting from ground
	585	do ilay=2,nlay
	586	do ig=1,ngrid
	587	zq(ig,ilay,iq)=zcq(ig,ilay)
	588	& +zdq(ig,ilay)*zq(ig,ilay-1,iq)
	589	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
	590	end do
	591	end do
	592
	593	do ig=1,ngrid
	594
	595	! --------------------------------------------------------------------------
	596	! On the ocean, if T > 0 C then the vapour tendency must replace the ice one
	597	! The surface vapour tracer is actually liquid. To make things difficult.
	598
[1297]	599	if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
[253]	600
	601	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
	602	pdqsdif(ig,iice)=0.0
	603
	604
[1297]	605	elseif (nint(rnat(ig)).eq.1) then ! (continent)
[253]	606
	607	! --------------------------------------------------------
	608	! Now check if we've taken too much water from the surface
	609	! This can only occur on the continent
	610
	611	! If water is evaporating / subliming, we take it from ice before liquid
	612	! -- is this valid??
	613	if(dqsdif_total(ig).lt.0)then
	614	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
	615	pdqsdif(ig,iice)=max(-pqsurf(ig,iice)/ptimestep
	616	& ,pdqsdif(ig,iice))
	617	endif
	618	! sublimation only greater than qsurf(ice)
	619	! ----------------------------------------
	620	! we just convert some liquid to vapour too
	621	! if latent heats are the same, no big deal
	622	if (-dqsdif_total(ig).gt.pqsurf(ig,iice))then
	623	pdqsdif(ig,iice) = -pqsurf(ig,iice)/ptimestep ! removes all the ice!
	624	pdqsdif(ig,ivap) = dqsdif_total(ig)/ptimestep
	625	& - pdqsdif(ig,iice) ! take the remainder from the liquid instead
	626	pdqsdif(ig,ivap) = max(-pqsurf(ig,ivap)/ptimestep
	627	& ,pdqsdif(ig,ivap))
	628	endif
	629
	630	endif ! if (rnat.ne.1)
	631
	632	! If water vapour is condensing, we must decide whether it forms ice or liquid.
	633	if(dqsdif_total(ig).gt.0)then ! a bug was here!
[650]	634	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
[253]	635	pdqsdif(ig,iice)=0.0
	636	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
	637	else
	638	pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
	639	pdqsdif(ig,ivap)=0.0
	640	endif
	641	endif
	642
	643	end do ! of DO ig=1,ngrid
	644	endif ! if (water et iq=ivap)
	645	end do ! of do iq=1,nq
	646	endif ! traceur
	647
	648
	649	!-----------------------------------------------------------------------
	650	! 8. Final calculation of the vertical diffusion tendencies
	651	! -----------------------------------------------------------------
	652
	653	do ilev = 1, nlay
	654	do ig=1,ngrid
	655	pdudif(ig,ilev)=(zu(ig,ilev)-
[1477]	656	& (pu(ig,ilev)))/ptimestep
[253]	657	pdvdif(ig,ilev)=(zv(ig,ilev)-
[1477]	658	& (pv(ig,ilev)))/ptimestep
[253]	659	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	660
[135]	661	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
[253]	662	enddo
	663	enddo
[594]	664
	665	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
	666	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
	667	ENDDO
	668
[253]	669	if (tracer) then
	670	do iq = 1, nq
	671	do ilev = 1, nlay
	672	do ig=1,ngrid
	673	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
	674	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
	675	& ptimestep
	676	enddo
	677	enddo
	678	enddo
[135]	679
[253]	680	if(water.and.forceWC)then ! force water conservation in model
	681	! we calculate the difference and add it to the ground
	682	! this is ugly and should be improved in the future
	683	do ig=1,ngrid
	684	Wtot=0.0
	685	do ilay=1,nlay
	686	masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
	687	! Wtot=Wtot+masse*(zq(ig,ilay,iice)-
	688	! & (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep))
	689	Wtot=Wtot+masse*(zq(ig,ilay,ivap)-
	690	& (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep))
	691	enddo
	692	Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
[135]	693
[650]	694	if(ztsrf2(ig).gt.T_h2O_ice_liq)then
[253]	695	pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
	696	else
	697	pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff
	698	endif
	699	enddo
[135]	700
[253]	701	endif
[135]	702
[253]	703	endif
[135]	704
[253]	705	if(water)then
	706	call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
	707	endif
	708
[303]	709
[2427]	710	end subroutine vdifc
[303]	711
[2427]	712	end module vdifc_mod

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

Download in other formats:

Original Format