/[lmdze]/trunk/phylmd/coefkz.f
ViewVC logotype

Diff of /trunk/phylmd/coefkz.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/coefkz.f revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/libf/phylmd/coefkz.f90 revision 40 by guez, Tue Feb 22 13:49:36 2011 UTC
# Line 1  Line 1 
1        SUBROUTINE coefkz(nsrf, knon, paprs, pplay,  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &
2  cIM 261103       t, q, qsurf, pcfm, pcfh)
3       .                  ksta, ksta_ter,  
4  cIM 261103    ! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 1993/09/22
5       .                  ts, rugos,    !           (une version strictement identique a l'ancien modele)
6       .                  u,v,t,q,    ! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
7       .                  qsurf,    !        coefficients d'echange turbulent dans l'atmosphere.
8       .                  pcfm, pcfh)  
9        use dimens_m    USE indicesol, ONLY : is_oce
10        use indicesol    USE dimphy, ONLY : klev, klon, max
11        use dimphy    USE iniprint, ONLY : prt_level
12        use iniprint    USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
13        use SUPHEC_M    USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2
14        use yoethf_m    USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
15        use fcttre    USE conf_phys_m, ONLY : iflag_pbl
16        use conf_phys_m  
17        IMPLICIT none    IMPLICIT none
18  c======================================================================  
19  c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922    ! Arguments:
20  c           (une version strictement identique a l'ancien modele)    ! nsrf-----input-I- indicateur de la nature du sol
21  c Objet: calculer le coefficient du frottement du sol (Cdrag) et les    ! knon-----input-I- nombre de points a traiter
22  c        coefficients d'echange turbulent dans l'atmosphere.    ! paprs----input-R- pression a chaque intercouche (en Pa)
23  c Arguments:    ! pplay----input-R- pression au milieu de chaque couche (en Pa)
24  c nsrf-----input-I- indicateur de la nature du sol    ! ts-------input-R- temperature du sol (en Kelvin)
25  c knon-----input-I- nombre de points a traiter    ! rugos----input-R- longeur de rugosite (en m)
26  c paprs----input-R- pression a chaque intercouche (en Pa)    ! u--------input-R- vitesse u
27  c pplay----input-R- pression au milieu de chaque couche (en Pa)    ! v--------input-R- vitesse v
28  c ts-------input-R- temperature du sol (en Kelvin)    ! q--------input-R- vapeur d'eau (kg/kg)
29  c rugos----input-R- longeur de rugosite (en m)  
30  c u--------input-R- vitesse u    ! itop-----output-I- numero de couche du sommet de la couche limite
31  c v--------input-R- vitesse v    ! pcfm-----output-R- coefficients a calculer (vitesse)
32  c t--------input-R- temperature (K)    ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
33  c q--------input-R- vapeur d'eau (kg/kg)  
34  c    INTEGER knon, nsrf
35  c itop-----output-I- numero de couche du sommet de la couche limite    REAL ts(klon)
36  c pcfm-----output-R- coefficients a calculer (vitesse)    REAL paprs(klon, klev+1), pplay(klon, klev)
37  c pcfh-----output-R- coefficients a calculer (chaleur et humidite)    REAL u(klon, klev), v(klon, klev), q(klon, klev)
38  c======================================================================    REAL, intent(in):: t(klon, klev) ! temperature (K)
39  c    REAL rugos(klon)
40  c Arguments:  
41  c    REAL pcfm(klon, klev), pcfh(klon, klev)
42        INTEGER knon, nsrf    INTEGER itop(klon)
43        REAL ts(klon)  
44        REAL paprs(klon,klev+1), pplay(klon,klev)    ! Quelques constantes et options:
45        REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)  
46        REAL rugos(klon)    REAL cepdu2, ckap, cb, cc, cd, clam
47  c    PARAMETER (cepdu2 =(0.1)**2)
48        REAL pcfm(klon,klev), pcfh(klon,klev)    PARAMETER (CKAP=0.4)
49        INTEGER itop(klon)    PARAMETER (cb=5.0)
50  c    PARAMETER (cc=5.0)
51  c Quelques constantes et options:    PARAMETER (cd=5.0)
52  c    PARAMETER (clam=160.0)
53        REAL cepdu2, ckap, cb, cc, cd, clam    REAL ratqs ! largeur de distribution de vapeur d'eau
54        PARAMETER (cepdu2 =(0.1)**2)    PARAMETER (ratqs=0.05)
55        PARAMETER (CKAP=0.4)    LOGICAL richum ! utilise le nombre de Richardson humide
56        PARAMETER (cb=5.0)    PARAMETER (richum=.TRUE.)
57        PARAMETER (cc=5.0)    REAL ric ! nombre de Richardson critique
58        PARAMETER (cd=5.0)    PARAMETER(ric=0.4)
59        PARAMETER (clam=160.0)    REAL prandtl
60        REAL ratqs ! largeur de distribution de vapeur d'eau    PARAMETER (prandtl=0.4)
61        PARAMETER (ratqs=0.05)    REAL kstable ! diffusion minimale (situation stable)
62        LOGICAL richum ! utilise le nombre de Richardson humide    ! GKtest
63        PARAMETER (richum=.TRUE.)    ! PARAMETER (kstable=1.0e-10)
64        REAL ric ! nombre de Richardson critique    REAL ksta, ksta_ter
65        PARAMETER(ric=0.4)    !IM: 261103     REAL kstable_ter, kstable_sinon
66        REAL prandtl    !IM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
67        PARAMETER (prandtl=0.4)    !IM: 261103     PARAMETER (kstable_ter = 1.0e-8)
68        REAL kstable ! diffusion minimale (situation stable)    !IM: 261103   PARAMETER (kstable_ter = 1.0e-10)
69        ! GKtest    !IM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
70        ! PARAMETER (kstable=1.0e-10)    ! fin GKtest
71        REAL ksta, ksta_ter    REAL mixlen ! constante controlant longueur de melange
72  cIM: 261103     REAL kstable_ter, kstable_sinon    PARAMETER (mixlen=35.0)
73  cIM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)    INTEGER isommet ! le sommet de la couche limite
74  cIM: 261103     PARAMETER (kstable_ter = 1.0e-8)    PARAMETER (isommet=klev)
75  cIM: 261103   PARAMETER (kstable_ter = 1.0e-10)    LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
76  cIM: 261103   PARAMETER (kstable_sinon = 1.0e-10)    PARAMETER (tvirtu=.TRUE.)
77        ! fin GKtest    LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
78        REAL mixlen ! constante controlant longueur de melange    PARAMETER (opt_ec=.FALSE.)
79        PARAMETER (mixlen=35.0)  
80        INTEGER isommet ! le sommet de la couche limite    ! Variables locales:
81        PARAMETER (isommet=klev)  
82        LOGICAL tvirtu ! calculer Ri d'une maniere plus performante    INTEGER i, k, kk !IM 120704
83        PARAMETER (tvirtu=.TRUE.)    REAL zgeop(klon, klev)
84        LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere    REAL zmgeom(klon)
85        PARAMETER (opt_ec=.FALSE.)    REAL zri(klon)
86      REAL zl2(klon)
87  c  
88  c Variables locales:    REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
89  c    REAL pcfm1(klon), pcfh1(klon)
90        INTEGER i, k, kk !IM 120704  
91        REAL zgeop(klon,klev)    REAL zdphi, zdu2, ztvd, ztvu, zcdn
92        REAL zmgeom(klon)    REAL zscf
93        REAL zri(klon)    REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
94        REAL zl2(klon)    REAL z2geomf, zalh2, zalm2, zscfh, zscfm
95      REAL t_coup
96        REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)    PARAMETER (t_coup=273.15)
97        REAL pcfm1(klon), pcfh1(klon)    !IM
98  c    LOGICAL check
99        REAL zdphi, zdu2, ztvd, ztvu, zcdn    PARAMETER (check=.false.)
100        REAL zscf  
101        REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs    ! contre-gradient pour la chaleur sensible: Kelvin/metre
102        REAL z2geomf, zalh2, zalm2, zscfh, zscfm    REAL gamt(2:klev)
103        REAL t_coup    real qsurf(klon)
104        PARAMETER (t_coup=273.15)  
105  cIM    LOGICAL appel1er
106        LOGICAL check    SAVE appel1er
107        PARAMETER (check=.false.)  
108  c    ! Fonctions thermodynamiques et fonctions d'instabilite
109  c contre-gradient pour la chaleur sensible: Kelvin/metre    REAL fsta, fins, x
110        REAL gamt(2:klev)    LOGICAL zxli ! utiliser un jeu de fonctions simples
111        real qsurf(klon)    PARAMETER (zxli=.FALSE.)
112  c  
113        LOGICAL appel1er    fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
114        SAVE appel1er    fins(x) = SQRT(1.0-18.0*x)
115  c  
116  c Fonctions thermodynamiques et fonctions d'instabilite    DATA appel1er /.TRUE./
117        REAL fsta, fins, x  
118        LOGICAL zxli ! utiliser un jeu de fonctions simples    !--------------------------------------------------------------------
119        PARAMETER (zxli=.FALSE.)  
120  c    IF (appel1er) THEN
121        fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))       if (prt_level > 9) THEN
122        fins(x) = SQRT(1.0-18.0*x)          print *, 'coefkz, opt_ec:', opt_ec
123  c          print *, 'coefkz, richum:', richum
124        DATA appel1er /.TRUE./          IF (richum) print *, 'coefkz, ratqs:', ratqs
125  c          print *, 'coefkz, isommet:', isommet
126        IF (appel1er) THEN          print *, 'coefkz, tvirtu:', tvirtu
127          if (prt_level > 9) THEN          appel1er = .FALSE.
128            print *,'coefkz, opt_ec:', opt_ec       endif
129            print *,'coefkz, richum:', richum    ENDIF
130            IF (richum) print *,'coefkz, ratqs:', ratqs  
131            print *,'coefkz, isommet:', isommet    ! Initialiser les sorties
132            print *,'coefkz, tvirtu:', tvirtu  
133            appel1er = .FALSE.    DO k = 1, klev
134          endif       DO i = 1, knon
135        ENDIF          pcfm(i, k) = 0.0
136  c          pcfh(i, k) = 0.0
137  c Initialiser les sorties       ENDDO
138  c    ENDDO
139        DO k = 1, klev    DO i = 1, knon
140        DO i = 1, knon       itop(i) = 0
141           pcfm(i,k) = 0.0    ENDDO
142           pcfh(i,k) = 0.0  
143        ENDDO    ! Prescrire la valeur de contre-gradient
144        ENDDO  
145        DO i = 1, knon    if (iflag_pbl.eq.1) then
146           itop(i) = 0       DO k = 3, klev
147        ENDDO          gamt(k) = -1.0E-03
148         ENDDO
149  c       gamt(2) = -2.5E-03
150  c Prescrire la valeur de contre-gradient    else
151  c       DO k = 2, klev
152        if (iflag_pbl.eq.1) then          gamt(k) = 0.0
153           DO k = 3, klev       ENDDO
154              gamt(k) = -1.0E-03    ENDIF
155           ENDDO  
156           gamt(2) = -2.5E-03    IF ( nsrf .NE. is_oce ) THEN
157        else       kstable = ksta_ter
158           DO k = 2, klev    ELSE
159              gamt(k) = 0.0       kstable = ksta
160           ENDDO    ENDIF
161        ENDIF  
162  cIM cf JLD/ GKtest    ! Calculer les géopotentiels de chaque couche
163        IF ( nsrf .NE. is_oce ) THEN  
164  cIM 261103     kstable = kstable_ter    DO i = 1, knon
165          kstable = ksta_ter       zgeop(i, 1) = RD * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) &
166        ELSE            * (paprs(i, 1) - pplay(i, 1))
167  cIM 261103     kstable = kstable_sinon    ENDDO
168          kstable = ksta    DO k = 2, klev
169        ENDIF       DO i = 1, knon
170  cIM cf JLD/ GKtest fin          zgeop(i, k) = zgeop(i, k-1) &
171  c               + RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
172  c Calculer les geopotentiels de chaque couche               * (pplay(i, k-1)-pplay(i, k))
173  c       ENDDO
174        DO i = 1, knon    ENDDO
175           zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))  
176       .                   * (paprs(i,1)-pplay(i,1))    ! Calculer le frottement au sol (Cdrag)
177        ENDDO  
178        DO k = 2, klev    DO i = 1, knon
179        DO i = 1, knon       u1(i) = u(i, 1)
180           zgeop(i,k) = zgeop(i,k-1)       v1(i) = v(i, 1)
181       .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)       t1(i) = t(i, 1)
182       .                   * (pplay(i,k-1)-pplay(i,k))       q1(i) = q(i, 1)
183        ENDDO       z1(i) = zgeop(i, 1)
184        ENDDO    ENDDO
185  c  
186  c Calculer le frottement au sol (Cdrag)    CALL clcdrag(klon, knon, nsrf, zxli,  &
187  c         u1, v1, t1, q1, z1, &
188        DO i = 1, knon         ts, qsurf, rugos, &
189         u1(i) = u(i,1)         pcfm1, pcfh1)
190         v1(i) = v(i,1)    !IM  $             ts, qsurf, rugos,
191         t1(i) = t(i,1)  
192         q1(i) = q(i,1)    DO i = 1, knon
193         z1(i) = zgeop(i,1)       pcfm(i, 1)=pcfm1(i)
194        ENDDO       pcfh(i, 1)=pcfh1(i)
195  c    ENDDO
196        CALL clcdrag(klon, knon, nsrf, zxli,  
197       $             u1, v1, t1, q1, z1,    ! Calculer les coefficients turbulents dans l'atmosphere
198       $             ts, qsurf, rugos,  
199       $             pcfm1, pcfh1)    DO i = 1, knon
200  cIM  $             ts, qsurf, rugos,       itop(i) = isommet
201  C    ENDDO
202        DO i = 1, knon  
203         pcfm(i,1)=pcfm1(i)    DO k = 2, isommet
204         pcfh(i,1)=pcfh1(i)       DO i = 1, knon
205        ENDDO          zdu2=MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &
206  c               +(v(i, k)-v(i, k-1))**2)
207  c Calculer les coefficients turbulents dans l'atmosphere          zmgeom(i)=zgeop(i, k)-zgeop(i, k-1)
208  c          zdphi =zmgeom(i) / 2.0
209        DO i = 1, knon          zt = (t(i, k)+t(i, k-1)) * 0.5
210           itop(i) = isommet          zq = (q(i, k)+q(i, k-1)) * 0.5
211        ENDDO  
212            !           calculer Qs et dQs/dT:
213    
214        DO k = 2, isommet          IF (thermcep) THEN
215        DO i = 1, knon             zdelta = MAX(0., SIGN(1., RTT-zt))
216              zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2             zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &
217       .                     +(v(i,k)-v(i,k-1))**2)                  + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
218              zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)             zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
219              zdphi =zmgeom(i) / 2.0             zqs = MIN(0.5, zqs)
220              zt = (t(i,k)+t(i,k-1)) * 0.5             zcor = 1./(1.-RETV*zqs)
221              zq = (q(i,k)+q(i,k-1)) * 0.5             zqs = zqs*zcor
222  c             zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
223  c           calculer Qs et dQs/dT:          ELSE
224  c             IF (zt .LT. t_coup) THEN
225              IF (thermcep) THEN                zqs = qsats(zt) / pplay(i, k)
226                zdelta = MAX(0.,SIGN(1.,RTT-zt))                zdqs = dqsats(zt, zqs)
227                zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)             ELSE
228       .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta                zqs = qsatl(zt) / pplay(i, k)
229                zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)                zdqs = dqsatl(zt, zqs)
230                zqs = MIN(0.5,zqs)             ENDIF
231                zcor = 1./(1.-RETV*zqs)          ENDIF
232                zqs = zqs*zcor  
233                zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)          !           calculer la fraction nuageuse (processus humide):
234              ELSE  
235                IF (zt .LT. t_coup) THEN          zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
236                   zqs = qsats(zt) / pplay(i,k)          zfr = MAX(0.0, MIN(1.0, zfr))
237                   zdqs = dqsats(zt,zqs)          IF (.NOT.richum) zfr = 0.0
238                ELSE  
239                   zqs = qsatl(zt) / pplay(i,k)          !           calculer le nombre de Richardson:
240                   zdqs = dqsatl(zt,zqs)  
241                ENDIF          IF (tvirtu) THEN
242              ENDIF             ztvd =( t(i, k) &
243  c                  + zdphi/RCPD/(1.+RVTMP2*zq) &
244  c           calculer la fraction nuageuse (processus humide):                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
245  c                  )*(1.+RETV*q(i, k))
246              zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)             ztvu =( t(i, k-1) &
247              zfr = MAX(0.0,MIN(1.0,zfr))                  - zdphi/RCPD/(1.+RVTMP2*zq) &
248              IF (.NOT.richum) zfr = 0.0                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
249  c                  )*(1.+RETV*q(i, k-1))
250  c           calculer le nombre de Richardson:             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
251  c             zri(i) = zri(i) &
252              IF (tvirtu) THEN                  + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
253              ztvd =( t(i,k)                  *(paprs(i, k)/101325.0)**RKAPPA &
254       .             + zdphi/RCPD/(1.+RVTMP2*zq)                  /(zdu2*0.5*(ztvd+ztvu))
255       .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )  
256       .            )*(1.+RETV*q(i,k))          ELSE ! calcul de Ridchardson compatible LMD5
257              ztvu =( t(i,k-1)  
258       .             - zdphi/RCPD/(1.+RVTMP2*zq)             zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
259       .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )                  -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) &
260       .            )*(1.+RETV*q(i,k-1))                  *(pplay(i, k)-pplay(i, k-1)) &
261              zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))                  )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k)))
262              zri(i) = zri(i)             zri(i) = zri(i) + &
263       .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)                  zmgeom(i)*zmgeom(i)*gamt(k)/RG &
264       .               *(paprs(i,k)/101325.0)**RKAPPA                  *(paprs(i, k)/101325.0)**RKAPPA &
265       .               /(zdu2*0.5*(ztvd+ztvu))                  /(zdu2*0.5*(t(i, k-1)+t(i, k)))
266  c          ENDIF
267              ELSE ! calcul de Ridchardson compatible LMD5  
268  c          !           finalement, les coefficients d'echange sont obtenus:
269              zri(i) =(RCPD*(t(i,k)-t(i,k-1))  
270       .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)          zcdn=SQRT(zdu2) / zmgeom(i) * RG
271       .               *(pplay(i,k)-pplay(i,k-1))  
272       .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))          IF (opt_ec) THEN
273              zri(i) = zri(i) +             z2geomf=zgeop(i, k-1)+zgeop(i, k)
274       .             zmgeom(i)*zmgeom(i)*gamt(k)/RG             zalm2=(0.5*ckap/RG*z2geomf &
275  cSB     .             /(paprs(i,k)/101325.0)**RKAPPA                  /(1.+0.5*ckap/rg/clam*z2geomf))**2
276       .             *(paprs(i,k)/101325.0)**RKAPPA             zalh2=(0.5*ckap/rg*z2geomf &
277       .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))                  /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
278              ENDIF             IF (zri(i).LT.0.0) THEN  ! situation instable
279  c                zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
280  c           finalement, les coefficients d'echange sont obtenus:                     / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
281  c                zscf = SQRT(-zri(i)*zscf)
282              zcdn=SQRT(zdu2) / zmgeom(i) * RG                zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
283  c                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
284            IF (opt_ec) THEN                pcfm(i, k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
285              z2geomf=zgeop(i,k-1)+zgeop(i,k)                pcfh(i, k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
286              zalm2=(0.5*ckap/RG*z2geomf             ELSE ! situation stable
287       .             /(1.+0.5*ckap/rg/clam*z2geomf))**2                zscf=SQRT(1.+cd*zri(i))
288              zalh2=(0.5*ckap/rg*z2geomf                pcfm(i, k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
289       .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                pcfh(i, k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
290              IF (zri(i).LT.0.0) THEN  ! situation instable             ENDIF
291                 zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3          ELSE
292       .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)             zl2(i)=(mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
293                 zscf = SQRT(-zri(i)*zscf)                  /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
294                 zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)             pcfm(i, k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
295                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)             pcfm(i, k)= zl2(i)* pcfm(i, k)
296                 pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)             pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different
297                 pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)          ENDIF
298              ELSE ! situation stable       ENDDO
299                 zscf=SQRT(1.+cd*zri(i))    ENDDO
300                 pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)  
301                 pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)    ! Au-dela du sommet, pas de diffusion turbulente:
302              ENDIF  
303            ELSE    DO i = 1, knon
304              zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))       IF (itop(i)+1 .LE. klev) THEN
305       .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2          DO k = itop(i)+1, klev
306              pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))             pcfh(i, k) = 0.0
307              pcfm(i,k)= zl2(i)* pcfm(i,k)             pcfm(i, k) = 0.0
308              pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different          ENDDO
309            ENDIF       ENDIF
310        ENDDO    ENDDO
311        ENDDO  
312  c  END SUBROUTINE coefkz
 c Au-dela du sommet, pas de diffusion turbulente:  
 c  
       DO i = 1, knon  
          IF (itop(i)+1 .LE. klev) THEN  
             DO k = itop(i)+1, klev  
                pcfh(i,k) = 0.0  
                pcfm(i,k) = 0.0  
             ENDDO  
          ENDIF  
       ENDDO  
 c  
       RETURN  
       END  
Legend:
Removed from v.38  
changed lines
  Added in v.40
  ViewVC Help
Powered by ViewVC 1.1.21