Context Navigation

vdifc.F @ 1542

Last change on this file since 1542 was 1542, checked in by emillour, 9 years ago

Generic GCM:

fix for 1D in writediagfi to enable writing at "ecritphy" rate.
move iniprint.h to "misc"
Some code cleanup in anticipation of future updates:
- changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
- removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

File size: 24.1 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,
[1477]	5	& 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 tracer_h
[1384]	13	use comcstfi_mod, only: g, r, cpp, rcp
[1397]	14	use callkeys_mod, only: water,tracer,nosurf
[135]	15
	16	implicit none
	17
[253]	18	!==================================================================
	19	!
	20	! Purpose
	21	! -------
	22	! Turbulent diffusion (mixing) for pot. T, U, V and tracers
	23	!
	24	! Implicit scheme
	25	! We start by adding to variables x the physical tendencies
	26	! already computed. We resolve the equation:
	27	!
	28	! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)
	29	!
	30	! Authors
	31	! -------
	32	! F. Hourdin, F. Forget, R. Fournier (199X)
	33	! R. Wordsworth, B. Charnay (2010)
	34	!
	35	!==================================================================
[135]	36
[253]	37	!-----------------------------------------------------------------------
	38	! declarations
	39	! ------------
[135]	40
	41
[253]	42	! arguments
	43	! ---------
[135]	44	INTEGER ngrid,nlay
	45	REAL ptimestep
	46	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	47	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	48	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	49	REAL ptsrf(ngrid),pemis(ngrid)
[1477]	50	REAL pdhfi(ngrid,nlay)
[135]	51	REAL pfluxsrf(ngrid)
	52	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
[594]	53	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
[135]	54	REAL pq2(ngrid,nlay+1)
[253]	55
[1297]	56	real rnat(ngrid)
[135]	57
[253]	58	! Arguments added for condensation
[135]	59	REAL ppopsk(ngrid,nlay)
	60	logical lecrit
	61	REAL pz0
	62
[253]	63	! Tracers
	64	! --------
[135]	65	integer nq
[253]	66	real pqsurf(ngrid,nq)
[135]	67	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	68	real pdqdif(ngrid,nlay,nq)
	69	real pdqsdif(ngrid,nq)
	70
[253]	71	! local
	72	! -----
	73	integer ilev,ig,ilay,nlev
[135]	74
[787]	75	REAL z4st,zdplanck(ngrid)
[1308]	76	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
[787]	77	REAL zcdv(ngrid),zcdh(ngrid)
	78	REAL zcdv_true(ngrid),zcdh_true(ngrid)
[1308]	79	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	80	REAL zh(ngrid,nlay)
[787]	81	REAL ztsrf2(ngrid)
	82	REAL z1(ngrid),z2(ngrid)
[1308]	83	REAL za(ngrid,nlay),zb(ngrid,nlay)
	84	REAL zb0(ngrid,nlay)
	85	REAL zc(ngrid,nlay),zd(ngrid,nlay)
[135]	86	REAL zcst1
[253]	87	REAL zu2!, a
[1308]	88	REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
[787]	89	REAL evap(ngrid)
	90	REAL zcq0(ngrid),zdq0(ngrid)
	91	REAL zx_alf1(ngrid),zx_alf2(ngrid)
[135]	92
	93	LOGICAL firstcall
	94	SAVE firstcall
[1315]	95	!$OMP THREADPRIVATE(firstcall)
[303]	96
	97	LOGICAL lastcall
[135]	98
[253]	99	! variables added for CO2 condensation
	100	! ------------------------------------
[1308]	101	REAL hh !, zhcond(ngrid,nlay)
[253]	102	! REAL latcond,tcond1mb
	103	! REAL acond,bcond
	104	! SAVE acond,bcond
[1315]	105	!!$OMP THREADPRIVATE(acond,bcond)
[253]	106	! DATA latcond,tcond1mb/5.9e5,136.27/
[135]	107
[253]	108	! Tracers
	109	! -------
[135]	110	INTEGER iq
[1308]	111	REAL zq(ngrid,nlay,nq)
[787]	112	REAL zq1temp(ngrid)
	113	REAL rho(ngrid) ! near-surface air density
	114	REAL qsat(ngrid)
[135]	115	DATA firstcall/.true./
	116	REAL kmixmin
	117
[253]	118	! Variables added for implicit latent heat inclusion
	119	! --------------------------------------------------
[787]	120	real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
	121	real z1_Tdry(ngrid), z2_Tdry(ngrid)
	122	real z1_T(ngrid), z2_T(ngrid)
	123	real zb_T(ngrid)
[1308]	124	real zc_T(ngrid,nlay)
	125	real zd_T(ngrid,nlay)
[787]	126	real lat1(ngrid), lat2(ngrid)
[253]	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
[1315]	133	!$OMP THREADPRIVATE(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
[1477]	206	zu(ig,ilev)=pu(ig,ilev)
	207	zv(ig,ilev)=pv(ig,ilev)
[135]	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
[1308]	251	call vdif_kc(ngrid,nlay,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
[1297]	514	if(nint(rnat(ig)).eq.1)then
[253]	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
[1297]	606	if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
[253]	607
	608	pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
	609	pdqsdif(ig,iice)=0.0
	610
	611
[1297]	612	elseif (nint(rnat(ig)).eq.1) then ! (continent)
[253]	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)-
[1477]	663	& (pu(ig,ilev)))/ptimestep
[253]	664	pdvdif(ig,ilev)=(zv(ig,ilev)-
[1477]	665	& (pv(ig,ilev)))/ptimestep
[253]	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