Context Navigation

yamada4.F90

Last change on this file was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 29.3 KB

Rev	Line
[3331]	1	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	2
	3	SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
	4	cd, q2, km, kn, kq, ustar, iflag_pbl)
	5
	6	USE dimphy
	7	USE ioipsl_getin_p_mod, ONLY : getin_p
	8
	9	IMPLICIT NONE
	10	include "iniprint.h"
	11	! .......................................................................
	12	! ym#include "dimensions.h"
	13	! ym#include "dimphy.h"
	14	! ************************************************************************************************
	15	!
	16	! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
	17	!
	18	! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
	19	! Yamada T.
	20	! J. Atmos. Sci, 40, 91-106, 1983
	21	!
	22	!************************************************************************************************
	23	! Input :
	24	!'======
	25	! ni: indice horizontal sur la grille de base, non restreinte
	26	! nsrf: type de surface
	27	! ngrid: nombre de mailles concern??es sur l'horizontal
	28	! dt : pas de temps
	29	! g : g
	30	! rconst: constante de l'air sec
	31	! zlev : altitude a chaque niveau (interface inferieure de la couche
	32	! de meme indice)
	33	! zlay : altitude au centre de chaque couche
	34	! u,v : vitesse au centre de chaque couche
	35	! (en entree : la valeur au debut du pas de temps)
	36	! teta : temperature potentielle virtuelle au centre de chaque couche
	37	! (en entree : la valeur au debut du pas de temps)
	38	! cd : cdrag pour la quantit?? de mouvement
	39	! (en entree : la valeur au debut du pas de temps)
	40	! ustar: vitesse de friction calcul??e par une formule diagnostique
	41	! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
	42	!
	43	! iflag_pbl doit valoir entre 6 et 9
	44	! l=6, on prend systematiquement une longueur d'equilibre
	45	! iflag_pbl=6 : MY 2.0
	46	! iflag_pbl=7 : MY 2.0.Fournier
	47	! iflag_pbl=8/9 : MY 2.5
	48	! iflag_pbl=8 with special obsolete treatments for convergence
	49	! with Cmpi5 NPv3.1 simulations
	50	! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
	51	! iflag_pbl=12 = 11 with vertical diffusion off q2
	52	!
	53	! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
	54	! Correction for very stable PBLs (iflag_pbl=10 and 11)
	55	! iflag_pbl=8 converges numerically with NPv3.1
	56	! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
	57	! -> the model can run with longer time-steps.
	58	!
	59	! Inpout/Output :
	60	!==============
	61	! q2 : $q^2$ au bas de chaque couche
	62	! (en entree : la valeur au debut du pas de temps)
	63	! (en sortie : la valeur a la fin du pas de temps)
	64
	65	! Outputs:
	66	!==========
	67	! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	68	! couche)
	69	! (en sortie : la valeur a la fin du pas de temps)
	70	! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	71	! (en sortie : la valeur a la fin du pas de temps)
	72	!
	73	!.......................................................................
	74
	75	!=======================================================================
	76	! Declarations:
	77	!=======================================================================
	78
	79
	80	! Inputs/Outputs
	81	!----------------
	82	REAL dt, g, rconst
	83	REAL plev(klon, klev+1), temp(klon, klev)
	84	REAL ustar(klon)
	85	REAL kmin, qmin, pblhmin(klon), coriol(klon)
	86	REAL zlev(klon, klev+1)
	87	REAL zlay(klon, klev)
	88	REAL u(klon, klev)
	89	REAL v(klon, klev)
	90	REAL teta(klon, klev)
	91	REAL cd(klon)
	92	REAL q2(klon, klev+1)
	93	REAL unsdz(klon, klev)
	94	REAL unsdzdec(klon, klev+1)
	95	REAL kn(klon, klev+1)
	96	REAL km(klon, klev+1)
	97	INTEGER iflag_pbl, ngrid, nsrf
	98	INTEGER ni(klon)
	99
	100
	101	! Local
	102	!-------
	103
	104	INCLUDE "clesphys.h"
	105
	106
	107	REAL kmpre(klon, klev+1), tmp2, qpre
	108	REAL mpre(klon, klev+1)
	109	REAL kq(klon, klev+1)
	110	REAL ff(klon, klev+1), delta(klon, klev+1)
	111	REAL aa(klon, klev+1), aa0, aa1
	112	INTEGER nlay, nlev
	113	LOGICAL first
	114	INTEGER ipas
	115	SAVE first, ipas
	116	! FH/IM data first,ipas/.true.,0/
	117	DATA first, ipas/.FALSE., 0/
	118	!$OMP THREADPRIVATE( first,ipas)
	119	LOGICAL hboville
	120	INTEGER ig, k
	121	REAL ri, zrif, zalpha, zsm, zsn
	122	REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
	123	REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
	124	REAL dtetadz(klon, klev+1)
	125	REAL m2cstat, mcstat, kmcstat
	126	REAL l(klon, klev+1)
	127	REAL zz(klon, klev+1)
	128	INTEGER iter
	129	REAL ric, rifc, b1, kap
	130	SAVE ric, rifc, b1, kap
	131	DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
	132	!$OMP THREADPRIVATE(ric,rifc,b1,kap)
	133	REAL seuilsm, seuilalpha
	134	REAL,SAVE :: lmixmin
	135	!$OMP THREADPRIVATE(lmixmin)
	136	REAL frif, falpha, fsm
	137	REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
	138	lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
	139	LOGICAL, SAVE :: firstcall = .TRUE.
	140	!$OMP THREADPRIVATE(firstcall)
	141
	142
	143	! Fonctions utilis??es
	144	!--------------------
	145
	146	frif(ri) = 0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	147	falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
	148	fsm(ri) = 1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	149
	150
	151	IF (firstcall) THEN
	152	firstcall = .FALSE.
	153	lmixmin=1.
	154	CALL getin_p('lmixmin',lmixmin)
	155	END IF
	156
	157
	158
	159	!===============================================================================
	160	! Flags, tests et ??valuations de constantes
	161	!===============================================================================
	162
	163	! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
	164	hboville=.TRUE.
	165
	166
	167	IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
	168	STOP 'probleme de coherence dans appel a MY'
	169	END IF
	170
	171	! Seuil dans le code de turbulence
	172	seuilalpha=1.12
	173	seuilsm=0.085
	174
	175	nlay = klev
	176	nlev = klev + 1
	177	ipas = ipas + 1
	178
	179
	180	!========================================================================
	181	! Calcul des increments verticaux
	182	!=========================================================================
	183
	184
	185	! Attention: zlev n'est pas declare a nlev
	186	DO ig = 1, ngrid
	187	zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
	188	END DO
	189
	190
	191	DO k = 1, nlay
	192	DO ig = 1, ngrid
	193	unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	194	END DO
	195	END DO
	196	DO ig = 1, ngrid
	197	unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
	198	END DO
	199	DO k = 2, nlay
	200	DO ig = 1, ngrid
	201	unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	202	END DO
	203	END DO
	204	DO ig = 1, ngrid
	205	unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	206	END DO
	207
	208	!=========================================================================
	209	! Richardson number and stability functions
	210	!=========================================================================
	211
	212	! initialize arrays:
	213
	214	m2(1:ngrid, :) = 0.0
	215	sm(1:ngrid, :) = 0.0
	216	rif(1:ngrid, :) = 0.0
	217
	218	!------------------------------------------------------------
	219	DO k = 2, klev
	220
	221	DO ig = 1, ngrid
	222	dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
	223	m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
	224	k-1))*2)/(dz(ig,k)dz(ig,k))
	225	dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
	226	n2(ig, k) = g2.dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
	227	ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
	228	IF (ri<ric) THEN
	229	rif(ig, k) = frif(ri)
	230	ELSE
	231	rif(ig, k) = rifc
	232	END IF
	233
	234	IF (rif(ig,k)<0.16) THEN
	235	alpha(ig, k) = falpha(rif(ig,k))
	236	sm(ig, k) = fsm(rif(ig,k))
	237	ELSE
	238	alpha(ig, k) = seuilalpha
	239	sm(ig, k) = seuilsm
	240	END IF
	241	zz(ig, k) = b1m2(ig, k)(1.-rif(ig,k))*sm(ig, k)
	242	END DO
	243	END DO
	244
	245
	246
	247	! ====================================================================
	248	! Computing the mixing length
	249	! ====================================================================
	250
	251
	252	CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
	253
	254
	255
	256	!=======================================================================
	257	! DIFFERENT TYPES DE SCHEMA de YAMADA
	258	!=======================================================================
	259
	260	!--------------
	261	! Yamada 2.0
	262	!--------------
	263	IF (iflag_pbl==6) THEN
	264
	265	DO k = 2, klev
	266	q2(1:ngrid, k) = l(1:ngrid, k)*2zz(1:ngrid, k)
	267	END DO
	268
	269
	270	!------------------
	271	! Yamada 2.Fournier
	272	!------------------
	273
	274	ELSE IF (iflag_pbl==7) THEN
	275
	276
	277	! Calcul de l, km, au pas precedent
	278	!....................................
	279	DO k = 2, klev
	280	DO ig = 1, ngrid
	281	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
	282	kmpre(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
	283	mpre(ig, k) = sqrt(m2(ig,k))
	284	END DO
	285	END DO
	286
	287	DO k = 2, klev - 1
	288	DO ig = 1, ngrid
	289	m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
	290	mcstat = sqrt(m2cstat)
	291
	292	! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
	293	!.........................................................................
	294
	295	IF (k==2) THEN
	296	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
	297	unsdz(ig,k-1)cd(ig)(sqrt(u(ig,3)2+v(ig,3)2)-mcstat/unsdzdec &
	298	(ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
	299	-1))
	300	ELSE
	301	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
	302	unsdz(ig,k-1)kmpre(ig,k-1)mpre(ig,k-1))/ &
	303	(unsdz(ig,k)+unsdz(ig,k-1))
	304	END IF
	305
	306	tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
	307	q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
	308
	309	END DO
	310	END DO
	311
	312
	313	! ------------------------
	314	! Yamada 2.5 a la Didi
	315	!-------------------------
	316
	317	ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
	318
	319	! Calcul de l, km, au pas precedent
	320	!....................................
	321	DO k = 2, klev
	322	DO ig = 1, ngrid
	323	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
	324	IF (delta(ig,k)<1.E-20) THEN
	325	delta(ig, k) = 1.E-20
	326	END IF
	327	km(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
	328	aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	329	aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	330	aa(ig, k) = aa1dt/(delta(ig,k)l(ig,k))
	331	qpre = sqrt(q2(ig,k))
	332	IF (aa(ig,k)>0.) THEN
	333	q2(ig, k) = (qpre+aa(ig,k)qpreqpre)**2
	334	ELSE
	335	q2(ig, k) = (qpre/(1.-aa(ig,k)qpre))*2
	336	END IF
	337	! else ! iflag_pbl=9
	338	! if (aa(ig,k)*qpre.gt.0.9) then
	339	! q2(ig,k)=(qpre10.)*2
	340	! else
	341	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	342	! endif
	343	! endif
	344	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
	345	END DO
	346	END DO
	347
	348	ELSE IF (iflag_pbl>=10) THEN
	349
	350	DO k = 2, klev - 1
	351	km(1:ngrid, k) = l(1:ngrid, k)sqrt(q2(1:ngrid,k))sm(1:ngrid, k)
	352	q2(1:ngrid, k) = q2(1:ngrid, k) + dtkm(1:ngrid, k)m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
	353	q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
	354	q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2l(1:ngrid,k)b1))
	355	q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
	356	END DO
	357
	358
	359	ELSE
	360	STOP 'Cas nom prevu dans yamada4'
	361
	362	END IF ! Fin du cas 8
	363
	364
	365	! ====================================================================
	366	! Calcul des coefficients de melange
	367	! ====================================================================
	368
	369	DO k = 2, klev
	370	DO ig = 1, ngrid
	371	zq = sqrt(q2(ig,k))
	372	km(ig, k) = l(ig, k)zqsm(ig, k) ! For momentum
	373	kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars
	374	kq(ig, k) = l(ig, k)zq0.2 ! For TKE
	375	END DO
	376	END DO
	377
	378
	379	!====================================================================
	380	! Transport diffusif vertical de la TKE par la TKE
	381	!====================================================================
	382
	383
	384	! initialize near-surface and top-layer mixing coefficients
	385	!...........................................................
	386
	387	kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
	388	kq(1:ngrid, klev+1) = 0 ! zero at the top
	389
	390	! Transport diffusif vertical de la TKE.
	391	!.......................................
	392
	393	IF (iflag_pbl>=12) THEN
	394	q2(1:ngrid, 1) = q2(1:ngrid, 2)
	395	CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
	396	END IF
	397
	398
	399	!====================================================================
	400	! Traitement particulier pour les cas tres stables, introduction d'une
	401	! longueur de m??lange minimale
	402	!====================================================================
	403	!
	404	! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
	405	! Holtslag A.A.M. and Boville B.A.
	406	! J. Clim., 6, 1825-1842, 1993
	407
	408
	409	IF (hboville) THEN
	410
	411
	412	IF (prt_level>1) THEN
	413	PRINT *, 'YAMADA4 0'
	414	END IF
	415
	416	DO ig = 1, ngrid
	417	coriol(ig) = 1.E-4
	418	pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
	419	END DO
	420
	421	IF (1==1) THEN
	422	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
	423
	424	DO k = 2, klev
	425	DO ig = 1, ngrid
	426	IF (teta(ig,2)>teta(ig,1)) THEN
	427	qmin = ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	428	kmin = kapzlev(ig, k)qmin
	429	ELSE
	430	kmin = -1. ! kmin n'est utilise que pour les SL stables.
	431	END IF
	432	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
	433
	434	kn(ig, k) = kmin
	435	km(ig, k) = kmin
	436	kq(ig, k) = kmin
	437
	438	! la longueur de melange est suposee etre l= kap z
	439	! K=l q Sm d'ou q2=(K/l Sm)**2
	440
	441	q2(ig, k) = (qmin/sm(ig,k))**2
	442	END IF
	443	END DO
	444	END DO
	445
	446	ELSE
	447	DO k = 2, klev
	448	DO ig = 1, ngrid
	449	IF (teta(ig,2)>teta(ig,1)) THEN
	450	qmin = ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	451	kmin = kapzlev(ig, k)qmin
	452	ELSE
	453	kmin = -1. ! kmin n'est utilise que pour les SL stables.
	454	END IF
	455	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
	456	kn(ig, k) = kmin
	457	km(ig, k) = kmin
	458	kq(ig, k) = kmin
	459	! la longueur de melange est suposee etre l= kap z
	460	! K=l q Sm d'ou q2=(K/l Sm)**2
	461	sm(ig, k) = 1.
	462	alpha(ig, k) = 1.
	463	q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
	464	zq = sqrt(q2(ig,k))
	465	km(ig, k) = l(ig, k)zqsm(ig, k)
	466	kn(ig, k) = km(ig, k)*alpha(ig, k)
	467	kq(ig, k) = l(ig, k)zq0.2
	468	END IF
	469	END DO
	470	END DO
	471	END IF
	472
	473	END IF
	474
	475	END IF ! hboville
	476
	477	IF (prt_level>1) THEN
	478	PRINT *, 'YAMADA4 1'
	479	END IF !(prt_level>1) THEN
	480
	481
	482	!======================================================
	483	! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	484	!======================================================
	485	!
	486	! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
	487	! Abdella K and McFarlane N
	488	! J. Atmos. Sci., 54, 1850-1867, 1997
	489
	490	! Diagnostique pour stokage
	491	!..........................
	492
	493	IF (1==0) THEN
	494	rino = rif
	495	smyam(1:ngrid, 1) = 0.
	496	styam(1:ngrid, 1) = 0.
	497	lyam(1:ngrid, 1) = 0.
	498	knyam(1:ngrid, 1) = 0.
	499	w2yam(1:ngrid, 1) = 0.
	500	t2yam(1:ngrid, 1) = 0.
	501
	502	smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
	503	styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
	504	lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
	505	knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
	506
	507
	508	! Calcul de w'2 et T'2
	509	!.......................
	510
	511	w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
	512	lyam(1:ngrid, 2:klev)5.17kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
	513	sqrt(q2(1:ngrid,2:klev))
	514
	515	t2yam(1:ngrid, 2:klev) = 9.1kn(1:ngrid, 2:klev) &
	516	dtetadz(1:ngrid, 2:klev)*2/sqrt(q2(1:ngrid,2:klev)) &
	517	lyam(1:ngrid, 2:klev)
	518	END IF
	519
	520	!============================================================================
	521
	522	first = .FALSE.
	523	RETURN
	524
	525
	526	END SUBROUTINE yamada4
	527
	528	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	529
	530	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	531	SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
	532
	533	USE dimphy
	534	IMPLICIT NONE
	535
	536	include "dimensions.h"
	537
	538	! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
	539	! avec un schema implicite en temps avec
	540	! inversion d'un syst??me tridiagonal
	541	!
	542	! Reference: Description of the interface with the surface and
	543	! the computation of the turbulet diffusion in LMDZ
	544	! Technical note on LMDZ
	545	! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
	546	!
	547	!============================================================================
	548	! Declarations
	549	!============================================================================
	550
	551	REAL plev(klon, klev+1)
	552	REAL temp(klon, klev)
	553	REAL timestep
	554	REAL gravity, rconst
	555	REAL kstar(klon, klev+1), zz
	556	REAL kmy(klon, klev+1)
	557	REAL q2(klon, klev+1)
	558	REAL deltap(klon, klev+1)
	559	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
	560	INTEGER ngrid
	561
	562	INTEGER i, k
	563
	564
	565	!=========================================================================
	566	! Calcul
	567	!=========================================================================
	568
	569	DO k = 1, klev
	570	DO i = 1, ngrid
	571	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	572	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
	573	(plev(i,k)-plev(i,k+1))*timestep
	574	END DO
	575	END DO
	576
	577	DO k = 2, klev
	578	DO i = 1, ngrid
	579	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
	580	END DO
	581	END DO
	582	DO i = 1, ngrid
	583	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
	584	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
	585	denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
	586	alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
	587	beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
	588	END DO
	589
	590	DO k = klev, 2, -1
	591	DO i = 1, ngrid
	592	denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
	593	kstar(i, k-1)
	594	alpha(i, k) = (q2(i,k)deltap(i,k)+kstar(i,k)alpha(i,k+1))/denom(i, k)
	595	beta(i, k) = kstar(i, k-1)/denom(i, k)
	596	END DO
	597	END DO
	598
	599	! Si on recalcule q2(1)
	600	!.......................
	601	IF (1==0) THEN
	602	DO i = 1, ngrid
	603	denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
	604	q2(i, 1) = (q2(i,1)deltap(i,1)+kstar(i,1)alpha(i,2))/denom(i, 1)
	605	END DO
	606	END IF
	607
	608
	609	DO k = 2, klev + 1
	610	DO i = 1, ngrid
	611	q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
	612	END DO
	613	END DO
	614
	615	RETURN
	616	END SUBROUTINE vdif_q2
	617	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	618
	619
	620
	621	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	622	SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
	623
	624	USE dimphy
	625	IMPLICIT NONE
	626
	627	include "dimensions.h"
	628	!
	629	! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
	630	! avec un schema explicite en temps
	631
	632
	633	!====================================================
	634	! Declarations
	635	!====================================================
	636
	637	REAL plev(klon, klev+1)
	638	REAL temp(klon, klev)
	639	REAL timestep
	640	REAL gravity, rconst
	641	REAL kstar(klon, klev+1), zz
	642	REAL kmy(klon, klev+1)
	643	REAL q2(klon, klev+1)
	644	REAL deltap(klon, klev+1)
	645	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
	646	INTEGER ngrid
	647	INTEGER i, k
	648
	649
	650	!==================================================
	651	! Calcul
	652	!==================================================
	653
	654	DO k = 1, klev
	655	DO i = 1, ngrid
	656	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	657	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
	658	(plev(i,k)-plev(i,k+1))*timestep
	659	END DO
	660	END DO
	661
	662	DO k = 2, klev
	663	DO i = 1, ngrid
	664	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
	665	END DO
	666	END DO
	667	DO i = 1, ngrid
	668	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
	669	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
	670	END DO
	671
	672	DO k = klev, 2, -1
	673	DO i = 1, ngrid
	674	q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
	675	k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
	676	END DO
	677	END DO
	678
	679	DO i = 1, ngrid
	680	q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
	681	q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
	682	klev)))/deltap(i, klev+1)
	683	END DO
	684
	685	RETURN
	686	END SUBROUTINE vdif_q2e
	687
	688	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	689
	690
	691	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	692
	693	SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
	694
	695
	696
	697	USE dimphy
	698	USE phys_state_var_mod, only: zstd, zsig, zmea
	699	USE phys_local_var_mod, only: l_mixmin, l_mix
	700
	701	! zstd: ecart type de la'altitud e sous-maille
	702	! zmea: altitude moyenne sous maille
	703	! zsig: pente moyenne de le maille
	704
	705	USE geometry_mod, only: cell_area
	706	! aire_cell: aire de la maille
	707
	708	IMPLICIT NONE
	709	!*************************************************************************
	710	! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
	711	! avec la formule de Blackadar
	712	! Calcul d'un minimum en fonction de l'orographie sous-maille:
	713	! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
	714	! induit par les circulations meso et submeso ??chelles.
	715	!
	716	! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
	717	! Blackadar A.K.
	718	! J. Geophys. Res., 64, No 8, 1962
	719	!
	720	! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
	721	! to large eddy simulations
	722	! Ayotte K et al
	723	! Boundary Layer Meteorology, 79, 131-175, 1996
	724	!
	725	!
	726	! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
	727	! Van de Wiel B.J.H et al
	728	! Boundary-Lay Meteorol, 128, 103-166, 2008
	729	!
	730	!
	731	! Histoire:
	732	!----------
	733	! * premi??re r??daction, Etienne et Frederic, 09/06/2016
	734	!
	735	! ***********************************************************************
	736
	737	!==================================================================
	738	! Declarations
	739	!==================================================================
	740
	741	! Inputs
	742	!-------
	743	INTEGER ni(klon) ! indice sur la grille original (non restreinte)
	744	INTEGER nsrf ! Type de surface
	745	INTEGER ngrid ! Nombre de points concern??s sur l'horizontal
	746	INTEGER iflag_pbl ! Choix du sch??ma de turbulence
	747	REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum
	748	REAL lmixmin ! Minimum absolu de la longueur de m??lange
	749	REAL zlay(klon, klev) ! altitude du centre de la couche
	750	REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
	751	REAL u(klon, klev) ! vitesse du vent zonal
	752	REAL v(klon, klev) ! vitesse du vent meridional
	753	REAL q2(klon, klev+1) ! energie cin??tique turbulente
	754	REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala
	755
	756	!In/out
	757	!-------
	758
	759	LOGICAL, SAVE :: firstcall = .TRUE.
	760	!$OMP THREADPRIVATE(firstcall)
	761
	762	! Outputs
	763	!---------
	764
	765	REAL lmix(klon, klev+1) ! Longueur de melange
	766
	767
	768	! Local
	769	!-------
	770
	771	INTEGER ig,jg, k
	772	REAL h_oro(klon)
	773	REAL hlim(klon)
	774	REAL, SAVE :: kap=0.4,kapb=0.4
	775	REAL zq
	776	REAL sq(klon), sqz(klon)
	777	REAL, ALLOCATABLE, SAVE :: l0(:)
	778	!$OMP THREADPRIVATE(l0)
	779	REAL fl, zzz, zl0, zq2, zn2
	780	REAL famorti, zzzz, zh_oro, zhlim
	781	REAL l1(klon, klev+1), l2(klon,klev+1)
	782	REAL winds(klon, klev)
	783	REAL xcell
	784	REAL zstdslope(klon)
	785	REAL lmax
	786	REAL l2strat, l2neutre, extent
	787	REAL l2limit(klon)
	788	!===============================================================
	789	! Fonctions utiles
	790	!===============================================================
	791
	792	! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
	793	!..........................................................................
	794
	795	fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)kapzlev(ig, &
	796	k)/(kapzlev(ig,k)+l0(ig)),0.5sqrt(q2(ig,k))/sqrt( &
	797	max(n2(ig,k),1.E-10))), 1.E-5)
	798
	799	! Fonction d'amortissement de la turbulence au dessus de la montagne
	800	! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
	801	!.....................................................................
	802
	803	famorti(zzzz, zh_oro, zhlim)=(-1.)ATAN((zzzz-zh_oro)/(zhlim-zh_oro))2./3.1416+1.
	804
	805	IF (ngrid==0) RETURN
	806
	807	IF (firstcall) THEN
	808	ALLOCATE (l0(klon))
	809	firstcall = .FALSE.
	810	END IF
	811
	812
	813	!=====================================================================
	814	! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
	815	!=====================================================================
	816
	817
	818	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
	819
	820
	821	! Iterative computation of l0
	822	! This version is kept for iflag_pbl only for convergence
	823	! with NPv3.1 Cmip5 simulations
	824	!...................................................................
	825
	826	DO ig = 1, ngrid
	827	sq(ig) = 1.E-10
	828	sqz(ig) = 1.E-10
	829	END DO
	830	DO k = 2, klev - 1
	831	DO ig = 1, ngrid
	832	zq = sqrt(q2(ig,k))
	833	sqz(ig) = sqz(ig) + zqzlev(ig, k)(zlay(ig,k)-zlay(ig,k-1))
	834	sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
	835	END DO
	836	END DO
	837	DO ig = 1, ngrid
	838	l0(ig) = 0.2*sqz(ig)/sq(ig)
	839	END DO
	840	DO k = 2, klev
	841	DO ig = 1, ngrid
	842	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
	843	END DO
	844	END DO
	845
	846	ELSE
	847
	848
	849	! In all other case, the assymptotic mixing length l0 is imposed (150m)
	850	!......................................................................
	851
	852	l0(1:ngrid) = 150.
	853	DO k = 2, klev
	854	DO ig = 1, ngrid
	855	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
	856	END DO
	857	END DO
	858
	859	END IF
	860
	861	!=================================================================================
	862	! CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
	863	! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
	864	! glacier, la glace de mer et les oc??ans)
	865	!=================================================================================
	866
	867	l2(1:ngrid,:)=0.0
	868	l_mixmin(1:ngrid,:,nsrf)=0.
	869	l_mix(1:ngrid,:,nsrf)=0.
	870
	871	IF (nsrf .EQ. 1) THEN
	872
	873	! coefficients
	874	!--------------
	875
	876	extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1.
	877	lmax=150. ! Longueur de m??lange max dans l'absolu
	878
	879	! calculs
	880	!---------
	881
	882	DO ig=1,ngrid
	883
	884	! On calcule la hauteur du relief
	885	!.................................
	886	! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
	887	! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
	888	! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
	889	! (en gros, une maille de taille xcell avec une pente constante zstdslope)
	890	jg=ni(ig)
	891	! IF (zsig(jg) .EQ. 0.) THEN
	892	! zstdslope(ig)=0.
	893	! ELSE
	894	! xcell=sqrt(cell_area(jg))
	895	! zstdslope(ig)=max((xcellzsig(jg)-zmea(jg))3 /(3.zsig(jg)),0.)
	896	! zstdslope(ig)=sqrt(zstdslope(ig))
	897	! END IF
	898
	899	! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief
	900	h_oro(ig)=zstd(jg)
	901	hlim(ig)=extent*h_oro(ig)
	902	ENDDO
	903
	904	l2limit(1:ngrid)=0.
	905
	906	DO k=2,klev
	907	DO ig=1,ngrid
	908	winds(ig,k)=sqrt(u(ig,k)2+v(ig,k)2)
	909	IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie
	910	l2strat= kapbpbl_lmixmin_alphawinds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
	911	l2neutre=kapzlev(ig,k)h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h
	912	l2neutre=MIN(l2neutre,lmax) ! On majore par lmax
	913	l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
	914	l2(ig,k)=l2limit(ig)
	915
	916	ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
	917
	918	! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
	919	! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
	920	! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
	921	l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
	922	ELSE ! Au dessus de extent*h, on prend l2=l0
	923	l2(ig,k)=0.
	924	END IF
	925	ENDDO
	926	ENDDO
	927	ENDIF ! pbl_lmixmin_alpha
	928
	929	!==================================================================================
	930	! On prend le max entre la longueur de melange de blackadar et celle calcul??e
	931	! en fonction de la topographie
	932	!===================================================================================
	933
	934
	935	DO k=2,klev
	936	DO ig=1,ngrid
	937	lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
	938	ENDDO
	939	ENDDO
	940
	941	! Diagnostics
	942
	943	DO k=2,klev
	944	DO ig=1,ngrid
	945	jg=ni(ig)
	946	l_mix(jg,k,nsrf)=lmix(ig,k)
	947	l_mixmin(jg,k,nsrf)=l2(ig,k)
	948	ENDDO
	949	ENDDO
	950	DO ig=1,ngrid
	951	jg=ni(ig)
	952	l_mix(jg,1,nsrf)=hlim(ig)
	953	ENDDO
	954
	955
	956
	957	END SUBROUTINE mixinglength

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

Download in other formats:

Original Format