Context Navigation

coef_diff_turb_mod.F90 @ 5160

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

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

Download in other formats:

Original Format