/[lmdze]/trunk/libf/phylmd/coefkz.f90
ViewVC logotype

Diff of /trunk/libf/phylmd/coefkz.f90

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

revision 61 by guez, Mon Jan 30 12:54:02 2012 UTC revision 62 by guez, Thu Jul 26 14:37:37 2012 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, max
17      USE conf_gcm_m, ONLY : prt_level      USE conf_gcm_m, ONLY: prt_level
18      USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
19      USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
20      USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
21      USE conf_phys_m, ONLY : iflag_pbl      USE conf_phys_m, ONLY: iflag_pbl
22        use clcdrag_m, only: clcdrag
23    
24      ! Arguments:      ! Arguments:
25    
# Line 38  contains Line 39  contains
39      REAL, intent(in):: t(klon, klev) ! temperature (K)      REAL, intent(in):: t(klon, klev) ! temperature (K)
40      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
41      real, intent(in):: qsurf(klon)      real, intent(in):: qsurf(klon)
42      REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse      REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
43      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité  
44        real, intent(out):: coefh(:, :) ! (knon, klev)
45        ! coefficient, chaleur et humidité
46    
47      ! Local:      ! Local:
48    
49      INTEGER itop(klon) ! numero de couche du sommet de la couche limite      INTEGER itop(knon) ! numero de couche du sommet de la couche limite
50    
51      ! Quelques constantes et options:      ! Quelques constantes et options:
52    
# Line 53  contains Line 56  contains
56      REAL, PARAMETER:: cc = 5.      REAL, PARAMETER:: cc = 5.
57      REAL, PARAMETER:: cd = 5.      REAL, PARAMETER:: cd = 5.
58      REAL, PARAMETER:: clam = 160.      REAL, PARAMETER:: clam = 160.
59      REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
60      LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide      
61        LOGICAL, PARAMETER:: richum = .TRUE.
62        ! utilise le nombre de Richardson humide
63    
64      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
65      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
66    
67      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
68      REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
69      INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
70    
71      LOGICAL, PARAMETER:: tvirtu = .TRUE.      LOGICAL, PARAMETER:: tvirtu = .TRUE.
72      ! calculer Ri d'une maniere plus performante      ! calculer Ri d'une maniere plus performante
# Line 71  contains Line 77  contains
77      INTEGER i, k, kk      INTEGER i, k, kk
78      REAL zgeop(klon, klev)      REAL zgeop(klon, klev)
79      REAL zmgeom(klon)      REAL zmgeom(klon)
80      REAL zri(klon)      REAL ri(klon)
81      REAL zl2(klon)      REAL l2(klon)
82    
83      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)  
84    
85      REAL zdphi, zdu2, ztvd, ztvu, zcdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
86      REAL zscf      REAL scf
87      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
88      REAL z2geomf, zalh2, zalm2, zscfh, zscfm      REAL z2geomf, zalh2, alm2, zscfh, scfm
89      REAL, PARAMETER:: t_coup = 273.15      REAL, PARAMETER:: t_coup = 273.15
90      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
91    
92      !--------------------------------------------------------------------      !--------------------------------------------------------------------
93    
     ! 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  
   
94      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
95      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
96         DO k = 3, klev         DO k = 3, klev
# Line 139  contains Line 133  contains
133      ENDDO      ENDDO
134    
135      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, &
136           rugos, pcfm1, pcfh1)           rugos, coefm(:, 1), coefh(:, 1))
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
137    
138      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
139    
140      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
141    
142      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
143         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
# Line 173  contains Line 160  contains
160               zqs = zqs*zcor               zqs = zqs*zcor
161               zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)               zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
162            ELSE            ELSE
163               IF (zt .LT. t_coup) THEN               IF (zt  <  t_coup) THEN
164                  zqs = qsats(zt) / pplay(i, k)                  zqs = qsats(zt) / pplay(i, k)
165                  zdqs = dqsats(zt, zqs)                  zdqs = dqsats(zt, zqs)
166               ELSE               ELSE
# Line 199  contains Line 186  contains
186                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1.+RVTMP2*zq) &
187                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
188                    )*(1.+RETV*q(i, k-1))                    )*(1.+RETV*q(i, k-1))
189               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
190               zri(i) = zri(i) &               ri(i) = ri(i) &
191                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
192                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
193                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2*0.5*(ztvd+ztvu))
194            ELSE            ELSE
195               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
196               zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
197                    -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) &
198                    *(pplay(i, k)-pplay(i, k-1)) &                    *(pplay(i, k)-pplay(i, k-1)) &
199                    )*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)))
200               zri(i) = zri(i) + &               ri(i) = ri(i) + &
201                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &
202                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
203                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))
# Line 218  contains Line 205  contains
205    
206            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
207    
208            zcdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
209    
210            IF (opt_ec) THEN            IF (opt_ec) THEN
211               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k-1)+zgeop(i, k)
212               zalm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5*ckap/RG*z2geomf &
213                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1.+0.5*ckap/rg/clam*z2geomf))**2
214               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5*ckap/rg*z2geomf &
215                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
216               IF (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
217                  ! situation instable                  ! situation instable
218                  zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
219                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
220                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(-ri(i)*scf)
221                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
222                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
223                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
224                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
225               ELSE               ELSE
226                  ! situation stable                  ! situation stable
227                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1.+cd*ri(i))
228                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
229                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
230               ENDIF               ENDIF
231            ELSE            ELSE
232               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)) &
233                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
234               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
235               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k)= l2(i) * coefm(i, k)
236               pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
237            ENDIF            ENDIF
238         ENDDO loop_horiz         ENDDO loop_horiz
239      ENDDO loop_vertical      ENDDO loop_vertical
240    
241      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
242        forall (i = 1: knon)
243      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
244         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
245            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
246    
247    END SUBROUTINE coefkz    END SUBROUTINE coefkz
248    
Legend:
Removed from v.61  
changed lines
  Added in v.62
  ViewVC Help
Powered by ViewVC 1.1.21