Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

yamada4.F90 @ 5225

Last change on this file since 5225 was 2924, checked in by fcheruy, 7 years ago

Creation of LMDZ branch to incorporate tree drag from ORCHIDEE.
Should merge in LMDZ trunk quickly

Cheruy with assistance from LF

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 36.2 KB

Rev	Line
[2561]	1	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	2
[2561]	3	SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
[2924]	4	cd, tke, km, kn, kq, ustar, iflag_pbl, drgpro)
[541]	5
[1992]	6	USE dimphy
[2573]	7	USE ioipsl_getin_p_mod, ONLY : getin_p
	8
[1992]	9	IMPLICIT NONE
[2561]	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
[1992]	28	! dt : pas de temps
	29	! g : g
[2561]	30	! rconst: constante de l'air sec
[1992]	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)
[2561]	36	! teta : temperature potentielle virtuelle au centre de chaque couche
[1992]	37	! (en entree : la valeur au debut du pas de temps)
[2561]	38	! cd : cdrag pour la quantit?? de mouvement
[1992]	39	! (en entree : la valeur au debut du pas de temps)
[2561]	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.
[2721]	58	! 2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
[2817]	59	! On met tke (=q2/2) en entr??e plut??t que q2
[2721]	60	! On corrige l'update de la tke
[2561]	61	!
	62	! Inpout/Output :
	63	!==============
[2721]	64	! tke : tke au bas de chaque couche
[1992]	65	! (en entree : la valeur au debut du pas de temps)
	66	! (en sortie : la valeur a la fin du pas de temps)
[2561]	67
	68	! Outputs:
	69	!==========
[1992]	70	! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	71	! couche)
	72	! (en sortie : la valeur a la fin du pas de temps)
	73	! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	74	! (en sortie : la valeur a la fin du pas de temps)
[2561]	75	!
	76	!.......................................................................
[541]	77
[2561]	78	!=======================================================================
	79	! Declarations:
	80	!=======================================================================
[541]	81
	82
[2561]	83	! Inputs/Outputs
	84	!----------------
[1992]	85	REAL dt, g, rconst
	86	REAL plev(klon, klev+1), temp(klon, klev)
	87	REAL ustar(klon)
	88	REAL kmin, qmin, pblhmin(klon), coriol(klon)
	89	REAL zlev(klon, klev+1)
	90	REAL zlay(klon, klev)
	91	REAL u(klon, klev)
	92	REAL v(klon, klev)
	93	REAL teta(klon, klev)
	94	REAL cd(klon)
[2721]	95	REAL tke(klon, klev+1)
[1992]	96	REAL unsdz(klon, klev)
	97	REAL unsdzdec(klon, klev+1)
[2561]	98	REAL kn(klon, klev+1)
	99	REAL km(klon, klev+1)
	100	INTEGER iflag_pbl, ngrid, nsrf
	101	INTEGER ni(klon)
[541]	102
[2924]	103	!FC
	104	REAL drgpro(klon,klev)
	105	REAL winds(klon,klev)
[2561]	106
	107	! Local
	108	!-------
	109
	110	INCLUDE "clesphys.h"
	111
[2721]	112	REAL q2(klon, klev+1)
[2561]	113	REAL kmpre(klon, klev+1), tmp2, qpre
[1992]	114	REAL mpre(klon, klev+1)
	115	REAL kq(klon, klev+1)
	116	REAL ff(klon, klev+1), delta(klon, klev+1)
	117	REAL aa(klon, klev+1), aa0, aa1
	118	INTEGER nlay, nlev
	119	LOGICAL first
	120	INTEGER ipas
	121	SAVE first, ipas
	122	! FH/IM data first,ipas/.true.,0/
	123	DATA first, ipas/.FALSE., 0/
	124	!$OMP THREADPRIVATE( first,ipas)
[2828]	125	LOGICAL,SAVE :: hboville=.TRUE.
	126	REAL,SAVE :: viscom,viscoh
	127	!$OMP THREADPRIVATE( hboville,viscom,viscoh)
[1992]	128	INTEGER ig, k
	129	REAL ri, zrif, zalpha, zsm, zsn
	130	REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
	131	REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
	132	REAL dtetadz(klon, klev+1)
	133	REAL m2cstat, mcstat, kmcstat
	134	REAL l(klon, klev+1)
[2561]	135	REAL zz(klon, klev+1)
[1992]	136	INTEGER iter
[2828]	137	REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
[2721]	138	REAL,SAVE :: ric0,ric,rifc, b1, kap
	139	!$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
	140	DATA b1, kap/16.6, 0.4/
	141	REAL,SAVE :: seuilsm, seuilalpha
	142	!$OMP THREADPRIVATE(seuilsm, seuilalpha)
[2573]	143	REAL,SAVE :: lmixmin
	144	!$OMP THREADPRIVATE(lmixmin)
[2721]	145	LOGICAL, SAVE :: new_yamada4
[2817]	146	INTEGER, SAVE :: yamada4_num
	147	!$OMP THREADPRIVATE(new_yamada4,yamada4_num)
[2891]	148	REAL, SAVE :: yun,ydeux
	149	REAL :: disseff
[2721]	150	!$OMP THREADPRIVATE(yun,ydeux)
[1992]	151	REAL frif, falpha, fsm
	152	REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
	153	lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
	154	LOGICAL, SAVE :: firstcall = .TRUE.
	155	!$OMP THREADPRIVATE(firstcall)
[2561]	156
	157
[2721]	158
[2561]	159	! Fonctions utilis??es
	160	!--------------------
	161
[1992]	162	frif(ri) = 0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	163	falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
	164	fsm(ri) = 1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
[2561]	165
[541]	166
[2573]	167	IF (firstcall) THEN
[2721]	168	! Seuil dans le code de turbulence
	169	new_yamada4=.false.
	170	CALL getin_p('new_yamada4',new_yamada4)
[2828]	171
[2721]	172	IF (new_yamada4) THEN
[2828]	173	! Corrections et reglages issus du travail de these d'Etienne Vignon.
[2721]	174	ric=0.143 ! qui donne des valeurs proches des seuils proposes
	175	! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
	176	! sm=1.1213 (1.12 dans Y83)
	177	CALL getin_p('yamada4_ric',ric)
	178	ric0=0.19489 ! ric=0.195 originalement, mais produisait sm<0
[2817]	179	ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
[2721]	180	rifc=frif(ric)
	181	seuilsm=fsm(frif(ric))
	182	seuilalpha=falpha(frif(ric))
	183	yun=1.
	184	ydeux=2.
[2828]	185	hboville=.FALSE.
	186	viscom=1.46E-5
	187	viscoh=2.06E-5
	188	lmixmin=0.
	189	yamada4_num=5
[2721]	190	ELSE
	191	ric=0.195
	192	rifc=0.191
	193	seuilalpha=1.12
	194	seuilsm=0.085
	195	yun=2.
	196	ydeux=1.
[2828]	197	hboville=.TRUE.
	198	viscom=0.
	199	viscoh=0.
	200	lmixmin=1.
	201	yamada4_num=0
[2721]	202	ENDIF
[2828]	203
[2721]	204	PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
[2573]	205	firstcall = .FALSE.
	206	CALL getin_p('lmixmin',lmixmin)
[2828]	207	CALL getin_p('yamada4_hboville',hboville)
	208	CALL getin_p('yamada4_num',yamada4_num)
[2573]	209	END IF
	210
	211
	212
[2561]	213	!===============================================================================
	214	! Flags, tests et ??valuations de constantes
	215	!===============================================================================
[541]	216
[2561]	217	! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
[1738]	218
[541]	219
[1992]	220	IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
	221	STOP 'probleme de coherence dans appel a MY'
	222	END IF
[541]	223
[2561]	224
	225	nlay = klev
	226	nlev = klev + 1
[1992]	227	ipas = ipas + 1
[541]	228
	229
[2561]	230	!========================================================================
	231	! Calcul des increments verticaux
	232	!=========================================================================
[541]	233
[2561]	234
	235	! Attention: zlev n'est pas declare a nlev
[1992]	236	DO ig = 1, ngrid
	237	zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
	238	END DO
[541]	239
[2561]	240
[1992]	241	DO k = 1, nlay
	242	DO ig = 1, ngrid
	243	unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	244	END DO
	245	END DO
	246	DO ig = 1, ngrid
	247	unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
	248	END DO
	249	DO k = 2, nlay
	250	DO ig = 1, ngrid
	251	unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	252	END DO
	253	END DO
	254	DO ig = 1, ngrid
	255	unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	256	END DO
[1738]	257
[2561]	258	!=========================================================================
	259	! Richardson number and stability functions
	260	!=========================================================================
	261
	262	! initialize arrays:
	263
[2574]	264	m2(1:ngrid, :) = 0.0
	265	sm(1:ngrid, :) = 0.0
	266	rif(1:ngrid, :) = 0.0
[1738]	267
[2561]	268	!------------------------------------------------------------
[1992]	269	DO k = 2, klev
[2561]	270
[1992]	271	DO ig = 1, ngrid
	272	dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
	273	m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
	274	k-1))*2)/(dz(ig,k)dz(ig,k))
	275	dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
	276	n2(ig, k) = g2.dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
	277	ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
	278	IF (ri<ric) THEN
	279	rif(ig, k) = frif(ri)
	280	ELSE
	281	rif(ig, k) = rifc
	282	END IF
[2721]	283	if (new_yamada4) then
	284	alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
	285	sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
	286	else
[1992]	287	IF (rif(ig,k)<0.16) THEN
	288	alpha(ig, k) = falpha(rif(ig,k))
	289	sm(ig, k) = fsm(rif(ig,k))
	290	ELSE
[2561]	291	alpha(ig, k) = seuilalpha
	292	sm(ig, k) = seuilsm
[1992]	293	END IF
[2721]	294
	295	end if
[1992]	296	zz(ig, k) = b1m2(ig, k)(1.-rif(ig,k))*sm(ig, k)
	297	END DO
	298	END DO
[1738]	299
	300
	301
[2721]	302
	303
	304	!=======================================================================
	305	! DIFFERENT TYPES DE SCHEMA de YAMADA
	306	!=======================================================================
	307
[2817]	308	! On commence par calculer q2 a partir de la tke
[2721]	309
	310	IF (new_yamada4) THEN
	311	DO k=1,klev+1
	312	q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
	313	ENDDO
	314	ELSE
	315	DO k=1,klev+1
	316	q2(1:ngrid,k)=tke(1:ngrid,k)
	317	ENDDO
	318	ENDIF
	319
[2561]	320	! ====================================================================
	321	! Computing the mixing length
	322	! ====================================================================
[541]	323
[2561]	324
[2573]	325	CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
[541]	326
	327
[2561]	328	!--------------
[1992]	329	! Yamada 2.0
[2561]	330	!--------------
[1992]	331	IF (iflag_pbl==6) THEN
[2561]	332
[1992]	333	DO k = 2, klev
[2574]	334	q2(1:ngrid, k) = l(1:ngrid, k)*2zz(1:ngrid, k)
[1992]	335	END DO
	336
	337
[2561]	338	!------------------
	339	! Yamada 2.Fournier
	340	!------------------
	341
[1992]	342	ELSE IF (iflag_pbl==7) THEN
	343
[2561]	344
[1992]	345	! Calcul de l, km, au pas precedent
[2561]	346	!....................................
[1992]	347	DO k = 2, klev
	348	DO ig = 1, ngrid
	349	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
	350	kmpre(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
	351	mpre(ig, k) = sqrt(m2(ig,k))
	352	END DO
	353	END DO
	354
	355	DO k = 2, klev - 1
	356	DO ig = 1, ngrid
	357	m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
	358	mcstat = sqrt(m2cstat)
	359
[2561]	360	! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
	361	!.........................................................................
[1992]	362
	363	IF (k==2) THEN
	364	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
	365	unsdz(ig,k-1)cd(ig)(sqrt(u(ig,3)2+v(ig,3)2)-mcstat/unsdzdec &
	366	(ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
	367	-1))
[541]	368	ELSE
[1992]	369	kmcstat = 1.E+0/mcstat(unsdz(ig,k)kmpre(ig,k+1)*mpre(ig,k+1)+ &
	370	unsdz(ig,k-1)kmpre(ig,k-1)mpre(ig,k-1))/ &
	371	(unsdz(ig,k)+unsdz(ig,k-1))
	372	END IF
[2561]	373
[1992]	374	tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
	375	q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
[541]	376
[1992]	377	END DO
	378	END DO
[541]	379
[2561]	380
	381	! ------------------------
[1992]	382	! Yamada 2.5 a la Didi
[2561]	383	!-------------------------
[541]	384
[2561]	385	ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
[541]	386
[2561]	387	! Calcul de l, km, au pas precedent
	388	!....................................
[1992]	389	DO k = 2, klev
	390	DO ig = 1, ngrid
	391	delta(ig, k) = q2(ig, k)/(l(ig,k)*2sm(ig,k))
	392	IF (delta(ig,k)<1.E-20) THEN
	393	delta(ig, k) = 1.E-20
	394	END IF
	395	km(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
	396	aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	397	aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	398	aa(ig, k) = aa1dt/(delta(ig,k)l(ig,k))
	399	qpre = sqrt(q2(ig,k))
	400	IF (aa(ig,k)>0.) THEN
	401	q2(ig, k) = (qpre+aa(ig,k)qpreqpre)**2
	402	ELSE
	403	q2(ig, k) = (qpre/(1.-aa(ig,k)qpre))*2
	404	END IF
	405	! else ! iflag_pbl=9
	406	! if (aa(ig,k)*qpre.gt.0.9) then
	407	! q2(ig,k)=(qpre10.)*2
	408	! else
	409	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	410	! endif
	411	! endif
	412	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
	413	END DO
	414	END DO
[1738]	415
[1992]	416	ELSE IF (iflag_pbl>=10) THEN
[1738]	417
[2817]	418	IF (yamada4_num>=1) THEN
	419
[1992]	420	DO k = 2, klev - 1
[2817]	421	DO ig=1,ngrid
	422	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
	423	km(ig, k) = l(ig, k)sqrt(q2(ig,k))sm(ig, k)
	424	shear(ig,k)=km(ig, k)*m2(ig, k)
	425	buoy(ig,k)=km(ig, k)m2(ig, k)(-1.*rif(ig,k))
	426	dissip(ig,k)=((sqrt(q2(ig,k)))*3)/(b1l(ig,k))
	427	ENDDO
	428	ENDDO
	429
	430	IF (yamada4_num==1) THEN ! Schema du MAR tel quel
	431	DO k = 2, klev - 1
	432	DO ig=1,ngrid
	433	tkeprov=q2(ig,k)/ydeux
[2828]	434	tkeprov= tkeprov* &
[2817]	435	& (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
	436	& (tkeprov+dt((-1.)min(0.,buoy(ig,k))+dissip(ig,k)))
[2828]	437	q2(ig,k)=tkeprov*ydeux
[2817]	438	ENDDO
	439	ENDDO
[2828]	440	ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
[2817]	441	DO k = 2, klev - 1
	442	DO ig=1,ngrid
	443	tkeprov=q2(ig,k)/ydeux
[2828]	444	disseff=dissip(ig,k)-min(0.,buoy(ig,k))
	445	tkeprov = tkeprov/(1.+dtdisseff/(2.tkeprov))**2
	446	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
	447	q2(ig,k)=tkeprov*ydeux
[2817]	448	! En cas stable, on traite la flotabilite comme la
	449	! dissipation, en supposant que buoy/q2^3 est constant.
	450	! Puis on prend la solution exacte
	451	ENDDO
	452	ENDDO
[2828]	453	ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
	454	DO k = 2, klev - 1
	455	DO ig=1,ngrid
	456	tkeprov=q2(ig,k)/ydeux
	457	disseff=dissip(ig,k)-min(0.,buoy(ig,k))
	458	tkeprov=tkeprovexp(-dtdisseff/tkeprov)
	459	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
	460	q2(ig,k)=tkeprov*ydeux
	461	! En cas stable, on traite la flotabilite comme la
	462	! dissipation, en supposant que buoy/q2^3 est constant.
	463	! Puis on prend la solution exacte
	464	ENDDO
	465	ENDDO
	466	ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
	467	DO k = 2, klev - 1
	468	DO ig=1,ngrid
	469	tkeprov=q2(ig,k)/ydeux
	470	tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
	471	tkeprov= tkeprov* &
	472	& tkeprov/ &
	473	& (tkeprov+dt((-1.)min(0.,buoy(ig,k))+dissip(ig,k)))
	474	q2(ig,k)=tkeprov*ydeux
	475	! En cas stable, on traite la flotabilite comme la
	476	! dissipation, en supposant que buoy/q2^3 est constant.
	477	! Puis on prend la solution exacte
	478	ENDDO
	479	ENDDO
	480	ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
	481	DO k = 2, klev - 1
	482	DO ig=1,ngrid
	483	tkeprov=q2(ig,k)/ydeux
[2924]	484
	485	! if(ifl_pbltree .eq. 0) then
	486	! disseff=dissip(ig,k)-min(0.,buoy(ig,k))
	487	! tkeexp=exp(-dt*disseff/tkeprov)
	488	! tkeprov= shear(ig,k)tkeprov/disseff(1.-tkeexp)+tkeprov*tkeexp
	489	! else
	490	!FC on ajoute la dissipation due aux arbres
	491	disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov
[2828]	492	tkeexp=exp(-dt*disseff/tkeprov)
[2924]	493	! on prend en compte la tke cree par les arbres
	494	winds(ig,k)=sqrt(u(ig,k)2+v(ig,k)2)
	495	tkeprov= (shear(ig,k)+ &
	496	& drgpro(ig,k)(winds(ig,k))3)tkeprov/disseff(1.-tkeexp)+tkeprovtkeexp
	497	! endif
	498
[2828]	499	q2(ig,k)=tkeprov*ydeux
[2924]	500
[2828]	501	! En cas stable, on traite la flotabilite comme la
	502	! dissipation, en supposant que buoy/q2^3 est constant.
	503	! Puis on prend la solution exacte
	504	ENDDO
	505	ENDDO
	506	ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
	507	DO k = 2, klev - 1
	508	DO ig=1,ngrid
	509	tkeprov=q2(ig,k)/ydeux
	510	tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
	511	disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
	512	tkeexp=exp(-dt*disseff/tkeprov)
	513	tkeprov= tkeprov*tkeexp
	514	q2(ig,k)=tkeprov*ydeux
	515	! En cas stable, on traite la flotabilite comme la
	516	! dissipation, en supposant que buoy/q2^3 est constant.
	517	! Puis on prend la solution exacte
	518	ENDDO
	519	ENDDO
[2817]	520	ENDIF
	521
	522	DO k = 2, klev - 1
	523	DO ig=1,ngrid
	524	q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
	525	ENDDO
	526	ENDDO
	527
	528	ELSE
	529
	530	DO k = 2, klev - 1
[2574]	531	km(1:ngrid, k) = l(1:ngrid, k)sqrt(q2(1:ngrid,k))sm(1:ngrid, k)
[2721]	532	q2(1:ngrid, k) = q2(1:ngrid, k) + ydeuxdtkm(1:ngrid, k)m2(1:ngrid, k)(1.-rif(1:ngrid,k))
	533	! q2(1:ngrid, k) = q2(1:ngrid, k) + dtkm(1:ngrid, k)m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
[2574]	534	q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
[2721]	535	q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yunl(1:ngrid,k)b1))
	536	! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2l(1:ngrid,k)b1))
[2574]	537	q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
[1992]	538	END DO
[1738]	539
[2817]	540	ENDIF
[1738]	541
[1992]	542	ELSE
	543	STOP 'Cas nom prevu dans yamada4'
[541]	544
[1992]	545	END IF ! Fin du cas 8
[541]	546
	547
[1992]	548	! ====================================================================
[2561]	549	! Calcul des coefficients de melange
[1992]	550	! ====================================================================
[2561]	551
[1992]	552	DO k = 2, klev
	553	DO ig = 1, ngrid
	554	zq = sqrt(q2(ig,k))
[2561]	555	km(ig, k) = l(ig, k)zqsm(ig, k) ! For momentum
	556	kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars
	557	kq(ig, k) = l(ig, k)zq0.2 ! For TKE
[1992]	558	END DO
	559	END DO
[2561]	560
	561
	562	!====================================================================
	563	! Transport diffusif vertical de la TKE par la TKE
	564	!====================================================================
	565
	566
[1992]	567	! initialize near-surface and top-layer mixing coefficients
[2561]	568	!...........................................................
[541]	569
[2561]	570	kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
	571	kq(1:ngrid, klev+1) = 0 ! zero at the top
	572
	573	! Transport diffusif vertical de la TKE.
	574	!.......................................
	575
[1992]	576	IF (iflag_pbl>=12) THEN
[2574]	577	q2(1:ngrid, 1) = q2(1:ngrid, 2)
[1992]	578	CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
	579	END IF
[541]	580
	581
[2561]	582	!====================================================================
	583	! Traitement particulier pour les cas tres stables, introduction d'une
	584	! longueur de m??lange minimale
	585	!====================================================================
	586	!
	587	! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
	588	! Holtslag A.A.M. and Boville B.A.
	589	! J. Clim., 6, 1825-1842, 1993
[541]	590
[2561]	591
	592	IF (hboville) THEN
	593
	594
[1992]	595	IF (prt_level>1) THEN
	596	PRINT *, 'YAMADA4 0'
[2561]	597	END IF
	598
[1992]	599	DO ig = 1, ngrid
	600	coriol(ig) = 1.E-4
	601	pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
	602	END DO
[1738]	603
[1992]	604	IF (1==1) THEN
	605	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
[1738]	606
[1992]	607	DO k = 2, klev
	608	DO ig = 1, ngrid
	609	IF (teta(ig,2)>teta(ig,1)) THEN
	610	qmin = ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	611	kmin = kapzlev(ig, k)qmin
	612	ELSE
	613	kmin = -1. ! kmin n'est utilise que pour les SL stables.
	614	END IF
	615	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
[2561]	616
[1992]	617	kn(ig, k) = kmin
	618	km(ig, k) = kmin
	619	kq(ig, k) = kmin
[2561]	620
	621	! la longueur de melange est suposee etre l= kap z
	622	! K=l q Sm d'ou q2=(K/l Sm)**2
	623
[1992]	624	q2(ig, k) = (qmin/sm(ig,k))**2
	625	END IF
	626	END DO
	627	END DO
[1738]	628
[1992]	629	ELSE
	630	DO k = 2, klev
	631	DO ig = 1, ngrid
	632	IF (teta(ig,2)>teta(ig,1)) THEN
	633	qmin = ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	634	kmin = kapzlev(ig, k)qmin
	635	ELSE
	636	kmin = -1. ! kmin n'est utilise que pour les SL stables.
	637	END IF
	638	IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
	639	kn(ig, k) = kmin
	640	km(ig, k) = kmin
	641	kq(ig, k) = kmin
[2561]	642	! la longueur de melange est suposee etre l= kap z
	643	! K=l q Sm d'ou q2=(K/l Sm)**2
[1992]	644	sm(ig, k) = 1.
	645	alpha(ig, k) = 1.
	646	q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
	647	zq = sqrt(q2(ig,k))
	648	km(ig, k) = l(ig, k)zqsm(ig, k)
	649	kn(ig, k) = km(ig, k)*alpha(ig, k)
	650	kq(ig, k) = l(ig, k)zq0.2
	651	END IF
	652	END DO
	653	END DO
	654	END IF
[1738]	655
[1992]	656	END IF
[541]	657
[2561]	658	END IF ! hboville
	659
[2828]	660	! Ajout d'une viscosite moleculaire
	661	km(:,:)=km(:,:)+viscom
	662	kn(:,:)=kn(:,:)+viscoh
	663	kq(:,:)=kq(:,:)+viscoh
	664
[1992]	665	IF (prt_level>1) THEN
	666	PRINT *, 'YAMADA4 1'
	667	END IF !(prt_level>1) THEN
[2561]	668
	669
	670	!======================================================
	671	! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	672	!======================================================
	673	!
	674	! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
	675	! Abdella K and McFarlane N
	676	! J. Atmos. Sci., 54, 1850-1867, 1997
	677
[1992]	678	! Diagnostique pour stokage
[2561]	679	!..........................
[541]	680
[1992]	681	IF (1==0) THEN
	682	rino = rif
	683	smyam(1:ngrid, 1) = 0.
	684	styam(1:ngrid, 1) = 0.
	685	lyam(1:ngrid, 1) = 0.
	686	knyam(1:ngrid, 1) = 0.
	687	w2yam(1:ngrid, 1) = 0.
	688	t2yam(1:ngrid, 1) = 0.
[878]	689
[1992]	690	smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
	691	styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
	692	lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
	693	knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
[541]	694
	695
[2561]	696	! Calcul de w'2 et T'2
	697	!.......................
	698
[1992]	699	w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
	700	lyam(1:ngrid, 2:klev)5.17kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
	701	sqrt(q2(1:ngrid,2:klev))
[541]	702
[1992]	703	t2yam(1:ngrid, 2:klev) = 9.1kn(1:ngrid, 2:klev) &
	704	dtetadz(1:ngrid, 2:klev)*2/sqrt(q2(1:ngrid,2:klev)) &
	705	lyam(1:ngrid, 2:klev)
	706	END IF
[1403]	707
[2721]	708
	709
[2561]	710	!============================================================================
[2817]	711	! Mise a jour de la tke
[2721]	712	!============================================================================
[2561]	713
[2721]	714	IF (new_yamada4) THEN
	715	DO k=1,klev+1
	716	tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
	717	ENDDO
	718	ELSE
	719	DO k=1,klev+1
	720	tke(1:ngrid,k)=q2(1:ngrid,k)
	721	ENDDO
	722	ENDIF
	723
	724
	725	!============================================================================
	726
[1992]	727	first = .FALSE.
	728	RETURN
[2561]	729
	730
[1992]	731	END SUBROUTINE yamada4
[2561]	732
	733	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	734
	735	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1992]	736	SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
[2561]	737
[1992]	738	USE dimphy
	739	IMPLICIT NONE
[2561]	740
	741	include "dimensions.h"
[1403]	742
[2561]	743	! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
	744	! avec un schema implicite en temps avec
	745	! inversion d'un syst??me tridiagonal
	746	!
	747	! Reference: Description of the interface with the surface and
	748	! the computation of the turbulet diffusion in LMDZ
	749	! Technical note on LMDZ
	750	! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
	751	!
	752	!============================================================================
	753	! Declarations
	754	!============================================================================
[1403]	755
[1992]	756	REAL plev(klon, klev+1)
	757	REAL temp(klon, klev)
	758	REAL timestep
	759	REAL gravity, rconst
	760	REAL kstar(klon, klev+1), zz
	761	REAL kmy(klon, klev+1)
	762	REAL q2(klon, klev+1)
	763	REAL deltap(klon, klev+1)
	764	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
	765	INTEGER ngrid
[1403]	766
[1992]	767	INTEGER i, k
[1403]	768
[2561]	769
	770	!=========================================================================
	771	! Calcul
	772	!=========================================================================
	773
[1992]	774	DO k = 1, klev
	775	DO i = 1, ngrid
	776	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	777	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
	778	(plev(i,k)-plev(i,k+1))*timestep
	779	END DO
	780	END DO
[1403]	781
[1992]	782	DO k = 2, klev
	783	DO i = 1, ngrid
	784	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
	785	END DO
	786	END DO
	787	DO i = 1, ngrid
	788	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
	789	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
	790	denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
	791	alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
	792	beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
	793	END DO
[1403]	794
[1992]	795	DO k = klev, 2, -1
	796	DO i = 1, ngrid
	797	denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
	798	kstar(i, k-1)
	799	alpha(i, k) = (q2(i,k)deltap(i,k)+kstar(i,k)alpha(i,k+1))/denom(i, k)
	800	beta(i, k) = kstar(i, k-1)/denom(i, k)
	801	END DO
	802	END DO
[1403]	803
[1992]	804	! Si on recalcule q2(1)
[2561]	805	!.......................
[1992]	806	IF (1==0) THEN
	807	DO i = 1, ngrid
	808	denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
	809	q2(i, 1) = (q2(i,1)deltap(i,1)+kstar(i,1)alpha(i,2))/denom(i, 1)
	810	END DO
	811	END IF
[1403]	812
[2561]	813
[1992]	814	DO k = 2, klev + 1
	815	DO i = 1, ngrid
	816	q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
	817	END DO
	818	END DO
[1403]	819
[1992]	820	RETURN
	821	END SUBROUTINE vdif_q2
[2561]	822	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	823
	824
	825
	826	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	827	SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
	828
	829	USE dimphy
[1992]	830	IMPLICIT NONE
[1403]	831
[2561]	832	include "dimensions.h"
	833	!
	834	! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
	835	! avec un schema explicite en temps
[1403]	836
[2561]	837
	838	!====================================================
	839	! Declarations
	840	!====================================================
	841
[1992]	842	REAL plev(klon, klev+1)
	843	REAL temp(klon, klev)
	844	REAL timestep
	845	REAL gravity, rconst
	846	REAL kstar(klon, klev+1), zz
	847	REAL kmy(klon, klev+1)
	848	REAL q2(klon, klev+1)
	849	REAL deltap(klon, klev+1)
	850	REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
	851	INTEGER ngrid
	852	INTEGER i, k
[1403]	853
[2561]	854
	855	!==================================================
	856	! Calcul
	857	!==================================================
	858
[1992]	859	DO k = 1, klev
	860	DO i = 1, ngrid
	861	zz = (plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	862	kstar(i, k) = 0.125(kmy(i,k+1)+kmy(i,k))zz*zz/ &
	863	(plev(i,k)-plev(i,k+1))*timestep
	864	END DO
	865	END DO
[1403]	866
[1992]	867	DO k = 2, klev
	868	DO i = 1, ngrid
	869	deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
	870	END DO
	871	END DO
	872	DO i = 1, ngrid
	873	deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
	874	deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
	875	END DO
	876
	877	DO k = klev, 2, -1
	878	DO i = 1, ngrid
	879	q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
	880	k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
	881	END DO
	882	END DO
	883
	884	DO i = 1, ngrid
	885	q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
	886	q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
	887	klev)))/deltap(i, klev+1)
	888	END DO
	889
	890	RETURN
	891	END SUBROUTINE vdif_q2e
[2561]	892
	893	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	894
	895
	896	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	897
[2573]	898	SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
[2561]	899
	900
	901
	902	USE dimphy
	903	USE phys_state_var_mod, only: zstd, zsig, zmea
	904	USE phys_local_var_mod, only: l_mixmin, l_mix
	905
	906	! zstd: ecart type de la'altitud e sous-maille
	907	! zmea: altitude moyenne sous maille
	908	! zsig: pente moyenne de le maille
	909
	910	USE geometry_mod, only: cell_area
	911	! aire_cell: aire de la maille
	912
	913	IMPLICIT NONE
	914	!*************************************************************************
	915	! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
	916	! avec la formule de Blackadar
	917	! Calcul d'un minimum en fonction de l'orographie sous-maille:
	918	! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
	919	! induit par les circulations meso et submeso ??chelles.
	920	!
	921	! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
	922	! Blackadar A.K.
	923	! J. Geophys. Res., 64, No 8, 1962
	924	!
	925	! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
	926	! to large eddy simulations
	927	! Ayotte K et al
	928	! Boundary Layer Meteorology, 79, 131-175, 1996
	929	!
	930	!
	931	! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
	932	! Van de Wiel B.J.H et al
	933	! Boundary-Lay Meteorol, 128, 103-166, 2008
	934	!
	935	!
	936	! Histoire:
	937	!----------
	938	! * premi??re r??daction, Etienne et Frederic, 09/06/2016
	939	!
	940	! ***********************************************************************
	941
	942	!==================================================================
	943	! Declarations
	944	!==================================================================
	945
	946	! Inputs
	947	!-------
	948	INTEGER ni(klon) ! indice sur la grille original (non restreinte)
	949	INTEGER nsrf ! Type de surface
	950	INTEGER ngrid ! Nombre de points concern??s sur l'horizontal
	951	INTEGER iflag_pbl ! Choix du sch??ma de turbulence
	952	REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum
	953	REAL lmixmin ! Minimum absolu de la longueur de m??lange
	954	REAL zlay(klon, klev) ! altitude du centre de la couche
	955	REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
	956	REAL u(klon, klev) ! vitesse du vent zonal
	957	REAL v(klon, klev) ! vitesse du vent meridional
	958	REAL q2(klon, klev+1) ! energie cin??tique turbulente
	959	REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala
	960
	961	!In/out
	962	!-------
	963
[2573]	964	LOGICAL, SAVE :: firstcall = .TRUE.
	965	!$OMP THREADPRIVATE(firstcall)
[2561]	966
	967	! Outputs
	968	!---------
	969
	970	REAL lmix(klon, klev+1) ! Longueur de melange
	971
	972
	973	! Local
	974	!-------
	975
	976	INTEGER ig,jg, k
	977	REAL h_oro(klon)
	978	REAL hlim(klon)
	979	REAL, SAVE :: kap=0.4,kapb=0.4
	980	REAL zq
	981	REAL sq(klon), sqz(klon)
	982	REAL, ALLOCATABLE, SAVE :: l0(:)
	983	!$OMP THREADPRIVATE(l0)
	984	REAL fl, zzz, zl0, zq2, zn2
	985	REAL famorti, zzzz, zh_oro, zhlim
	986	REAL l1(klon, klev+1), l2(klon,klev+1)
	987	REAL winds(klon, klev)
	988	REAL xcell
	989	REAL zstdslope(klon)
	990	REAL lmax
	991	REAL l2strat, l2neutre, extent
	992	REAL l2limit(klon)
	993	!===============================================================
	994	! Fonctions utiles
	995	!===============================================================
	996
	997	! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
	998	!..........................................................................
	999
	1000	fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)kapzlev(ig, &
	1001	k)/(kapzlev(ig,k)+l0(ig)),0.5sqrt(q2(ig,k))/sqrt( &
	1002	max(n2(ig,k),1.E-10))), 1.E-5)
	1003
	1004	! Fonction d'amortissement de la turbulence au dessus de la montagne
	1005	! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
	1006	!.....................................................................
	1007
	1008	famorti(zzzz, zh_oro, zhlim)=(-1.)ATAN((zzzz-zh_oro)/(zhlim-zh_oro))2./3.1416+1.
	1009
[2574]	1010	IF (ngrid==0) RETURN
[2561]	1011
	1012	IF (firstcall) THEN
	1013	ALLOCATE (l0(klon))
	1014	firstcall = .FALSE.
	1015	END IF
	1016
	1017
	1018	!=====================================================================
	1019	! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
	1020	!=====================================================================
	1021
	1022
	1023	IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
	1024
	1025
	1026	! Iterative computation of l0
	1027	! This version is kept for iflag_pbl only for convergence
	1028	! with NPv3.1 Cmip5 simulations
	1029	!...................................................................
	1030
	1031	DO ig = 1, ngrid
	1032	sq(ig) = 1.E-10
	1033	sqz(ig) = 1.E-10
	1034	END DO
	1035	DO k = 2, klev - 1
	1036	DO ig = 1, ngrid
	1037	zq = sqrt(q2(ig,k))
	1038	sqz(ig) = sqz(ig) + zqzlev(ig, k)(zlay(ig,k)-zlay(ig,k-1))
	1039	sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
	1040	END DO
	1041	END DO
	1042	DO ig = 1, ngrid
	1043	l0(ig) = 0.2*sqz(ig)/sq(ig)
	1044	END DO
	1045	DO k = 2, klev
	1046	DO ig = 1, ngrid
	1047	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
	1048	END DO
	1049	END DO
	1050
	1051	ELSE
	1052
	1053
	1054	! In all other case, the assymptotic mixing length l0 is imposed (150m)
	1055	!......................................................................
	1056
[2574]	1057	l0(1:ngrid) = 150.
[2561]	1058	DO k = 2, klev
	1059	DO ig = 1, ngrid
	1060	l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
	1061	END DO
	1062	END DO
	1063
	1064	END IF
	1065
	1066	!=================================================================================
	1067	! CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
	1068	! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
	1069	! glacier, la glace de mer et les oc??ans)
	1070	!=================================================================================
	1071
[2574]	1072	l2(1:ngrid,:)=0.0
	1073	l_mixmin(1:ngrid,:,nsrf)=0.
	1074	l_mix(1:ngrid,:,nsrf)=0.
[2561]	1075
	1076	IF (nsrf .EQ. 1) THEN
	1077
	1078	! coefficients
	1079	!--------------
	1080
[2574]	1081	extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1.
	1082	lmax=150. ! Longueur de m??lange max dans l'absolu
[2561]	1083
	1084	! calculs
	1085	!---------
	1086
[2574]	1087	DO ig=1,ngrid
[2561]	1088
	1089	! On calcule la hauteur du relief
	1090	!.................................
	1091	! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
	1092	! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
	1093	! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
	1094	! (en gros, une maille de taille xcell avec une pente constante zstdslope)
	1095	jg=ni(ig)
	1096	! IF (zsig(jg) .EQ. 0.) THEN
	1097	! zstdslope(ig)=0.
	1098	! ELSE
	1099	! xcell=sqrt(cell_area(jg))
	1100	! zstdslope(ig)=max((xcellzsig(jg)-zmea(jg))3 /(3.zsig(jg)),0.)
	1101	! zstdslope(ig)=sqrt(zstdslope(ig))
	1102	! END IF
	1103
	1104	! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief
	1105	h_oro(ig)=zstd(jg)
	1106	hlim(ig)=extent*h_oro(ig)
[2574]	1107	ENDDO
[2561]	1108
[2574]	1109	l2limit(1:ngrid)=0.
[2561]	1110
[2574]	1111	DO k=2,klev
	1112	DO ig=1,ngrid
	1113	winds(ig,k)=sqrt(u(ig,k)2+v(ig,k)2)
	1114	IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie
	1115	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)
	1116	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
	1117	l2neutre=MIN(l2neutre,lmax) ! On majore par lmax
	1118	l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
	1119	l2(ig,k)=l2limit(ig)
[2561]	1120
[2574]	1121	ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
[2561]	1122
	1123	! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
	1124	! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
	1125	! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
[2574]	1126	l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
	1127	ELSE ! Au dessus de extent*h, on prend l2=l0
	1128	l2(ig,k)=0.
	1129	END IF
	1130	ENDDO
	1131	ENDDO
	1132	ENDIF ! pbl_lmixmin_alpha
[2561]	1133
	1134	!==================================================================================
	1135	! On prend le max entre la longueur de melange de blackadar et celle calcul??e
	1136	! en fonction de la topographie
	1137	!===================================================================================
	1138
	1139
	1140	DO k=2,klev
	1141	DO ig=1,ngrid
[2574]	1142	lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
	1143	ENDDO
	1144	ENDDO
[2561]	1145
[2574]	1146	! Diagnostics
[2561]	1147
	1148	DO k=2,klev
[2574]	1149	DO ig=1,ngrid
	1150	jg=ni(ig)
	1151	l_mix(jg,k,nsrf)=lmix(ig,k)
	1152	l_mixmin(jg,k,nsrf)=l2(ig,k)
	1153	ENDDO
[2561]	1154	ENDDO
[2574]	1155	DO ig=1,ngrid
	1156	jg=ni(ig)
	1157	l_mix(jg,1,nsrf)=hlim(ig)
	1158	ENDDO
[2561]	1159
	1160
	1161
	1162	END SUBROUTINE mixinglength

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

Download in other formats:

Original Format