Context Navigation

turbdiff.F90 @ 2276

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

29/08/2018 == JL

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

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

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

Download in other formats:

Original Format