Context Navigation

vdifc.F @ 2176

Last change on this file since 2176 was 1993, checked in by jleconte, 6 years ago

29/08/2018 == JL

-watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...) which used Psat_water from watercommon.
This is now harmonized. ALl routines use Psat_water. Watersat.F has been removed, but the routine is now in watercommon for archival purpose. It is not used anymore.
-also changed the number of chars for tname in the dyn3D/infotrac.F90 to be able to run rcm1d.

File size: 23.7 KB

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

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

Download in other formats:

Original Format