Context Navigation

vdifc.F @ 293

Last change on this file since 293 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

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

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

Download in other formats:

Original Format