Context Navigation

vdifc.F @ 1448

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

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

Download in other formats:

Original Format