/[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

revision 207 by guez, Thu Sep 1 10:30:53 2016 UTC revision 250 by guez, Fri Jan 5 18:18:53 2018 UTC
# Line 4  module coefkz_m Line 4  module coefkz_m
4    
5  contains  contains
6    
7    SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &    SUBROUTINE coefkz(nsrf, paprs, pplay, ts, u, v, t, q, zgeop, coefm, coefh)
        t, q, qsurf, coefm, coefh)  
8    
9      ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)      ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
10      ! date: 1993/09/22      ! Date: September 22nd, 1993
     ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les  
     ! coefficients d'échange turbulent dans l'atmosphère.  
11    
12        ! Objet : calculer les coefficients d'échange turbulent dans
13        ! l'atmosphère.
14    
15        USE clesphys, ONLY: ksta, ksta_ter
16        USE conf_phys_m, ONLY: iflag_pbl
17        USE dimphy, ONLY: klev
18        USE fcttre, ONLY: foede, foeew
19      USE indicesol, ONLY: is_oce      USE indicesol, ONLY: is_oce
     USE dimphy, ONLY: klev, klon  
20      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
21      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
     USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats  
     USE conf_phys_m, ONLY: iflag_pbl  
     use clcdrag_m, only: clcdrag  
   
     ! Arguments:  
22    
23      integer, intent(in):: nsrf ! indicateur de la nature du sol      integer, intent(in):: nsrf ! indicateur de la nature du sol
     INTEGER, intent(in):: knon ! nombre de points a traiter  
24    
25      REAL, intent(in):: paprs(klon, klev+1)      REAL, intent(in):: paprs(:, :) ! (knon, klev + 1)
26      ! pression a chaque intercouche (en Pa)      ! pression a chaque intercouche (en Pa)
27    
28      real, intent(in):: pplay(klon, klev)      real, intent(in):: pplay(:, :) ! (knon, klev)
29      ! pression au milieu de chaque couche (en Pa)      ! pression au milieu de chaque couche (en Pa)
30    
31      REAL, intent(in):: ksta, ksta_ter      REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
32      REAL, intent(in):: ts(klon) ! temperature du sol (en Kelvin)      REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev) wind
33      REAL, intent(in):: rugos(klon) ! longeur de rugosite (en m)      REAL, intent(in):: t(:, :) ! (knon, klev) temperature (K)
34      REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind      real, intent(in):: q(:, :) ! (knon, klev) vapeur d'eau (kg/kg)
35      REAL, intent(in):: t(klon, klev) ! temperature (K)      REAL, intent(in):: zgeop(:, :) ! (knon, klev)
36      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)      REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse
     real, intent(in):: qsurf(klon)  
     REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse  
37    
38      real, intent(out):: coefh(:, :) ! (knon, klev)      real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
39      ! coefficient, chaleur et humidité      ! coefficient, chaleur et humidité
40    
41      ! Local:      ! Local:
42    
43      INTEGER itop(knon) ! numero de couche du sommet de la couche limite      INTEGER knon ! nombre de points a traiter
44    
45        INTEGER itop(size(ts)) ! (knon)
46        ! numero de couche du sommet de la couche limite
47    
48      ! Quelques constantes et options:      ! Quelques constantes et options:
49    
# Line 74  contains Line 72  contains
72      ! formule du Centre Europeen dans l'atmosphere      ! formule du Centre Europeen dans l'atmosphere
73    
74      INTEGER i, k      INTEGER i, k
75      REAL zgeop(klon, klev)      REAL zmgeom(size(ts))
76      REAL zmgeom(klon)      REAL ri(size(ts))
77      REAL ri(klon)      REAL l2(size(ts))
     REAL l2(klon)  
   
     REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)  
   
78      REAL zdphi, zdu2, ztvd, ztvu, cdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
79      REAL scf      REAL scf
80      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
# Line 90  contains Line 84  contains
84    
85      !--------------------------------------------------------------------      !--------------------------------------------------------------------
86    
87        knon = size(ts)
88    
89      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
90      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
91         DO k = 3, klev         DO k = 3, klev
92            gamt(k) = -1.0E-03            gamt(k) = - 1E-3
93         ENDDO         ENDDO
94         gamt(2) = -2.5E-03         gamt(2) = - 2.5E-3
95      else      else
96         DO k = 2, klev         DO k = 2, klev
97            gamt(k) = 0.0            gamt(k) = 0.0
98         ENDDO         ENDDO
99      ENDIF      ENDIF
100    
101      IF ( nsrf .NE. is_oce ) THEN      IF (nsrf .NE. is_oce ) THEN
102         kstable = ksta_ter         kstable = ksta_ter
103      ELSE      ELSE
104         kstable = ksta         kstable = ksta
105      ENDIF      ENDIF
106    
     ! Calculer les géopotentiels de chaque couche  
     DO i = 1, knon  
        zgeop(i, 1) = RD * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) &  
             * (paprs(i, 1) - pplay(i, 1))  
     ENDDO  
     DO k = 2, klev  
        DO i = 1, knon  
           zgeop(i, k) = zgeop(i, k-1) &  
                + RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &  
                * (pplay(i, k-1)-pplay(i, k))  
        ENDDO  
     ENDDO  
   
     ! Calculer le frottement au sol (Cdrag)  
   
     DO i = 1, knon  
        u1(i) = u(i, 1)  
        v1(i) = v(i, 1)  
        t1(i) = t(i, 1)  
        q1(i) = q(i, 1)  
        z1(i) = zgeop(i, 1)  
     ENDDO  
   
     CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &  
          rugos, coefm(:, 1), coefh(:, 1))  
   
107      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
108    
109      itop = isommet      itop = isommet
110    
111      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
112         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
113            zdu2 = MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &            zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
114                 +(v(i, k)-v(i, k-1))**2)                 + (v(i, k) - v(i, k - 1))**2)
115            zmgeom(i) = zgeop(i, k)-zgeop(i, k-1)            zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
116            zdphi =zmgeom(i) / 2.0            zdphi = zmgeom(i) / 2.0
117            zt = (t(i, k)+t(i, k-1)) * 0.5            zt = (t(i, k) + t(i, k - 1)) * 0.5
118            zq = (q(i, k)+q(i, k-1)) * 0.5            zq = (q(i, k) + q(i, k - 1)) * 0.5
119    
120            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
121    
122            zdelta = RTT >=zt            zdelta = RTT >=zt
123            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
124                 / (1. + RVTMP2*zq)                 / (1. + RVTMP2 * zq)
125            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
126            zqs = MIN(0.5, zqs)            zqs = MIN(0.5, zqs)
127            zcor = 1./(1.-RETV*zqs)            zcor = 1./(1. - RETV * zqs)
128            zqs = zqs*zcor            zqs = zqs * zcor
129            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
130    
131            ! calculer la fraction nuageuse (processus humide):            ! calculer la fraction nuageuse (processus humide):
132    
133            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)            zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
134            zfr = MAX(0.0, MIN(1.0, zfr))            zfr = MAX(0.0, MIN(1.0, zfr))
135            IF (.NOT.richum) zfr = 0.0            IF (.NOT.richum) zfr = 0.0
136    
137            !  calculer le nombre de Richardson:            !  calculer le nombre de Richardson:
138    
139            IF (tvirtu) THEN            IF (tvirtu) THEN
140               ztvd =( t(i, k) &               ztvd = (t(i, k) &
141                    + zdphi/RCPD/(1.+RVTMP2*zq) &                    + zdphi/RCPD/(1. + RVTMP2 * zq) &
142                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs) ) &
143                    )*(1.+RETV*q(i, k))                    ) * (1. + RETV * q(i, k))
144               ztvu =( t(i, k-1) &               ztvu = (t(i, k - 1) &
145                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1. + RVTMP2 * zq) &
146                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs) ) &
147                    )*(1.+RETV*q(i, k-1))                    ) * (1. + RETV * q(i, k - 1))
148               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
149               ri(i) = ri(i) &               ri(i) = ri(i) &
150                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i) * zmgeom(i)/RG * gamt(k) &
151                    *(paprs(i, k)/101325.0)**RKAPPA &                    * (paprs(i, k)/101325.0)**RKAPPA &
152                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2 * 0.5 * (ztvd + ztvu))
153            ELSE            ELSE
154               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
155               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) = (RCPD * (t(i, k) - t(i, k - 1)) &
156                    -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) &
157                    *(pplay(i, k)-pplay(i, k-1)) &                    * (pplay(i, k) - pplay(i, k - 1)) &
158                    )*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)))
159               ri(i) = ri(i) + &               ri(i) = ri(i) + &
160                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i) * zmgeom(i) * gamt(k)/RG &
161                    *(paprs(i, k)/101325.0)**RKAPPA &                    * (paprs(i, k)/101325.0)**RKAPPA &
162                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2 * 0.5 * (t(i, k - 1) + t(i, k)))
163            ENDIF            ENDIF
164    
165            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
# Line 197  contains Line 167  contains
167            cdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
168    
169            IF (opt_ec) THEN            IF (opt_ec) THEN
170               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k - 1) + zgeop(i, k)
171               alm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5 * ckap/RG * z2geomf &
172                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1. + 0.5 * ckap/rg/clam * z2geomf))**2
173               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5 * ckap/rg * z2geomf &
174                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1. + 0.5 * ckap/RG/(clam * SQRT(1.5 * cd)) * z2geomf))**2
175               IF (ri(i) < 0.) THEN               IF (ri(i) < 0.) THEN
176                  ! situation instable                  ! situation instable
177                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k - 1))**(1./3.) - 1.)**3 &
178                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k - 1)/RG)
179                  scf = SQRT(-ri(i)*scf)                  scf = SQRT(- ri(i) * scf)
180                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)                  scfm = 1.0 / (1.0 + 3.0 * cb * cc * alm2 * scf)
181                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)                  zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * scf)
182                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
183                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)                  coefh(i, k) = cdn * zalh2 * (1. - 3.0 * cb * ri(i) * zscfh)
184               ELSE               ELSE
185                  ! situation stable                  ! situation stable
186                  scf = SQRT(1.+cd*ri(i))                  scf = SQRT(1. + cd * ri(i))
187                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
188                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)                  coefh(i, k) = cdn * zalh2/(1. + 3.0 * cb * ri(i) * scf)
189               ENDIF               ENDIF
190            ELSE            ELSE
191               l2(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)) &
192                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2) - paprs(i, itop(i) + 1)) ))**2
193               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
194               coefm(i, k)= l2(i) * coefm(i, k)               coefm(i, k) = l2(i) * coefm(i, k)
195               coefh(i, k) = coefm(i, k) / prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
196            ENDIF            ENDIF
197         ENDDO loop_horiz         ENDDO loop_horiz
Legend:
Removed from v.207  
changed lines
  Added in v.250
  ViewVC Help
Powered by ViewVC 1.1.21