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

trunk/libf/phylmd/coefkz.f90 revision 47 by guez, Fri Jul 1 15:00:48 2011 UTC trunk/phylmd/coefkz.f revision 254 by guez, Mon Feb 5 10:39:38 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, pcfm, pcfh)  
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.  
   
     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  
11    
12      ! Arguments:      ! 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
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):: 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
37      real, intent(in):: qsurf(klon)  
38      REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse      real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
39      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité      ! coefficient, chaleur et humidité
40    
41      ! Local:      ! Local:
42    
43      INTEGER itop(klon) ! 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 53  contains Line 53  contains
53      REAL, PARAMETER:: cc = 5.      REAL, PARAMETER:: cc = 5.
54      REAL, PARAMETER:: cd = 5.      REAL, PARAMETER:: cd = 5.
55      REAL, PARAMETER:: clam = 160.      REAL, PARAMETER:: clam = 160.
56      REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
57      LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide  
58        LOGICAL, PARAMETER:: richum = .TRUE.
59        ! utilise le nombre de Richardson humide
60    
61      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
62      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
63    
64      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
65      REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
66      INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
67    
68      LOGICAL, PARAMETER:: tvirtu = .TRUE.      LOGICAL, PARAMETER:: tvirtu = .TRUE.
69      ! calculer Ri d'une maniere plus performante      ! calculer Ri d'une maniere plus performante
# Line 68  contains Line 71  contains
71      LOGICAL, PARAMETER:: opt_ec = .FALSE.      LOGICAL, PARAMETER:: opt_ec = .FALSE.
72      ! formule du Centre Europeen dans l'atmosphere      ! formule du Centre Europeen dans l'atmosphere
73    
74      INTEGER i, k, kk      INTEGER i, k
75      REAL zgeop(klon, klev)      REAL zmgeom(size(ts))
76      REAL zmgeom(klon)      REAL ri(size(ts))
77      REAL zri(klon)      REAL l2(size(ts))
78      REAL zl2(klon)      REAL zdphi, zdu2, ztvd, ztvu, cdn
79        REAL scf
80      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
81      REAL pcfm1(klon), pcfh1(klon)      logical zdelta
82        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  
83      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
84    
85      !--------------------------------------------------------------------      !--------------------------------------------------------------------
86    
87      ! 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  
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, pcfm1, pcfh1)  
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
   
107      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
108    
109      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
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            IF (thermcep) THEN            zdelta = RTT >=zt
123               zdelta = MAX(0., SIGN(1., RTT-zt))            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
124               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &                 / (1. + RVTMP2 * zq)
125                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
126               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = MIN(0.5, zqs)
127               zqs = MIN(0.5, zqs)            zcor = 1./(1. - RETV * zqs)
128               zcor = 1./(1.-RETV*zqs)            zqs = zqs * zcor
129               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  
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               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
149               zri(i) = zri(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               zri(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               zri(i) = zri(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:
166    
167            zcdn = 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               zalm2 = (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 (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
176                  ! situation instable                  ! situation instable
177                  zscf = ((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                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(- ri(i) * scf)
180                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0 + 3.0 * cb * cc * alm2 * scf)
181                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * scf)
182                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
183                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn * zalh2 * (1. - 3.0 * cb * ri(i) * zscfh)
184               ELSE               ELSE
185                  ! situation stable                  ! situation stable
186                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1. + cd * ri(i))
187                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
188                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn * zalh2/(1. + 3.0 * cb * ri(i) * scf)
189               ENDIF               ENDIF
190            ELSE            ELSE
191               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)) &
192                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2) - paprs(i, itop(i) + 1)) ))**2
193               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
194               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k) = l2(i) * coefm(i, k)
195               pcfh(i, k) = pcfm(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
198      ENDDO loop_vertical      ENDDO loop_vertical
199    
200      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
201        forall (i = 1: knon)
202      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
203         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
204            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
205    
206    END SUBROUTINE coefkz    END SUBROUTINE coefkz
207    
Legend:
Removed from v.47  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21