Context Navigation

vdifc.F @ 1243

Last change on this file since 1243 was 952, checked in by aslmd, 12 years ago
LMDZ.GENERIC. added the possibility to remove surface (nosurf option) and to add an internal heat flux (intheat value, default is 0.)
File size: 24.2 KB

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

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

Download in other formats:

Original Format