/[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 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, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &    SUBROUTINE coefkz(nsrf, paprs, pplay, ksta, ksta_ter, ts, u, v, t, q, &
8         t, q, qsurf, pcfm, pcfh)         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: 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.  
   
     USE indicesol, ONLY : is_oce  
     USE dimphy, ONLY : klev, klon, max  
     USE iniprint, ONLY : prt_level  
     USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt  
     USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2  
     USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep  
     USE conf_phys_m, ONLY : iflag_pbl  
12    
13      ! Arguments:      ! Objet : calculer les coefficients d'échange turbulent dans
14        ! l'atmosphère.
15    
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
20        USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
21        USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
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):: ksta, ksta_ter
32      REAL, intent(in):: ts(klon) ! 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(klon, klev), 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)
37      real, intent(in):: qsurf(klon)      REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse
38      REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse  
39      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité      real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
40        ! coefficient, chaleur et humidité
41    
42      ! Local:      ! Local:
43    
44      INTEGER itop(klon) ! numero de couche du sommet de la couche limite      INTEGER knon ! nombre de points a traiter
45    
46        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 53  contains Line 54  contains
54      REAL, PARAMETER:: cc = 5.      REAL, PARAMETER:: cc = 5.
55      REAL, PARAMETER:: cd = 5.      REAL, PARAMETER:: cd = 5.
56      REAL, PARAMETER:: clam = 160.      REAL, PARAMETER:: clam = 160.
57      REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
58      LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide  
59        LOGICAL, PARAMETER:: richum = .TRUE.
60        ! utilise le nombre de Richardson humide
61    
62      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
63      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
64    
65      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
66      REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
67      INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
68    
69      LOGICAL, PARAMETER:: tvirtu = .TRUE.      LOGICAL, PARAMETER:: tvirtu = .TRUE.
70      ! calculer Ri d'une maniere plus performante      ! calculer Ri d'une maniere plus performante
# Line 68  contains Line 72  contains
72      LOGICAL, PARAMETER:: opt_ec = .FALSE.      LOGICAL, PARAMETER:: opt_ec = .FALSE.
73      ! formule du Centre Europeen dans l'atmosphere      ! formule du Centre Europeen dans l'atmosphere
74    
75      INTEGER i, k, kk      INTEGER i, k
76      REAL zgeop(klon, klev)      REAL zmgeom(size(ts))
77      REAL zmgeom(klon)      REAL ri(size(ts))
78      REAL zri(klon)      REAL l2(size(ts))
79      REAL zl2(klon)      REAL zdphi, zdu2, ztvd, ztvu, cdn
80        REAL scf
81      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
82      REAL pcfm1(klon), pcfh1(klon)      logical zdelta
83        REAL z2geomf, zalh2, alm2, zscfh, scfm
     REAL zdphi, zdu2, ztvd, ztvu, zcdn  
     REAL zscf  
     REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs  
     REAL z2geomf, zalh2, zalm2, zscfh, zscfm  
     REAL, PARAMETER:: t_coup = 273.15  
84      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
85    
86      !--------------------------------------------------------------------      !--------------------------------------------------------------------
87    
88      ! Initialiser les sorties      knon = size(ts)
     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  
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(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &  
          rugos, pcfm1, pcfh1)  
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
   
108      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
109    
110      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
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            IF (thermcep) THEN            zdelta = RTT >=zt
124               zdelta = MAX(0., SIGN(1., RTT-zt))            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
125               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &                 / (1. + RVTMP2 * zq)
126                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
127               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = MIN(0.5, zqs)
128               zqs = MIN(0.5, zqs)            zcor = 1./(1. - RETV * zqs)
129               zcor = 1./(1.-RETV*zqs)            zqs = zqs * zcor
130               zqs = zqs*zcor            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
              zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)  
           ELSE  
              IF (zt .LT. t_coup) THEN  
                 zqs = qsats(zt) / pplay(i, k)  
                 zdqs = dqsats(zt, zqs)  
              ELSE  
                 zqs = qsatl(zt) / pplay(i, k)  
                 zdqs = dqsatl(zt, zqs)  
              ENDIF  
           ENDIF  
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               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
150               zri(i) = zri(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               zri(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               zri(i) = zri(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:
167    
168            zcdn = 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               zalm2 = (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 (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
177                  ! situation instable                  ! situation instable
178                  zscf = ((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                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(- ri(i) * scf)
181                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0 + 3.0 * cb * cc * alm2 * scf)
182                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * scf)
183                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
184                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn * zalh2 * (1. - 3.0 * cb * ri(i) * zscfh)
185               ELSE               ELSE
186                  ! situation stable                  ! situation stable
187                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1. + cd * ri(i))
188                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
189                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn * zalh2/(1. + 3.0 * cb * ri(i) * scf)
190               ENDIF               ENDIF
191            ELSE            ELSE
192               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)) &
193                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2) - paprs(i, itop(i) + 1)) ))**2
194               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
195               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k) = l2(i) * coefm(i, k)
196               pcfh(i, k) = pcfm(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
199      ENDDO loop_vertical      ENDDO loop_vertical
200    
201      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
202        forall (i = 1: knon)
203      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
204         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
205            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
206    
207    END SUBROUTINE coefkz    END SUBROUTINE coefkz
208    
Legend:
Removed from v.47  
changed lines
  Added in v.248
  ViewVC Help
Powered by ViewVC 1.1.21