Context Navigation

turbdiff.F90 @ 1202

Last change on this file since 1202 was 1194, checked in by sglmd, 11 years ago
Latitude-dependent gravity field added. Option oblate = .true. in callphys.def, and three additional variables needed: J2, Rmean and MassPlanet?.
File size: 27.0 KB

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

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

Download in other formats:

Original Format