Context Navigation

vdifc.F @ 1862

Last change on this file since 1862 was 1668, checked in by jvatant, 8 years ago
Another round of cleaning of dust and dummy tracers JVO
File size: 14.8 KB

Rev	Line
[1647]	1	subroutine vdifc(ngrid,nlay,nq,ppopsk,
[253]	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
[600]	9	use radcommon_h, only : sigma
[787]	10	USE surfdat_h
	11	USE tracer_h
[1384]	12	use comcstfi_mod, only: g, r, cpp, rcp
[1647]	13	use callkeys_mod, only: tracer,nosurf
[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
[253]	41	! arguments
	42	! ---------
[135]	43	INTEGER ngrid,nlay
	44	REAL ptimestep
	45	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	46	REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
	47	REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
	48	REAL ptsrf(ngrid),pemis(ngrid)
[1477]	49	REAL pdhfi(ngrid,nlay)
[135]	50	REAL pfluxsrf(ngrid)
	51	REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
[594]	52	REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
[135]	53	REAL pq2(ngrid,nlay+1)
[1647]	54
[135]	55
[253]	56	! Arguments added for condensation
[135]	57	REAL ppopsk(ngrid,nlay)
	58	logical lecrit
	59	REAL pz0
	60
[253]	61	! Tracers
	62	! --------
[135]	63	integer nq
[253]	64	real pqsurf(ngrid,nq)
[135]	65	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
	66	real pdqdif(ngrid,nlay,nq)
	67	real pdqsdif(ngrid,nq)
	68
[253]	69	! local
	70	! -----
	71	integer ilev,ig,ilay,nlev
[135]	72
[787]	73	REAL z4st,zdplanck(ngrid)
[1308]	74	REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
[787]	75	REAL zcdv(ngrid),zcdh(ngrid)
	76	REAL zcdv_true(ngrid),zcdh_true(ngrid)
[1308]	77	REAL zu(ngrid,nlay),zv(ngrid,nlay)
	78	REAL zh(ngrid,nlay)
[787]	79	REAL ztsrf2(ngrid)
	80	REAL z1(ngrid),z2(ngrid)
[1308]	81	REAL za(ngrid,nlay),zb(ngrid,nlay)
	82	REAL zb0(ngrid,nlay)
	83	REAL zc(ngrid,nlay),zd(ngrid,nlay)
[135]	84	REAL zcst1
[253]	85	REAL zu2!, a
[1308]	86	REAL zcq(ngrid,nlay),zdq(ngrid,nlay)
[787]	87	REAL evap(ngrid)
	88	REAL zcq0(ngrid),zdq0(ngrid)
	89	REAL zx_alf1(ngrid),zx_alf2(ngrid)
[135]	90
	91	LOGICAL firstcall
	92	SAVE firstcall
[1315]	93	!$OMP THREADPRIVATE(firstcall)
[303]	94
	95	LOGICAL lastcall
[135]	96
[253]	97	! variables added for CO2 condensation
	98	! ------------------------------------
[1668]	99	REAL hh
[135]	100
[253]	101	! Tracers
	102	! -------
[135]	103	INTEGER iq
[1308]	104	REAL zq(ngrid,nlay,nq)
[787]	105	REAL zq1temp(ngrid)
	106	REAL rho(ngrid) ! near-surface air density
	107	REAL qsat(ngrid)
[135]	108	DATA firstcall/.true./
	109	REAL kmixmin
	110
[253]	111	real, parameter :: karman=0.4
	112	real cd0, roughratio
[135]	113
[253]	114	real masse, Wtot, Wdiff
[135]	115
[253]	116	real dqsdif_total(ngrid)
	117	real zq0(ngrid)
[135]	118
	119
[253]	120	! Coherence test
	121	! --------------
[135]	122
[253]	123	IF (firstcall) THEN
	124	firstcall=.false.
	125	ENDIF
	126
	127	!-----------------------------------------------------------------------
	128	! 1. Initialisation
	129	! -----------------
	130
[135]	131	nlev=nlay+1
	132
[253]	133	! Calculate rhodz and dtrho/dz=dtrho*2 g/dp
	134	! with rho=p/RT=p/ (R Theta) (p/ps)**kappa
	135	! ---------------------------------------------
[135]	136
	137	DO ilay=1,nlay
	138	DO ig=1,ngrid
	139	za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
	140	ENDDO
	141	ENDDO
	142
[253]	143	zcst1=4.gptimestep/(R*R)
[135]	144	DO ilev=2,nlev-1
	145	DO ig=1,ngrid
	146	zb0(ig,ilev)=pplev(ig,ilev)*
[253]	147	s (pplev(ig,1)/pplev(ig,ilev))**rcp /
	148	s (ph(ig,ilev-1)+ph(ig,ilev))
[135]	149	zb0(ig,ilev)=zcst1zb0(ig,ilev)zb0(ig,ilev)/
[253]	150	s (pplay(ig,ilev-1)-pplay(ig,ilev))
[135]	151	ENDDO
	152	ENDDO
	153	DO ig=1,ngrid
[253]	154	zb0(ig,1)=ptimesteppplev(ig,1)/(Rptsrf(ig))
[135]	155	ENDDO
	156
[253]	157	dqsdif_total(:)=0.0
[135]	158
[253]	159	!-----------------------------------------------------------------------
	160	! 2. Add the physical tendencies computed so far
	161	! ----------------------------------------------
[135]	162
	163	DO ilev=1,nlay
	164	DO ig=1,ngrid
[1477]	165	zu(ig,ilev)=pu(ig,ilev)
	166	zv(ig,ilev)=pv(ig,ilev)
[135]	167	zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	168	ENDDO
	169	ENDDO
	170	if(tracer) then
[253]	171	DO iq =1, nq
	172	DO ilev=1,nlay
	173	DO ig=1,ngrid
	174	zq(ig,ilev,iq)=pq(ig,ilev,iq) +
	175	& pdqfi(ig,ilev,iq)*ptimestep
	176	ENDDO
	177	ENDDO
[135]	178	ENDDO
	179	end if
	180
[253]	181	!-----------------------------------------------------------------------
	182	! 3. Turbulence scheme
	183	! --------------------
	184	!
	185	! Source of turbulent kinetic energy at the surface
	186	! -------------------------------------------------
	187	! Formula is Cd_0 = (karman / log[1+z1/z0])^2
[135]	188
[253]	189	DO ig=1,ngrid
	190	roughratio = 1.E+0 + pzlay(ig,1)/pz0
	191	cd0 = karman/log(roughratio)
	192	cd0 = cd0*cd0
	193	zcdv_true(ig) = cd0
	194	zcdh_true(ig) = cd0
[952]	195	if (nosurf) then
	196	zcdv_true(ig) = 0. !! disable sensible momentum flux
	197	zcdh_true(ig) = 0. !! disable sensible heat flux
	198	endif
[253]	199	ENDDO
[135]	200
	201	DO ig=1,ngrid
[253]	202	zu2=pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1)
	203	zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
	204	zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
[135]	205	ENDDO
	206
[253]	207	! Turbulent diffusion coefficients in the boundary layer
	208	! ------------------------------------------------------
[135]	209
[1308]	210	call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
[253]	211	& ,pu,pv,ph,zcdv_true
	212	& ,pq2,zkv,zkh)
[135]	213
[253]	214	! Adding eddy mixing to mimic 3D general circulation in 1D
	215	! R. Wordsworth & F. Forget (2010)
[135]	216	if ((ngrid.eq.1)) then
[253]	217	kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D
	218	do ilev=1,nlay
	219	do ig=1,ngrid
	220	!zkh(ig,ilev) = 1.0
	221	zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
	222	zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
	223	end do
	224	end do
[135]	225	end if
	226
[253]	227	!-----------------------------------------------------------------------
	228	! 4. Implicit inversion of u
	229	! --------------------------
[135]	230
[253]	231	! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1)
	232	! avec
	233	! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.)
	234	! et
	235	! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
	236	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	237	! et /zkv/ = Ku
	238
[135]	239	CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
	240	CALL multipl(ngrid,zcdv,zb0,zb)
	241
	242	DO ig=1,ngrid
	243	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	244	zc(ig,nlay)=za(ig,nlay)zu(ig,nlay)z1(ig)
	245	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	246	ENDDO
	247
	248	DO ilay=nlay-1,1,-1
	249	DO ig=1,ngrid
	250	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	251	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	252	zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
[253]	253	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	254	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	255	ENDDO
	256	ENDDO
	257
	258	DO ig=1,ngrid
	259	zu(ig,1)=zc(ig,1)
	260	ENDDO
	261	DO ilay=2,nlay
	262	DO ig=1,ngrid
	263	zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
	264	ENDDO
	265	ENDDO
	266
[253]	267	!-----------------------------------------------------------------------
	268	! 5. Implicit inversion of v
	269	! --------------------------
[135]	270
[253]	271	! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1)
	272	! avec
	273	! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.)
	274	! et
	275	! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
	276	! donc les entrees sont /zcdv/ pour la condition a la limite sol
	277	! et /zkv/ = Kv
[135]	278
	279	DO ig=1,ngrid
	280	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	281	zc(ig,nlay)=za(ig,nlay)zv(ig,nlay)z1(ig)
	282	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	283	ENDDO
	284
	285	DO ilay=nlay-1,1,-1
	286	DO ig=1,ngrid
	287	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	288	$ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	289	zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
[253]	290	$ zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	291	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	292	ENDDO
	293	ENDDO
	294
	295	DO ig=1,ngrid
	296	zv(ig,1)=zc(ig,1)
	297	ENDDO
	298	DO ilay=2,nlay
	299	DO ig=1,ngrid
	300	zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
	301	ENDDO
	302	ENDDO
	303
[253]	304	!----------------------------------------------------------------------------
	305	! 6. Implicit inversion of h, not forgetting the coupling with the ground
[135]	306
[253]	307	! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1)
	308	! avec
	309	! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.)
	310	! et
	311	! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
	312	! donc les entrees sont /zcdh/ pour la condition de raccord au sol
	313	! et /zkh/ = Kh
[135]	314
[253]	315	! Using the wind modified by friction for lifting and sublimation
	316	! ---------------------------------------------------------------
	317	DO ig=1,ngrid
	318	zu2 = zu(ig,1)zu(ig,1)+zv(ig,1)zv(ig,1)
	319	zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
	320	zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
	321	ENDDO
	322
[135]	323	CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
	324	CALL multipl(ngrid,zcdh,zb0,zb)
	325
	326	DO ig=1,ngrid
	327	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	328	zc(ig,nlay)=za(ig,nlay)zh(ig,nlay)z1(ig)
	329	zd(ig,nlay)=zb(ig,nlay)*z1(ig)
	330	ENDDO
	331
[253]	332	DO ilay=nlay-1,2,-1
[135]	333	DO ig=1,ngrid
	334	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
[253]	335	& zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[135]	336	zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
[253]	337	& zb(ig,ilay+1)zc(ig,ilay+1))z1(ig)
[135]	338	zd(ig,ilay)=zb(ig,ilay)*z1(ig)
	339	ENDDO
	340	ENDDO
	341
	342	DO ig=1,ngrid
[253]	343	z1(ig)=1./(za(ig,1)+zb(ig,1)+
	344	& zb(ig,2)*(1.-zd(ig,2)))
	345	zc(ig,1)=(za(ig,1)*zh(ig,1)+
	346	& zb(ig,2)zc(ig,2))z1(ig)
	347	zd(ig,1)=zb(ig,1)*z1(ig)
[135]	348	ENDDO
	349
[253]	350	! Calculate (d Planck / dT) at the interface temperature
	351	! ------------------------------------------------------
[135]	352
[253]	353	z4st=4.0sigmaptimestep
[135]	354	DO ig=1,ngrid
[253]	355	zdplanck(ig)=z4stpemis(ig)ptsrf(ig)ptsrf(ig)ptsrf(ig)
[135]	356	ENDDO
	357
[253]	358	! Calculate temperature tendency at the interface (dry case)
	359	! ----------------------------------------------------------
	360	! Sum of fluxes at interface at time t + \delta t gives change in T:
	361	! radiative fluxes
	362	! turbulent convective (sensible) heat flux
	363	! flux (if any) from subsurface
[135]	364
[253]	365
[135]	366	DO ig=1,ngrid
[253]	367
	368	z1(ig) = pcapcal(ig)ptsrf(ig) + cppzb(ig,1)*zc(ig,1)
	369	& + zdplanck(ig)ptsrf(ig) + pfluxsrf(ig)ptimestep
	370	z2(ig) = pcapcal(ig) + cppzb(ig,1)(1.-zd(ig,1))
	371	& +zdplanck(ig)
	372	ztsrf2(ig) = z1(ig) / z2(ig)
	373	pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
	374	zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
[135]	375	ENDDO
	376
[253]	377	! Recalculate temperature to top of atmosphere, starting from ground
	378	! ------------------------------------------------------------------
[135]	379
[253]	380	DO ilay=2,nlay
	381	DO ig=1,ngrid
	382	hh = zh(ig,ilay-1)
	383	zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
	384	ENDDO
	385	ENDDO
[135]	386
	387
[253]	388	!-----------------------------------------------------------------------
	389	! TRACERS (no vapour)
	390	! -------
[135]	391
[253]	392	if(tracer) then
[135]	393
[253]	394	! Calculate vertical flux from the bottom to the first layer (dust)
	395	! -----------------------------------------------------------------
[787]	396	do ig=1,ngrid
[253]	397	rho(ig) = zb0(ig,1) /ptimestep
	398	end do
[135]	399
[253]	400	call zerophys(ngrid*nq,pdqsdif)
[135]	401
[253]	402	! Implicit inversion of q
	403	! -----------------------
	404	do iq=1,nq
[135]	405
	406	DO ig=1,ngrid
[253]	407	z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
	408	zcq(ig,nlay)=za(ig,nlay)zq(ig,nlay,iq)z1(ig)
	409	zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
	410	ENDDO
	411
	412	DO ilay=nlay-1,2,-1
	413	DO ig=1,ngrid
	414	z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
	415	& zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
	416	zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
	417	& zb(ig,ilay+1)zcq(ig,ilay+1))z1(ig)
	418	zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
	419	ENDDO
[135]	420	ENDDO
	421
[1647]	422
	423
[253]	424	DO ig=1,ngrid
	425	z1(ig)=1./(za(ig,1)+
	426	& zb(ig,2)*(1.-zdq(ig,2)))
	427	zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
	428	& zb(ig,2)*zcq(ig,2)
	429	& +(-pdqsdif(ig,iq))ptimestep)z1(ig)
	430	! tracer flux from surface
	431	! currently pdqsdif always zero here,
	432	! so last line is superfluous
	433	enddo
[135]	434
	435
[253]	436	! Starting upward calculations for simple tracer mixing (e.g., dust)
	437	do ig=1,ngrid
	438	zq(ig,1,iq)=zcq(ig,1)
	439	end do
[135]	440
[253]	441	do ilay=2,nlay
	442	do ig=1,ngrid
	443	zq(ig,ilay,iq)=zcq(ig,ilay)+
	444	$ zdq(ig,ilay)*zq(ig,ilay-1,iq)
	445	end do
	446	end do
[135]	447
	448
[253]	449	end do ! of do iq=1,nq
	450	endif ! traceur
	451
	452
	453	!-----------------------------------------------------------------------
	454	! 8. Final calculation of the vertical diffusion tendencies
	455	! -----------------------------------------------------------------
	456
	457	do ilev = 1, nlay
	458	do ig=1,ngrid
	459	pdudif(ig,ilev)=(zu(ig,ilev)-
[1477]	460	& (pu(ig,ilev)))/ptimestep
[253]	461	pdvdif(ig,ilev)=(zv(ig,ilev)-
[1477]	462	& (pv(ig,ilev)))/ptimestep
[253]	463	hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
	464
[135]	465	pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
[253]	466	enddo
	467	enddo
[594]	468
	469	DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
	470	sensibFlux(ig)=cppzb(ig,1)/ptimestep(zh(ig,1)-ztsrf2(ig))
	471	ENDDO
	472
[253]	473	if (tracer) then
	474	do iq = 1, nq
	475	do ilev = 1, nlay
	476	do ig=1,ngrid
	477	pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
	478	& (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/
	479	& ptimestep
	480	enddo
	481	enddo
	482	enddo
[135]	483
[253]	484	endif
[135]	485
[303]	486	! if(lastcall)then
	487	! if(ngrid.eq.1)then
	488	! print*,'Saving k.out...'
	489	! OPEN(12,file='k.out',form='formatted')
	490	! DO ilay=1,nlay
	491	! write(12,*) zkh(1,ilay), pplay(1,ilay)
	492	! ENDDO
	493	! CLOSE(12)
	494	! endif
	495	! endif
	496
	497
[253]	498	return
	499	end

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

Download in other formats:

Original Format