Context Navigation

vdifc.F @ 724

Last change on this file since 724 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

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

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

Download in other formats:

Original Format