Context Navigation

coef_diff_turb_mod.F90 @ 4270

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

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

Download in other formats:

Original Format