Context Navigation

vdifc_mod.F @ 2987

Last change on this file since 2987 was 2427, checked in by emillour, 4 years ago
Generic GCM: Bug fix on call arguments sent from physiq to vdifc (probably not as bad as it sounds, as turbdiff is usually used instead of vdifc). In the process turned vdifc into a module, as well as turbdiff, for better control, and removed unused arguments. EM+YJ
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
[1308]	248	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
[253]	249	& ,pu,pv,ph,zcdv_true
	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