Context Navigation

coef_diff_turb_mod.F90 @ 2031

Last change on this file since 2031 was 1932, checked in by lguez, 11 years ago

Declaration of variables appearing in array bounds must be specified
before the declaration of the array.

Variables used in FCTTRE.h are declared in YOMCST.h and YOETHF.h so
YOMCST.h and YOETHF.h must be included before FCTTRE.h.

(There were compilation errors with gfortran and debugging options.)

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:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 18.8 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
[1785]	19	USE indice_sol_mod
[1067]	20	!
[782]	21	! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
[1067]	22	! atmosphere
	23	! NB! No values are calculated between surface and the first model layer.
	24	! ycoefm(:,1) and ycoefh(:,1) are not valid !!!
[782]	25	!
	26	!
	27	! Input arguments
	28	!****************************************************************************************
	29	REAL, INTENT(IN) :: dtime
	30	INTEGER, INTENT(IN) :: nsrf, knon
	31	INTEGER, DIMENSION(klon), INTENT(IN) :: ni
	32	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs
	33	REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay
	34	REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv
	35	REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt
	36	REAL, DIMENSION(klon), INTENT(IN) :: yts, yrugos, yqsurf
[1067]	37	REAL, DIMENSION(klon), INTENT(IN) :: ycdragm
	38
	39	! InOutput arguments
	40	!****************************************************************************************
	41	REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
[782]	42
	43	! Output arguments
	44	!****************************************************************************************
	45	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh
	46	REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm
	47
	48	! Other local variables
	49	!****************************************************************************************
	50	INTEGER :: k, i, j
	51	REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta
[878]	52	REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq
[782]	53	REAL, DIMENSION(klon) :: yustar
	54
	55	! Include
	56	!****************************************************************************************
[793]	57	INCLUDE "clesphys.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
[1761]	160	IF (iflag_pbl.GE.31) THEN
[782]	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)
[1761]	165	ELSE IF (iflag_pbl<20) THEN
[782]	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
[1785]	189	USE indice_sol_mod
[1067]	190
[782]	191	!======================================================================
	192	! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
	193	! (une version strictement identique a l'ancien modele)
	194	! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
	195	! coefficients d'echange turbulent dans l'atmosphere.
	196	! Arguments:
	197	! nsrf-----input-I- indicateur de la nature du sol
	198	! knon-----input-I- nombre de points a traiter
	199	! paprs----input-R- pregssion a chaque intercouche (en Pa)
	200	! pplay----input-R- pression au milieu de chaque couche (en Pa)
	201	! ts-------input-R- temperature du sol (en Kelvin)
	202	! rugos----input-R- longeur de rugosite (en m)
	203	! u--------input-R- vitesse u
	204	! v--------input-R- vitesse v
	205	! t--------input-R- temperature (K)
	206	! q--------input-R- vapeur d'eau (kg/kg)
	207	!
	208	! pcfm-----output-R- coefficients a calculer (vitesse)
	209	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	210	!======================================================================
[793]	211	INCLUDE "YOETHF.h"
[1932]	212	INCLUDE "YOMCST.h"
[793]	213	INCLUDE "FCTTRE.h"
[782]	214	INCLUDE "iniprint.h"
	215	INCLUDE "compbl.h"
	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	!
[1604]	391	if (zq /= 0.) then
	392	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
	393	else
	394	zfr = 0.
	395	end if
[782]	396	zfr = MAX(0.0,MIN(1.0,zfr))
	397	IF (.NOT.richum) zfr = 0.0
	398	!
	399	! calculer le nombre de Richardson:
	400	!
	401	IF (tvirtu) THEN
	402	ztvd =( t(i,k) &
	403	+ zdphi/RCPD/(1.+RVTMP2*zq) &
	404	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
	405	)(1.+RETVq(i,k))
	406	ztvu =( t(i,k-1) &
	407	- zdphi/RCPD/(1.+RVTMP2*zq) &
	408	( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
	409	)(1.+RETVq(i,k-1))
	410	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
	411	zri(i) = zri(i) &
	412	+ zmgeom(i)zmgeom(i)/RGgamt(k) &
	413	(paprs(i,k)/101325.0)*RKAPPA &
	414	/(zdu20.5(ztvd+ztvu))
	415
	416	ELSE ! calcul de Ridchardson compatible LMD5
	417
	418	zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
	419	-RD0.5(t(i,k)+t(i,k-1))/paprs(i,k) &
	420	*(pplay(i,k)-pplay(i,k-1)) &
	421	)zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
	422	zri(i) = zri(i) + &
	423	zmgeom(i)zmgeom(i)gamt(k)/RG &
	424	(paprs(i,k)/101325.0)*RKAPPA &
	425	/(zdu20.5(t(i,k-1)+t(i,k)))
	426	ENDIF
	427	!
	428	! finalement, les coefficients d'echange sont obtenus:
	429	!
	430	zcdn=SQRT(zdu2) / zmgeom(i) * RG
	431
	432	IF (opt_ec) THEN
	433	z2geomf=zgeop(i,k-1)+zgeop(i,k)
	434	zalm2=(0.5ckap/RGz2geomf &
	435	/(1.+0.5ckap/rg/clamz2geomf))**2
	436	zalh2=(0.5ckap/rgz2geomf &
	437	/(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
	438	IF (zri(i).LT.0.0) THEN ! situation instable
	439	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3 &
	440	/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
	441	zscf = SQRT(-zri(i)*zscf)
	442	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
	443	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
	444	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
	445	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
	446	ELSE ! situation stable
	447	zscf=SQRT(1.+cd*zri(i))
	448	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
	449	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
	450	ENDIF
	451	ELSE
	452	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
	453	/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
	454	pcfm(i,k)=SQRT(MAX(zcdnzcdn(ric-zri(i))/ric, kstable))
	455	pcfm(i,k)= zl2(i)* pcfm(i,k)
	456	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	457	ENDIF
	458	ENDDO
	459	ENDDO
	460
	461	!
	462	! Au-dela du sommet, pas de diffusion turbulente:
	463	!
	464	DO i = 1, knon
	465	IF (itop(i)+1 .LE. klev) THEN
	466	DO k = itop(i)+1, klev
	467	pcfh(i,k) = 0.0
	468	pcfm(i,k) = 0.0
	469	ENDDO
	470	ENDIF
	471	ENDDO
	472
	473	END SUBROUTINE coefkz
	474	!
	475	!****************************************************************************************
	476	!
	477	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
	478	pcfm, pcfh)
	479
[1067]	480	USE dimphy
[1785]	481	USE indice_sol_mod
[1067]	482
[782]	483	!======================================================================
	484	! J'introduit un peu de diffusion sauf dans les endroits
	485	! ou une forte inversion est presente
	486	! On peut dire qu'il represente la convection peu profonde
	487	!
	488	! Arguments:
	489	! nsrf-----input-I- indicateur de la nature du sol
	490	! knon-----input-I- nombre de points a traiter
	491	! paprs----input-R- pression a chaque intercouche (en Pa)
	492	! pplay----input-R- pression au milieu de chaque couche (en Pa)
	493	! t--------input-R- temperature (K)
	494	!
	495	! pcfm-----output-R- coefficients a calculer (vitesse)
	496	! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	497	!======================================================================
	498	!
	499	! Arguments:
	500	!
	501	INTEGER, INTENT(IN) :: knon, nsrf
	502	REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
	503	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
	504	REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev)
	505
	506	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
	507	!
	508	! Quelques constantes et options:
	509	!
[1067]	510	REAL, PARAMETER :: prandtl=0.4
	511	REAL, PARAMETER :: kstable=0.002
	512	! REAL, PARAMETER :: kstable=0.001
	513	REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
	514	REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
[782]	515	! PARAMETER (seuil=-0.04)
	516	! PARAMETER (seuil=-0.06)
	517	! PARAMETER (seuil=-0.09)
	518
	519	!
	520	! Variables locales:
	521	!
	522	INTEGER i, k, invb(knon)
	523	REAL zl2(knon)
	524	REAL zdthmin(knon), zdthdp
	525
[793]	526	INCLUDE "YOMCST.h"
[782]	527	!
	528	! Initialiser les sorties
	529	!
	530	DO k = 1, klev
	531	DO i = 1, knon
	532	pcfm(i,k) = 0.0
	533	pcfh(i,k) = 0.0
	534	ENDDO
	535	ENDDO
	536
	537	!
	538	! Chercher la zone d'inversion forte
	539	!
	540	DO i = 1, knon
	541	invb(i) = klev
	542	zdthmin(i)=0.0
	543	ENDDO
	544	DO k = 2, klev/2-1
	545	DO i = 1, knon
	546	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
	547	- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
	548	zdthdp = zdthdp * 100.0
	549	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. &
	550	zdthdp.LT.zdthmin(i) ) THEN
	551	zdthmin(i) = zdthdp
	552	invb(i) = k
	553	ENDIF
	554	ENDDO
	555	ENDDO
	556
	557	!
	558	! Introduire une diffusion:
	559	!
	560	IF ( nsrf.EQ.is_oce ) THEN
	561	DO k = 2, klev
	562	DO i = 1, knon
	563	!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
	564	!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
	565	!IM cf JLD/ GKtest TERkz2
	566	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
	567	! fin GKtest
	568
	569
	570	! s'il n'y a pas d'inversion ou si l'inversion est trop faible
	571	! IF ( (nsrf.EQ.is_oce) .AND. &
	572	IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN
	573	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
	574	/(paprs(i,2)-paprs(i,klev+1)) ))**2
	575	pcfm(i,k)= zl2(i)* kstable
	576	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	577	ENDIF
	578	ENDDO
	579	ENDDO
	580	ENDIF
	581
	582	END SUBROUTINE coefkz2
	583	!
	584	!****************************************************************************************
	585	!
	586	END MODULE coef_diff_turb_mod

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

Download in other formats:

Original Format