/[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.f90 revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC trunk/Sources/phylmd/coefkz.f revision 207 by guez, Thu Sep 1 10:30:53 2016 UTC
# Line 5  module coefkz_m Line 5  module coefkz_m
5  contains  contains
6    
7    SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &    SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &
8         t, q, qsurf, pcfm, pcfh)         t, q, qsurf, coefm, coefh)
9    
10      ! Authors: F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS)      ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
11      ! date: 1993/09/22      ! date: 1993/09/22
12      ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les      ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
13      ! coefficients d'échange turbulent dans l'atmosphère.      ! coefficients d'échange turbulent dans l'atmosphère.
14    
15      USE indicesol, ONLY : is_oce      USE indicesol, ONLY: is_oce
16      USE dimphy, ONLY : klev, klon, max      USE dimphy, ONLY: klev, klon
17      USE conf_gcm_m, ONLY : prt_level      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
18      USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
19      USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2      USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats
20      USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep      USE conf_phys_m, ONLY: iflag_pbl
21      USE conf_phys_m, ONLY : iflag_pbl      use clcdrag_m, only: clcdrag
22    
23      ! Arguments:      ! Arguments:
24    
# Line 38  contains Line 38  contains
38      REAL, intent(in):: t(klon, klev) ! temperature (K)      REAL, intent(in):: t(klon, klev) ! temperature (K)
39      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
40      real, intent(in):: qsurf(klon)      real, intent(in):: qsurf(klon)
41      REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse      REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
42      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité  
43        real, intent(out):: coefh(:, :) ! (knon, klev)
44        ! coefficient, chaleur et humidité
45    
46      ! Local:      ! Local:
47    
48      INTEGER itop(klon) ! numero de couche du sommet de la couche limite      INTEGER itop(knon) ! numero de couche du sommet de la couche limite
49    
50      ! Quelques constantes et options:      ! Quelques constantes et options:
51    
# Line 53  contains Line 55  contains
55      REAL, PARAMETER:: cc = 5.      REAL, PARAMETER:: cc = 5.
56      REAL, PARAMETER:: cd = 5.      REAL, PARAMETER:: cd = 5.
57      REAL, PARAMETER:: clam = 160.      REAL, PARAMETER:: clam = 160.
58      REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
59      LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide  
60        LOGICAL, PARAMETER:: richum = .TRUE.
61        ! utilise le nombre de Richardson humide
62    
63      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
64      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
65    
66      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
67      REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
68      INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
69    
70      LOGICAL, PARAMETER:: tvirtu = .TRUE.      LOGICAL, PARAMETER:: tvirtu = .TRUE.
71      ! calculer Ri d'une maniere plus performante      ! calculer Ri d'une maniere plus performante
# Line 68  contains Line 73  contains
73      LOGICAL, PARAMETER:: opt_ec = .FALSE.      LOGICAL, PARAMETER:: opt_ec = .FALSE.
74      ! formule du Centre Europeen dans l'atmosphere      ! formule du Centre Europeen dans l'atmosphere
75    
76      INTEGER i, k, kk      INTEGER i, k
77      REAL zgeop(klon, klev)      REAL zgeop(klon, klev)
78      REAL zmgeom(klon)      REAL zmgeom(klon)
79      REAL zri(klon)      REAL ri(klon)
80      REAL zl2(klon)      REAL l2(klon)
81    
82      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
     REAL pcfm1(klon), pcfh1(klon)  
83    
84      REAL zdphi, zdu2, ztvd, ztvu, zcdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
85      REAL zscf      REAL scf
86      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
87      REAL z2geomf, zalh2, zalm2, zscfh, zscfm      logical zdelta
88      REAL, PARAMETER:: t_coup = 273.15      REAL z2geomf, zalh2, alm2, zscfh, scfm
89      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
90    
91      !--------------------------------------------------------------------      !--------------------------------------------------------------------
92    
     ! Initialiser les sorties  
     DO k = 1, klev  
        DO i = 1, knon  
           pcfm(i, k) = 0.  
           pcfh(i, k) = 0.  
        ENDDO  
     ENDDO  
     DO i = 1, knon  
        itop(i) = 0  
     ENDDO  
   
93      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
94      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
95         DO k = 3, klev         DO k = 3, klev
# Line 139  contains Line 132  contains
132      ENDDO      ENDDO
133    
134      CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &      CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &
135           rugos, pcfm1, pcfh1)           rugos, coefm(:, 1), coefh(:, 1))
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
136    
137      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
138    
139      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
140    
141      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
142         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
# Line 163  contains Line 149  contains
149    
150            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
151    
152            IF (thermcep) THEN            zdelta = RTT >=zt
153               zdelta = MAX(0., SIGN(1., RTT-zt))            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
154               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &                 / (1. + RVTMP2*zq)
155                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
156               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = MIN(0.5, zqs)
157               zqs = MIN(0.5, zqs)            zcor = 1./(1.-RETV*zqs)
158               zcor = 1./(1.-RETV*zqs)            zqs = zqs*zcor
159               zqs = zqs*zcor            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
              zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)  
           ELSE  
              IF (zt .LT. t_coup) THEN  
                 zqs = qsats(zt) / pplay(i, k)  
                 zdqs = dqsats(zt, zqs)  
              ELSE  
                 zqs = qsatl(zt) / pplay(i, k)  
                 zdqs = dqsatl(zt, zqs)  
              ENDIF  
           ENDIF  
160    
161            ! calculer la fraction nuageuse (processus humide):            ! calculer la fraction nuageuse (processus humide):
162    
# Line 199  contains Line 175  contains
175                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1.+RVTMP2*zq) &
176                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
177                    )*(1.+RETV*q(i, k-1))                    )*(1.+RETV*q(i, k-1))
178               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
179               zri(i) = zri(i) &               ri(i) = ri(i) &
180                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
181                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
182                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2*0.5*(ztvd+ztvu))
183            ELSE            ELSE
184               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
185               zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
186                    -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) &                    -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) &
187                    *(pplay(i, k)-pplay(i, k-1)) &                    *(pplay(i, k)-pplay(i, k-1)) &
188                    )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k)))                    )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k)))
189               zri(i) = zri(i) + &               ri(i) = ri(i) + &
190                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &
191                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
192                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))
# Line 218  contains Line 194  contains
194    
195            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
196    
197            zcdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
198    
199            IF (opt_ec) THEN            IF (opt_ec) THEN
200               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k-1)+zgeop(i, k)
201               zalm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5*ckap/RG*z2geomf &
202                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1.+0.5*ckap/rg/clam*z2geomf))**2
203               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5*ckap/rg*z2geomf &
204                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
205               IF (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
206                  ! situation instable                  ! situation instable
207                  zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
208                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
209                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(-ri(i)*scf)
210                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
211                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
212                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
213                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
214               ELSE               ELSE
215                  ! situation stable                  ! situation stable
216                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1.+cd*ri(i))
217                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
218                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
219               ENDIF               ENDIF
220            ELSE            ELSE
221               zl2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &               l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
222                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
223               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
224               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k)= l2(i) * coefm(i, k)
225               pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
226            ENDIF            ENDIF
227         ENDDO loop_horiz         ENDDO loop_horiz
228      ENDDO loop_vertical      ENDDO loop_vertical
229    
230      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
231        forall (i = 1: knon)
232      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
233         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
234            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
235    
236    END SUBROUTINE coefkz    END SUBROUTINE coefkz
237    
Legend:
Removed from v.57  
changed lines
  Added in v.207
  ViewVC Help
Powered by ViewVC 1.1.21