Context Navigation

coef_diff_turb_mod.F90 @ 2197

Last change on this file since 2197 was 1067, checked in by Laurent Fairhead, 16 years ago

Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

Preparation pour un couplage des courrant oceaniques.

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 18.7 KB

Rev	Line
[782]	1	!
	2	MODULE coef_diff_turb_mod
	3	!
	4	! This module contains some procedures for calculation of the coefficients of the
	5	! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion
	6	! at surface(cdrag)
	7	!
	8	IMPLICIT NONE
	9
	10	CONTAINS
	11	!
	12	!****************************************************************************************
	13	!
	14	SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
[1067]	15	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
[878]	16	ycoefm, ycoefh ,yq2)
[1067]	17
	18	USE dimphy
	19	!
[782]	20	! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
[1067]	21	! atmosphere
	22	! NB! No values are calculated between surface and the first model layer.
	23	! ycoefm(:,1) and ycoefh(:,1) are not valid !!!
[782]	24	!
	25	!
	26	! Input arguments
	27	!****************************************************************************************
	28	REAL, INTENT(IN) :: dtime
	29	INTEGER, INTENT(IN) :: nsrf, knon
	30	INTEGER, DIMENSION(klon), INTENT(IN) :: ni
	31	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs
	32	REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay
	33	REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv
	34	REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt
	35	REAL, DIMENSION(klon), INTENT(IN) :: yts, yrugos, yqsurf
[1067]	36	REAL, DIMENSION(klon), INTENT(IN) :: ycdragm
	37
	38	! InOutput arguments
	39	!****************************************************************************************
	40	REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
[782]	41
	42	! Output arguments
	43	!****************************************************************************************
	44	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh
	45	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm
	46
	47	! Other local variables
	48	!****************************************************************************************
	49	INTEGER :: k, i, j
	50	REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta
[878]	51	REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq
[782]	52	REAL, DIMENSION(klon) :: yustar
	53
	54	! Include
	55	!****************************************************************************************
[793]	56	INCLUDE "clesphys.h"
	57	INCLUDE "indicesol.h"
[782]	58	INCLUDE "iniprint.h"
	59	INCLUDE "compbl.h"
[793]	60	INCLUDE "YOETHF.h"
	61	INCLUDE "YOMCST.h"
[782]	62
	63
	64	!****************************************************************************************
[1067]	65	! Calcul de coefficients de diffusion turbulent de l'atmosphere :
	66	! ycoefm(:,2:klev), ycoefh(:,2:klev)
[782]	67	!
	68	!****************************************************************************************
	69
	70	CALL coefkz(nsrf, knon, ypaprs, ypplay, &
	71	ksta, ksta_ter, &
	72	yts, yrugos, yu, yv, yt, yq, &
	73	yqsurf, &
	74	ycoefm, ycoefh)
	75
	76	!****************************************************************************************
[1067]	77	! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere :
	78	! ycoefm(:,2:klev), ycoefh(:,2:klev)
[782]	79	!
	80	!****************************************************************************************
	81
	82	IF (iflag_pbl.EQ.1) THEN
	83	CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
	84	ycoefm0, ycoefh0)
	85
[878]	86	DO k = 2, klev
[782]	87	DO i = 1, knon
	88	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	89	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	90	ENDDO
	91	ENDDO
	92	ENDIF
	93
	94
	95	!****************************************************************************************
	96	! Calcul d'une diffusion minimale pour les conditions tres stables
	97	!
	98	!****************************************************************************************
	99	IF (ok_kzmin) THEN
[1067]	100	CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, &
[782]	101	ycoefm0,ycoefh0)
	102
[878]	103	DO k = 2, klev
[782]	104	DO i = 1, knon
	105	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	106	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	107	ENDDO
	108	ENDDO
	109
	110	ENDIF
	111
	112
	113	!****************************************************************************************
	114	! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
[1067]	115	!
[782]	116	!****************************************************************************************
	117
	118	IF (iflag_pbl.GE.3) THEN
	119
	120	yzlay(1:knon,1)= &
	121	RDyt(1:knon,1)/(0.5(ypaprs(1:knon,1)+ypplay(1:knon,1))) &
	122	*(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
	123	DO k=2,klev
	124	DO i = 1, knon
	125	yzlay(i,k)= &
	126	yzlay(i,k-1)+RD0.5(yt(i,k-1)+yt(i,k)) &
	127	/ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
	128	END DO
	129	END DO
	130
	131	DO k=1,klev
	132	DO i = 1, knon
	133	yteta(i,k)= &
	134	yt(i,k)(ypaprs(i,1)/ypplay(i,k))*RKAPPA &
	135	(1.+0.61yq(i,k))
	136	END DO
	137	END DO
	138
	139	yzlev(1:knon,1)=0.
	140	yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
	141	DO k=2,klev
	142	DO i = 1, knon
	143	yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1))
	144	END DO
	145	END DO
	146
[1067]	147	!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	148	!!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
	149	!!$! bug sur les coefficients de surface :
	150	!!$! ycdragh(1:knon) = ycoefm(1:knon,1)
	151	!!$! ycdragm(1:knon) = ycoefh(1:knon,1)
	152	!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	153	CALL ustarhb(knon,yu,yv,ycdragm, yustar)
[782]	154
	155	IF (prt_level > 9) THEN
	156	WRITE(lunout,*) 'USTAR = ',yustar
	157	ENDIF
	158
	159	! iflag_pbl peut etre utilise comme longuer de melange
	160	IF (iflag_pbl.GE.11) THEN
	161	CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
	162	yzlev,yzlay,yu,yv,yteta, &
[1067]	163	ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
[782]	164	iflag_pbl)
	165	ELSE
	166	CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
	167	yzlev,yzlay,yu,yv,yteta, &
[1067]	168	ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &
[782]	169	iflag_pbl)
	170	ENDIF
	171
	172	ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
	173	ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
	174
	175	ENDIF !(iflag_pbl.ge.3)
	176
	177	END SUBROUTINE coef_diff_turb
	178	!
	179	!****************************************************************************************
	180	!
	181	SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
	182	ksta, ksta_ter, &
	183	ts, rugos, &
	184	u,v,t,q, &
	185	qsurf, &
	186	pcfm, pcfh)
	187
[1067]	188	USE dimphy
	189
[782]	190	!======================================================================
	191	! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
	192	! (une version strictement identique a l'ancien modele)
	193	! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
	194	! coefficients d'echange turbulent dans l'atmosphere.
	195	! Arguments:
	196	! nsrf-----input-I- indicateur de la nature du sol
	197	! knon-----input-I- nombre de points a traiter
	198	! paprs----input-R- pregssion a chaque intercouche (en Pa)
	199	! pplay----input-R- pression au milieu de chaque couche (en Pa)
	200	! ts-------input-R- temperature du sol (en Kelvin)
	201	! rugos----input-R- longeur de rugosite (en m)
	202	! u--------input-R- vitesse u
	203	! v--------input-R- vitesse v
	204	! t--------input-R- temperature (K)
	205	! q--------input-R- vapeur d'eau (kg/kg)
	206	!
	207	! pcfm-----output-R- coefficients a calculer (vitesse)
	208	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	209	!======================================================================
[793]	210	INCLUDE "YOETHF.h"
	211	INCLUDE "FCTTRE.h"
[782]	212	INCLUDE "iniprint.h"
[793]	213	INCLUDE "indicesol.h"
[782]	214	INCLUDE "compbl.h"
[793]	215	INCLUDE "YOMCST.h"
[782]	216	!
	217	! Arguments:
	218	!
	219	INTEGER, INTENT(IN) :: knon, nsrf
[1067]	220	REAL, INTENT(IN) :: ksta, ksta_ter
[782]	221	REAL, DIMENSION(klon), INTENT(IN) :: ts
[1067]	222	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
	223	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
[782]	224	REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q
	225	REAL, DIMENSION(klon), INTENT(IN) :: rugos
[1067]	226	REAL, DIMENSION(klon), INTENT(IN) :: qsurf
[782]	227
	228	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh
	229
	230	!
[1067]	231	! Local variables:
	232	!
	233	INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite
	234	!
[782]	235	! Quelques constantes et options:
	236	!
[1067]	237	REAL, PARAMETER :: cepdu2=0.1**2
	238	REAL, PARAMETER :: CKAP=0.4
	239	REAL, PARAMETER :: cb=5.0
	240	REAL, PARAMETER :: cc=5.0
	241	REAL, PARAMETER :: cd=5.0
	242	REAL, PARAMETER :: clam=160.0
	243	REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
	244	LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
	245	REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
	246	REAL, PARAMETER :: prandtl=0.4
[782]	247	REAL kstable ! diffusion minimale (situation stable)
	248	! GKtest
	249	! PARAMETER (kstable=1.0e-10)
	250	!IM: 261103 REAL kstable_ter, kstable_sinon
	251	!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
	252	!IM: 261103 PARAMETER (kstable_ter = 1.0e-8)
	253	!IM: 261103 PARAMETER (kstable_ter = 1.0e-10)
	254	!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
	255	! fin GKtest
[1067]	256	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
[782]	257	INTEGER isommet ! le sommet de la couche limite
[1067]	258	LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
	259	LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
[782]	260
	261	!
	262	! Variables locales:
	263	INTEGER i, k !IM 120704
	264	REAL zgeop(klon,klev)
	265	REAL zmgeom(klon)
	266	REAL zri(klon)
	267	REAL zl2(klon)
	268	REAL zdphi, zdu2, ztvd, ztvu, zcdn
	269	REAL zscf
	270	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
	271	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
[1067]	272	REAL, PARAMETER :: t_coup=273.15
	273	LOGICAL, PARAMETER :: check=.FALSE.
[782]	274	!
	275	! contre-gradient pour la chaleur sensible: Kelvin/metre
	276	REAL gamt(2:klev)
	277
[1067]	278	LOGICAL, SAVE :: appel1er=.TRUE.
[782]	279	!$OMP THREADPRIVATE(appel1er)
	280	!
	281	! Fonctions thermodynamiques et fonctions d'instabilite
	282	REAL fsta, fins, x
	283
	284	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
	285	fins(x) = SQRT(1.0-18.0*x)
	286
	287	isommet=klev
	288
	289	IF (appel1er) THEN
	290	IF (prt_level > 9) THEN
	291	WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
	292	WRITE(lunout,*)'coefkz, richum:', richum
	293	IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
	294	WRITE(lunout,*)'coefkz, isommet:', isommet
	295	WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
	296	appel1er = .FALSE.
	297	ENDIF
	298	ENDIF
	299	!
	300	! Initialiser les sorties
	301	!
	302	DO k = 1, klev
	303	DO i = 1, knon
	304	pcfm(i,k) = 0.0
	305	pcfh(i,k) = 0.0
	306	ENDDO
	307	ENDDO
	308	DO i = 1, knon
	309	itop(i) = 0
	310	ENDDO
	311
	312	!
	313	! Prescrire la valeur de contre-gradient
	314	!
	315	IF (iflag_pbl.EQ.1) THEN
	316	DO k = 3, klev
	317	gamt(k) = -1.0E-03
	318	ENDDO
	319	gamt(2) = -2.5E-03
	320	ELSE
	321	DO k = 2, klev
	322	gamt(k) = 0.0
	323	ENDDO
	324	ENDIF
	325	!IM cf JLD/ GKtest
	326	IF ( nsrf .NE. is_oce ) THEN
	327	!IM 261103 kstable = kstable_ter
	328	kstable = ksta_ter
	329	ELSE
	330	!IM 261103 kstable = kstable_sinon
	331	kstable = ksta
	332	ENDIF
	333	!IM cf JLD/ GKtest fin
	334
	335	!
	336	! Calculer les geopotentiels de chaque couche
	337	!
	338	DO i = 1, knon
	339	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
	340	* (paprs(i,1)-pplay(i,1))
	341	ENDDO
	342	DO k = 2, klev
	343	DO i = 1, knon
	344	zgeop(i,k) = zgeop(i,k-1) &
	345	+ RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
	346	* (pplay(i,k-1)-pplay(i,k))
	347	ENDDO
	348	ENDDO
	349
	350	!
	351	! Calculer les coefficients turbulents dans l'atmosphere
	352	!
	353	DO i = 1, knon
	354	itop(i) = isommet
	355	ENDDO
	356
	357
	358	DO k = 2, isommet
	359	DO i = 1, knon
	360	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
	361	+(v(i,k)-v(i,k-1))**2)
	362	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
	363	zdphi =zmgeom(i) / 2.0
	364	zt = (t(i,k)+t(i,k-1)) * 0.5
	365	zq = (q(i,k)+q(i,k-1)) * 0.5
	366
	367	!
	368	! Calculer Qs et dQs/dT:
	369	!
	370	IF (thermcep) THEN
	371	zdelta = MAX(0.,SIGN(1.,RTT-zt))
	372	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta) &
	373	+ R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
	374	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
	375	zqs = MIN(0.5,zqs)
	376	zcor = 1./(1.-RETV*zqs)
	377	zqs = zqs*zcor
	378	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
	379	ELSE
	380	IF (zt .LT. t_coup) THEN
	381	zqs = qsats(zt) / pplay(i,k)
	382	zdqs = dqsats(zt,zqs)
	383	ELSE
	384	zqs = qsatl(zt) / pplay(i,k)
	385	zdqs = dqsatl(zt,zqs)
	386	ENDIF
	387	ENDIF
	388	!
	389	! calculer la fraction nuageuse (processus humide):
	390	!
	391	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
	392	zfr = MAX(0.0,MIN(1.0,zfr))
	393	IF (.NOT.richum) zfr = 0.0
	394	!
	395	! calculer le nombre de Richardson:
	396	!
	397	IF (tvirtu) THEN
	398	ztvd =( t(i,k) &
	399	+ zdphi/RCPD/(1.+RVTMP2*zq) &
	400	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
	401	)(1.+RETVq(i,k))
	402	ztvu =( t(i,k-1) &
	403	- zdphi/RCPD/(1.+RVTMP2*zq) &
	404	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
	405	)(1.+RETVq(i,k-1))
	406	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
	407	zri(i) = zri(i) &
	408	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
	409	(paprs(i,k)/101325.0)*RKAPPA &
	410	/(zdu20.5(ztvd+ztvu))
	411
	412	ELSE ! calcul de Ridchardson compatible LMD5
	413
	414	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
	415	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
	416	*(pplay(i,k)-pplay(i,k-1)) &
	417	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
	418	zri(i) = zri(i) + &
	419	zmgeom(i)zmgeom(i)gamt(k)/RG &
	420	(paprs(i,k)/101325.0)*RKAPPA &
	421	/(zdu20.5(t(i,k-1)+t(i,k)))
	422	ENDIF
	423	!
	424	! finalement, les coefficients d'echange sont obtenus:
	425	!
	426	zcdn=SQRT(zdu2) / zmgeom(i) * RG
	427
	428	IF (opt_ec) THEN
	429	z2geomf=zgeop(i,k-1)+zgeop(i,k)
	430	zalm2=(0.5ckap/RGz2geomf &
	431	/(1.+0.5ckap/rg/clamz2geomf))**2
	432	zalh2=(0.5ckap/rgz2geomf &
	433	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
	434	IF (zri(i).LT.0.0) THEN ! situation instable
	435	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
	436	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
	437	zscf = SQRT(-zri(i)*zscf)
	438	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
	439	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
	440	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
	441	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
	442	ELSE ! situation stable
	443	zscf=SQRT(1.+cd*zri(i))
	444	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
	445	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
	446	ENDIF
	447	ELSE
	448	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
	449	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
	450	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
	451	pcfm(i,k)= zl2(i)* pcfm(i,k)
	452	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	453	ENDIF
	454	ENDDO
	455	ENDDO
	456
	457	!
	458	! Au-dela du sommet, pas de diffusion turbulente:
	459	!
	460	DO i = 1, knon
	461	IF (itop(i)+1 .LE. klev) THEN
	462	DO k = itop(i)+1, klev
	463	pcfh(i,k) = 0.0
	464	pcfm(i,k) = 0.0
	465	ENDDO
	466	ENDIF
	467	ENDDO
	468
	469	END SUBROUTINE coefkz
	470	!
	471	!****************************************************************************************
	472	!
	473	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
	474	pcfm, pcfh)
	475
[1067]	476	USE dimphy
	477
[782]	478	!======================================================================
	479	! J'introduit un peu de diffusion sauf dans les endroits
	480	! ou une forte inversion est presente
	481	! On peut dire qu'il represente la convection peu profonde
	482	!
	483	! Arguments:
	484	! nsrf-----input-I- indicateur de la nature du sol
	485	! knon-----input-I- nombre de points a traiter
	486	! paprs----input-R- pression a chaque intercouche (en Pa)
	487	! pplay----input-R- pression au milieu de chaque couche (en Pa)
	488	! t--------input-R- temperature (K)
	489	!
	490	! pcfm-----output-R- coefficients a calculer (vitesse)
	491	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	492	!======================================================================
	493	!
	494	! Arguments:
	495	!
	496	INTEGER, INTENT(IN) :: knon, nsrf
	497	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
	498	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
	499	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
	500
	501	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
	502	!
	503	! Quelques constantes et options:
	504	!
[1067]	505	REAL, PARAMETER :: prandtl=0.4
	506	REAL, PARAMETER :: kstable=0.002
	507	! REAL, PARAMETER :: kstable=0.001
	508	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
	509	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
[782]	510	! PARAMETER (seuil=-0.04)
	511	! PARAMETER (seuil=-0.06)
	512	! PARAMETER (seuil=-0.09)
	513
	514	!
	515	! Variables locales:
	516	!
	517	INTEGER i, k, invb(knon)
	518	REAL zl2(knon)
	519	REAL zdthmin(knon), zdthdp
	520
[793]	521	INCLUDE "indicesol.h"
	522	INCLUDE "YOMCST.h"
[782]	523	!
	524	! Initialiser les sorties
	525	!
	526	DO k = 1, klev
	527	DO i = 1, knon
	528	pcfm(i,k) = 0.0
	529	pcfh(i,k) = 0.0
	530	ENDDO
	531	ENDDO
	532
	533	!
	534	! Chercher la zone d'inversion forte
	535	!
	536	DO i = 1, knon
	537	invb(i) = klev
	538	zdthmin(i)=0.0
	539	ENDDO
	540	DO k = 2, klev/2-1
	541	DO i = 1, knon
	542	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
	543	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
	544	zdthdp = zdthdp * 100.0
	545	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
	546	zdthdp.LT.zdthmin(i) ) THEN
	547	zdthmin(i) = zdthdp
	548	invb(i) = k
	549	ENDIF
	550	ENDDO
	551	ENDDO
	552
	553	!
	554	! Introduire une diffusion:
	555	!
	556	IF ( nsrf.EQ.is_oce ) THEN
	557	DO k = 2, klev
	558	DO i = 1, knon
	559	!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
	560	!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
	561	!IM cf JLD/ GKtest TERkz2
	562	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
	563	! fin GKtest
	564
	565
	566	! s'il n'y a pas d'inversion ou si l'inversion est trop faible
	567	! IF ( (nsrf.EQ.is_oce) .AND. &
	568	IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
	569	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
	570	/(paprs(i,2)-paprs(i,klev+1)) ))**2
	571	pcfm(i,k)= zl2(i)* kstable
	572	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	573	ENDIF
	574	ENDDO
	575	ENDDO
	576	ENDIF
	577
	578	END SUBROUTINE coefkz2
	579	!
	580	!****************************************************************************************
	581	!
	582	END MODULE coef_diff_turb_mod

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

Download in other formats:

Original Format