Context Navigation

vdifc.F @ 486

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

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

Download in other formats:

Original Format