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

Diff of /trunk/phylmd/Interface_surf/coefkz.f

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

revision 287 by guez, Tue Jul 24 15:22:48 2018 UTC revision 288 by guez, Tue Jul 24 16:27:12 2018 UTC
# Line 48  contains Line 48  contains
48      ! Quelques constantes et options:      ! Quelques constantes et options:
49    
50      REAL, PARAMETER:: cepdu2 =0.1**2      REAL, PARAMETER:: cepdu2 =0.1**2
     REAL, PARAMETER:: CKAP = 0.4  
     REAL, PARAMETER:: cb = 5.  
     REAL, PARAMETER:: cc = 5.  
     REAL, PARAMETER:: cd = 5.  
     REAL, PARAMETER:: clam = 160.  
51      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
   
     LOGICAL, PARAMETER:: richum = .TRUE.  
     ! utilise le nombre de Richardson humide  
   
52      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
53      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
54    
55      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
56      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
57      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
   
     LOGICAL, PARAMETER:: tvirtu = .TRUE.  
     ! calculer Ri d'une maniere plus performante  
   
     LOGICAL, PARAMETER:: opt_ec = .FALSE.  
     ! formule du Centre Europeen dans l'atmosphere  
   
58      INTEGER i, k      INTEGER i, k
59      REAL zmgeom(size(ts))      REAL zmgeom(size(ts))
60      REAL ri(size(ts))      REAL ri(size(ts))
61      REAL l2(size(ts))      REAL l2(size(ts))
62      REAL zdphi, zdu2, ztvd, ztvu, cdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
     REAL scf  
63      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
64      logical zdelta      logical zdelta
     REAL z2geomf, zalh2, alm2, zscfh, scfm  
65      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
66    
67      !--------------------------------------------------------------------      !--------------------------------------------------------------------
# Line 98  contains Line 80  contains
80         ENDDO         ENDDO
81      ENDIF      ENDIF
82    
83      IF (nsrf /= is_oce) THEN      kstable = merge(ksta, ksta_ter, nsrf == is_oce)
        kstable = ksta_ter  
     ELSE  
        kstable = ksta  
     ENDIF  
84    
85      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
86    
87      itop = isommet      itop = isommet
88    
89      loop_vertical: DO k = 2, isommet      DO k = 2, isommet
90         loop_horiz: DO i = 1, knon         DO i = 1, knon
91            zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &            zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
92                 + (v(i, k) - v(i, k - 1))**2)                 + (v(i, k) - v(i, k - 1))**2)
93            zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)            zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
# Line 118  contains Line 96  contains
96            zq = (q(i, k) + q(i, k - 1)) * 0.5            zq = (q(i, k) + q(i, k - 1)) * 0.5
97    
98            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
   
99            zdelta = RTT >=zt            zdelta = RTT >=zt
100            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
101                 / (1. + RVTMP2 * zq)                 / (1. + RVTMP2 * zq)
# Line 129  contains Line 106  contains
106            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
107    
108            ! calculer la fraction nuageuse (processus humide):            ! calculer la fraction nuageuse (processus humide):
   
109            zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)            zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
110            zfr = MAX(0.0, MIN(1.0, zfr))            zfr = MAX(0.0, MIN(1.0, zfr))
           IF (.NOT.richum) zfr = 0.0  
111    
112            ! calculer le nombre de Richardson:            ! calculer le nombre de Richardson:
113              ztvd = (t(i, k) &
114            IF (tvirtu) THEN                 + zdphi/RCPD/(1. + RVTMP2 * zq) &
115               ztvd = (t(i, k) &                 * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
116                    + zdphi/RCPD/(1. + RVTMP2 * zq) &                 ) * (1. + RETV * q(i, k))
117                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &            ztvu = (t(i, k - 1) &
118                    ) * (1. + RETV * q(i, k))                 - zdphi/RCPD/(1. + RVTMP2 * zq) &
119               ztvu = (t(i, k - 1) &                 * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
120                    - zdphi/RCPD/(1. + RVTMP2 * zq) &                 ) * (1. + RETV * q(i, k - 1))
121                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &            ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
122                    ) * (1. + RETV * q(i, k - 1))            ri(i) = ri(i) &
123               ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))                 + zmgeom(i) * zmgeom(i)/RG * gamt(k) &
124               ri(i) = ri(i) &                 * (paprs(i, k)/101325.0)**RKAPPA &
125                    + zmgeom(i) * zmgeom(i)/RG * gamt(k) &                 /(zdu2 * 0.5 * (ztvd + ztvu))
                   * (paprs(i, k)/101325.0)**RKAPPA &  
                   /(zdu2 * 0.5 * (ztvd + ztvu))  
           ELSE  
              ! calcul de Ridchardson compatible LMD5  
              ri(i) = (RCPD * (t(i, k) - t(i, k - 1)) &  
                   - RD * 0.5 * (t(i, k) + t(i, k - 1))/paprs(i, k) &  
                   * (pplay(i, k) - pplay(i, k - 1)) &  
                   ) * zmgeom(i)/(zdu2 * 0.5 * RCPD * (t(i, k - 1) + t(i, k)))  
              ri(i) = ri(i) + &  
                   zmgeom(i) * zmgeom(i) * gamt(k)/RG &  
                   * (paprs(i, k)/101325.0)**RKAPPA &  
                   /(zdu2 * 0.5 * (t(i, k - 1) + t(i, k)))  
           ENDIF  
126    
127            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
128    
129            cdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
130    
131            IF (opt_ec) THEN            l2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, itop(i) + 1)) &
132               z2geomf = zgeop(i, k - 1) + zgeop(i, k)                 /(paprs(i, 2) - paprs(i, itop(i) + 1))))**2
133               alm2 = (0.5 * ckap/RG * z2geomf &            coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
134                    /(1. + 0.5 * ckap/rg/clam * z2geomf))**2            coefm(i, k) = l2(i) * coefm(i, k)
135               zalh2 = (0.5 * ckap/rg * z2geomf &            coefh(i, k) = coefm(i, k) / prandtl ! h et m different
136                    /(1. + 0.5 * ckap/RG/(clam * SQRT(1.5 * cd)) * z2geomf))**2         ENDDO
137               IF (ri(i) < 0.) THEN      ENDDO
                 ! situation instable  
                 scf = ((zgeop(i, k)/zgeop(i, k - 1))**(1./3.) - 1.)**3 &  
                      / (zmgeom(i)/RG)**3 / (zgeop(i, k - 1)/RG)  
                 scf = SQRT(- ri(i) * scf)  
                 scfm = 1.0 / (1.0 + 3.0 * cb * cc * alm2 * scf)  
                 zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * scf)  
                 coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)  
                 coefh(i, k) = cdn * zalh2 * (1. - 3.0 * cb * ri(i) * zscfh)  
              ELSE  
                 ! situation stable  
                 scf = SQRT(1. + cd * ri(i))  
                 coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)  
                 coefh(i, k) = cdn * zalh2/(1. + 3.0 * cb * ri(i) * scf)  
              ENDIF  
           ELSE  
              l2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, itop(i) + 1)) &  
                   /(paprs(i, 2) - paprs(i, itop(i) + 1))))**2  
              coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))  
              coefm(i, k) = l2(i) * coefm(i, k)  
              coefh(i, k) = coefm(i, k) / prandtl ! h et m different  
           ENDIF  
        ENDDO loop_horiz  
     ENDDO loop_vertical  
138    
139      ! Au-delà du sommet, pas de diffusion turbulente :      ! Au-delà du sommet, pas de diffusion turbulente :
140      forall (i = 1: knon)      forall (i = 1: knon)
Legend:
Removed from v.287  
changed lines
  Added in v.288
  ViewVC Help
Powered by ViewVC 1.1.21