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