/[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/phylmd/coefkz.f revision 103 by guez, Fri Aug 29 13:00:05 2014 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 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, zcvm5, zcor, zqs, zfr, zdqs
88      REAL z2geomf, zalh2, zalm2, zscfh, zscfm      logical zdelta
89        REAL z2geomf, zalh2, alm2, zscfh, scfm
90      REAL, PARAMETER:: t_coup = 273.15      REAL, PARAMETER:: t_coup = 273.15
91      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
92    
93      !--------------------------------------------------------------------      !--------------------------------------------------------------------
94    
     ! 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  
   
95      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
96      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
97         DO k = 3, klev         DO k = 3, klev
# Line 139  contains Line 134  contains
134      ENDDO      ENDDO
135    
136      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, &
137           rugos, pcfm1, pcfh1)           rugos, coefm(:, 1), coefh(:, 1))
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
138    
139      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
140    
141      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
142    
143      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
144         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
# Line 164  contains Line 152  contains
152            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
153    
154            IF (thermcep) THEN            IF (thermcep) THEN
155               zdelta = MAX(0., SIGN(1., RTT-zt))               zdelta = RTT >=zt
156               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &               zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
157                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta                    / (1. + RVTMP2*zq)
158               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
159               zqs = MIN(0.5, zqs)               zqs = MIN(0.5, zqs)
160               zcor = 1./(1.-RETV*zqs)               zcor = 1./(1.-RETV*zqs)
161               zqs = zqs*zcor               zqs = zqs*zcor
162               zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)               zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
163            ELSE            ELSE
164               IF (zt .LT. t_coup) THEN               IF (zt  <  t_coup) THEN
165                  zqs = qsats(zt) / pplay(i, k)                  zqs = qsats(zt) / pplay(i, k)
166                  zdqs = dqsats(zt, zqs)                  zdqs = dqsats(zt, zqs)
167               ELSE               ELSE
# Line 199  contains Line 187  contains
187                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1.+RVTMP2*zq) &
188                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
189                    )*(1.+RETV*q(i, k-1))                    )*(1.+RETV*q(i, k-1))
190               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
191               zri(i) = zri(i) &               ri(i) = ri(i) &
192                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
193                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
194                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2*0.5*(ztvd+ztvu))
195            ELSE            ELSE
196               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
197               zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
198                    -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) &
199                    *(pplay(i, k)-pplay(i, k-1)) &                    *(pplay(i, k)-pplay(i, k-1)) &
200                    )*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)))
201               zri(i) = zri(i) + &               ri(i) = ri(i) + &
202                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &
203                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
204                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))
# Line 218  contains Line 206  contains
206    
207            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
208    
209            zcdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
210    
211            IF (opt_ec) THEN            IF (opt_ec) THEN
212               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k-1)+zgeop(i, k)
213               zalm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5*ckap/RG*z2geomf &
214                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1.+0.5*ckap/rg/clam*z2geomf))**2
215               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5*ckap/rg*z2geomf &
216                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
217               IF (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
218                  ! situation instable                  ! situation instable
219                  zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
220                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
221                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(-ri(i)*scf)
222                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
223                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
224                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
225                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
226               ELSE               ELSE
227                  ! situation stable                  ! situation stable
228                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1.+cd*ri(i))
229                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
230                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
231               ENDIF               ENDIF
232            ELSE            ELSE
233               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)) &
234                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
235               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
236               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k)= l2(i) * coefm(i, k)
237               pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
238            ENDIF            ENDIF
239         ENDDO loop_horiz         ENDDO loop_horiz
240      ENDDO loop_vertical      ENDDO loop_vertical
241    
242      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
243        forall (i = 1: knon)
244      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
245         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
246            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
247    
248    END SUBROUTINE coefkz    END SUBROUTINE coefkz
249    
Legend:
Removed from v.57  
changed lines
  Added in v.103
  ViewVC Help
Powered by ViewVC 1.1.21