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

Diff of /trunk/Sources/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 42 by guez, Thu Mar 24 11:52:41 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, u1, v1, t1, q1, z1, ts, qsurf, rugos, &
187  c         pcfm1, pcfh1)
188        DO i = 1, knon  
189         u1(i) = u(i,1)    DO i = 1, knon
190         v1(i) = v(i,1)       pcfm(i, 1)=pcfm1(i)
191         t1(i) = t(i,1)       pcfh(i, 1)=pcfh1(i)
192         q1(i) = q(i,1)    ENDDO
193         z1(i) = zgeop(i,1)  
194        ENDDO    ! Calculer les coefficients turbulents dans l'atmosphere
195  c  
196        CALL clcdrag(klon, knon, nsrf, zxli,    DO i = 1, knon
197       $             u1, v1, t1, q1, z1,       itop(i) = isommet
198       $             ts, qsurf, rugos,    ENDDO
199       $             pcfm1, pcfh1)  
200  cIM  $             ts, qsurf, rugos,    DO k = 2, isommet
201  C       DO i = 1, knon
202        DO i = 1, knon          zdu2=MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &
203         pcfm(i,1)=pcfm1(i)               +(v(i, k)-v(i, k-1))**2)
204         pcfh(i,1)=pcfh1(i)          zmgeom(i)=zgeop(i, k)-zgeop(i, k-1)
205        ENDDO          zdphi =zmgeom(i) / 2.0
206  c          zt = (t(i, k)+t(i, k-1)) * 0.5
207  c Calculer les coefficients turbulents dans l'atmosphere          zq = (q(i, k)+q(i, k-1)) * 0.5
208  c  
209        DO i = 1, knon          !           calculer Qs et dQs/dT:
210           itop(i) = isommet  
211        ENDDO          IF (thermcep) THEN
212               zdelta = MAX(0., SIGN(1., RTT-zt))
213               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &
214        DO k = 2, isommet                  + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
215        DO i = 1, knon             zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
216              zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2             zqs = MIN(0.5, zqs)
217       .                     +(v(i,k)-v(i,k-1))**2)             zcor = 1./(1.-RETV*zqs)
218              zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)             zqs = zqs*zcor
219              zdphi =zmgeom(i) / 2.0             zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
220              zt = (t(i,k)+t(i,k-1)) * 0.5          ELSE
221              zq = (q(i,k)+q(i,k-1)) * 0.5             IF (zt .LT. t_coup) THEN
222  c                zqs = qsats(zt) / pplay(i, k)
223  c           calculer Qs et dQs/dT:                zdqs = dqsats(zt, zqs)
224  c             ELSE
225              IF (thermcep) THEN                zqs = qsatl(zt) / pplay(i, k)
226                zdelta = MAX(0.,SIGN(1.,RTT-zt))                zdqs = dqsatl(zt, zqs)
227                zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)             ENDIF
228       .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta          ENDIF
229                zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)  
230                zqs = MIN(0.5,zqs)          !           calculer la fraction nuageuse (processus humide):
231                zcor = 1./(1.-RETV*zqs)  
232                zqs = zqs*zcor          zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
233                zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)          zfr = MAX(0.0, MIN(1.0, zfr))
234              ELSE          IF (.NOT.richum) zfr = 0.0
235                IF (zt .LT. t_coup) THEN  
236                   zqs = qsats(zt) / pplay(i,k)          !           calculer le nombre de Richardson:
237                   zdqs = dqsats(zt,zqs)  
238                ELSE          IF (tvirtu) THEN
239                   zqs = qsatl(zt) / pplay(i,k)             ztvd =( t(i, k) &
240                   zdqs = dqsatl(zt,zqs)                  + zdphi/RCPD/(1.+RVTMP2*zq) &
241                ENDIF                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
242              ENDIF                  )*(1.+RETV*q(i, k))
243  c             ztvu =( t(i, k-1) &
244  c           calculer la fraction nuageuse (processus humide):                  - zdphi/RCPD/(1.+RVTMP2*zq) &
245  c                  *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
246              zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)                  )*(1.+RETV*q(i, k-1))
247              zfr = MAX(0.0,MIN(1.0,zfr))             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
248              IF (.NOT.richum) zfr = 0.0             zri(i) = zri(i) &
249  c                  + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
250  c           calculer le nombre de Richardson:                  *(paprs(i, k)/101325.0)**RKAPPA &
251  c                  /(zdu2*0.5*(ztvd+ztvu))
252              IF (tvirtu) THEN  
253              ztvd =( t(i,k)          ELSE ! calcul de Ridchardson compatible LMD5
254       .             + zdphi/RCPD/(1.+RVTMP2*zq)  
255       .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )             zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
256       .            )*(1.+RETV*q(i,k))                  -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) &
257              ztvu =( t(i,k-1)                  *(pplay(i, k)-pplay(i, k-1)) &
258       .             - zdphi/RCPD/(1.+RVTMP2*zq)                  )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k)))
259       .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )             zri(i) = zri(i) + &
260       .            )*(1.+RETV*q(i,k-1))                  zmgeom(i)*zmgeom(i)*gamt(k)/RG &
261              zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))                  *(paprs(i, k)/101325.0)**RKAPPA &
262              zri(i) = zri(i)                  /(zdu2*0.5*(t(i, k-1)+t(i, k)))
263       .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)          ENDIF
264       .               *(paprs(i,k)/101325.0)**RKAPPA  
265       .               /(zdu2*0.5*(ztvd+ztvu))          !           finalement, les coefficients d'echange sont obtenus:
266  c  
267              ELSE ! calcul de Ridchardson compatible LMD5          zcdn=SQRT(zdu2) / zmgeom(i) * RG
268  c  
269              zri(i) =(RCPD*(t(i,k)-t(i,k-1))          IF (opt_ec) THEN
270       .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)             z2geomf=zgeop(i, k-1)+zgeop(i, k)
271       .               *(pplay(i,k)-pplay(i,k-1))             zalm2=(0.5*ckap/RG*z2geomf &
272       .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))                  /(1.+0.5*ckap/rg/clam*z2geomf))**2
273              zri(i) = zri(i) +             zalh2=(0.5*ckap/rg*z2geomf &
274       .             zmgeom(i)*zmgeom(i)*gamt(k)/RG                  /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
275  cSB     .             /(paprs(i,k)/101325.0)**RKAPPA             IF (zri(i).LT.0.0) THEN  ! situation instable
276       .             *(paprs(i,k)/101325.0)**RKAPPA                zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
277       .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))                     / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
278              ENDIF                zscf = SQRT(-zri(i)*zscf)
279  c                zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
280  c           finalement, les coefficients d'echange sont obtenus:                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
281  c                pcfm(i, k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
282              zcdn=SQRT(zdu2) / zmgeom(i) * RG                pcfh(i, k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
283  c             ELSE ! situation stable
284            IF (opt_ec) THEN                zscf=SQRT(1.+cd*zri(i))
285              z2geomf=zgeop(i,k-1)+zgeop(i,k)                pcfm(i, k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
286              zalm2=(0.5*ckap/RG*z2geomf                pcfh(i, k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
287       .             /(1.+0.5*ckap/rg/clam*z2geomf))**2             ENDIF
288              zalh2=(0.5*ckap/rg*z2geomf          ELSE
289       .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2             zl2(i)=(mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
290              IF (zri(i).LT.0.0) THEN  ! situation instable                  /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
291                 zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3             pcfm(i, k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
292       .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)             pcfm(i, k)= zl2(i)* pcfm(i, k)
293                 zscf = SQRT(-zri(i)*zscf)             pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different
294                 zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)          ENDIF
295                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)       ENDDO
296                 pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)    ENDDO
297                 pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)  
298              ELSE ! situation stable    ! Au-dela du sommet, pas de diffusion turbulente:
299                 zscf=SQRT(1.+cd*zri(i))  
300                 pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)    DO i = 1, knon
301                 pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)       IF (itop(i)+1 .LE. klev) THEN
302              ENDIF          DO k = itop(i)+1, klev
303            ELSE             pcfh(i, k) = 0.0
304              zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))             pcfm(i, k) = 0.0
305       .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2          ENDDO
306              pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))       ENDIF
307              pcfm(i,k)= zl2(i)* pcfm(i,k)    ENDDO
308              pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different  
309            ENDIF  END SUBROUTINE coefkz
       ENDDO  
       ENDDO  
 c  
 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.42
  ViewVC Help
Powered by ViewVC 1.1.21