Context Navigation

turbdiff_mod.F90 @ 3469

Last change on this file since 3469 was 3236, checked in by gmilcareck, 10 months ago
Add virtual correction for convective adjustment and turbulent diffusion (turbdiff and vdifc) + correction of allocated tables in physiq_mod for varspec.
File size: 28.4 KB

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

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

Download in other formats:

Original Format