/[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 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/Sources/phylmd/coefkz.f revision 178 by guez, Fri Mar 11 18:47:26 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 iniprint, 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, thermcep
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 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 164  contains Line 151  contains
151            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
152    
153            IF (thermcep) THEN            IF (thermcep) THEN
154               zdelta = MAX(0., SIGN(1., RTT-zt))               zdelta = RTT >=zt
155               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &               zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
156                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta                    / (1. + RVTMP2*zq)
157               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
158               zqs = MIN(0.5, zqs)               zqs = MIN(0.5, zqs)
159               zcor = 1./(1.-RETV*zqs)               zcor = 1./(1.-RETV*zqs)
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.47  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21