Context Navigation

turbdiff.F90 @ 601

Last change on this file since 601 was 600, checked in by jleconte, 13 years ago

Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

OLR now equal to LW net heating/cooling at equilibrium!
As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

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

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

Download in other formats:

Original Format