Context Navigation

coef_diff_turb_mod.F90 @ 3333

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

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

Download in other formats:

Original Format