/[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 57 by guez, Mon Jan 30 12:54:02 2012 UTC trunk/Sources/phylmd/coefkz.f revision 233 by guez, Tue Nov 7 10:52:46 2017 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, rugos, u, v, t, &
8         t, q, qsurf, pcfm, pcfh)         q, qsurf, coefm, coefh, ycdragm, ycdragh)
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
12      ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les      ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
13      ! coefficients d'échange turbulent dans l'atmosphère.      ! coefficients d'échange turbulent dans l'atmosphère.
14    
15      USE indicesol, ONLY : is_oce      use clcdrag_m, only: clcdrag
16      USE dimphy, ONLY : klev, klon, max      USE conf_phys_m, ONLY: iflag_pbl
17      USE conf_gcm_m, ONLY : prt_level      USE dimphy, ONLY: klev, klon
18      USE suphec_m, ONLY : rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt      USE fcttre, ONLY: foede, foeew
19      USE yoethf_m, ONLY : r2es, r5ies, r5les, rvtmp2      USE indicesol, ONLY: is_oce
20      USE fcttre, ONLY : dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
21      USE conf_phys_m, ONLY : iflag_pbl      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
   
     ! 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(:, :) ! (klon, 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(:, :) ! (klon, 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):: rugos(:) ! (klon) longeur de rugosite (en m)
34      REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind      REAL, intent(in):: u(:, :), v(:, :) ! (klon, klev) wind
35      REAL, intent(in):: t(klon, klev) ! temperature (K)      REAL, intent(in):: t(:, :) ! (klon, klev) temperature (K)
36      real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)      real, intent(in):: q(:, :) ! (klon, klev) vapeur d'eau (kg/kg)
37      real, intent(in):: qsurf(klon)      real, intent(in):: qsurf(:) ! (knon)
38      REAL, intent(inout):: pcfm(klon, klev) ! coefficient, vitesse      REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse
39      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité  
40        real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
41        ! coefficient, chaleur et humidité
42    
43        real, intent(out):: ycdragm(:), ycdragh(:) ! (knon)
44    
45      ! Local:      ! Local:
46    
47      INTEGER itop(klon) ! numero de couche du sommet de la couche limite      INTEGER knon ! nombre de points a traiter
48        INTEGER itop(size(coefm, 1)) ! (knon) numero de couche du sommet
49                                     ! de la couche limite
50    
51      ! Quelques constantes et options:      ! Quelques constantes et options:
52    
# Line 53  contains Line 56  contains
56      REAL, PARAMETER:: cc = 5.      REAL, PARAMETER:: cc = 5.
57      REAL, PARAMETER:: cd = 5.      REAL, PARAMETER:: cd = 5.
58      REAL, PARAMETER:: clam = 160.      REAL, PARAMETER:: clam = 160.
59      REAL, PARAMETER:: ratqs = 0.05 ! ! largeur de distribution de vapeur d'eau      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
60      LOGICAL, PARAMETER:: richum = .TRUE. ! utilise le nombre de Richardson humide  
61        LOGICAL, PARAMETER:: richum = .TRUE.
62        ! utilise le nombre de Richardson humide
63    
64      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
65      REAL, PARAMETER:: prandtl = 0.4      REAL, PARAMETER:: prandtl = 0.4
66    
67      REAL kstable ! diffusion minimale (situation stable)      REAL kstable ! diffusion minimale (situation stable)
68      REAL, PARAMETER:: mixlen = 35. ! constante controlant longueur de melange      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
69      INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
70    
71      LOGICAL, PARAMETER:: tvirtu = .TRUE.      LOGICAL, PARAMETER:: tvirtu = .TRUE.
72      ! calculer Ri d'une maniere plus performante      ! calculer Ri d'une maniere plus performante
# Line 68  contains Line 74  contains
74      LOGICAL, PARAMETER:: opt_ec = .FALSE.      LOGICAL, PARAMETER:: opt_ec = .FALSE.
75      ! formule du Centre Europeen dans l'atmosphere      ! formule du Centre Europeen dans l'atmosphere
76    
77      INTEGER i, k, kk      INTEGER i, k
78      REAL zgeop(klon, klev)      REAL zgeop(klon, klev)
79      REAL zmgeom(klon)      REAL zmgeom(klon)
80      REAL zri(klon)      REAL ri(klon)
81      REAL zl2(klon)      REAL l2(klon)
82    
83      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
     REAL pcfm1(klon), pcfh1(klon)  
84    
85      REAL zdphi, zdu2, ztvd, ztvu, zcdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
86      REAL zscf      REAL scf
87      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
88      REAL z2geomf, zalh2, zalm2, zscfh, zscfm      logical zdelta
89      REAL, PARAMETER:: t_coup = 273.15      REAL z2geomf, zalh2, alm2, zscfh, scfm
90      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
91    
92      !--------------------------------------------------------------------      !--------------------------------------------------------------------
93    
94      ! Initialiser les sorties      knon = size(coefm, 1)
     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  
95    
96      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
97      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
# Line 138  contains Line 134  contains
134         z1(i) = zgeop(i, 1)         z1(i) = zgeop(i, 1)
135      ENDDO      ENDDO
136    
137      CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &      CALL clcdrag(nsrf, u1, v1, t1, q1, z1, ts, qsurf, rugos, ycdragm, ycdragh)
          rugos, pcfm1, pcfh1)  
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
138    
139      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
140    
141      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
142    
143      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
144         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
# Line 163  contains Line 151  contains
151    
152            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
153    
154            IF (thermcep) THEN            zdelta = RTT >=zt
155               zdelta = MAX(0., SIGN(1., RTT-zt))            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
156               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &                 / (1. + RVTMP2*zq)
157                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
158               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = MIN(0.5, zqs)
159               zqs = MIN(0.5, zqs)            zcor = 1./(1.-RETV*zqs)
160               zcor = 1./(1.-RETV*zqs)            zqs = zqs*zcor
161               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  
162    
163            ! calculer la fraction nuageuse (processus humide):            ! calculer la fraction nuageuse (processus humide):
164    
# Line 199  contains Line 177  contains
177                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1.+RVTMP2*zq) &
178                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
179                    )*(1.+RETV*q(i, k-1))                    )*(1.+RETV*q(i, k-1))
180               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
181               zri(i) = zri(i) &               ri(i) = ri(i) &
182                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
183                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
184                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2*0.5*(ztvd+ztvu))
185            ELSE            ELSE
186               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
187               zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
188                    -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) &
189                    *(pplay(i, k)-pplay(i, k-1)) &                    *(pplay(i, k)-pplay(i, k-1)) &
190                    )*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)))
191               zri(i) = zri(i) + &               ri(i) = ri(i) + &
192                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &
193                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
194                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))
# Line 218  contains Line 196  contains
196    
197            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
198    
199            zcdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
200    
201            IF (opt_ec) THEN            IF (opt_ec) THEN
202               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k-1)+zgeop(i, k)
203               zalm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5*ckap/RG*z2geomf &
204                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1.+0.5*ckap/rg/clam*z2geomf))**2
205               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5*ckap/rg*z2geomf &
206                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
207               IF (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
208                  ! situation instable                  ! situation instable
209                  zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
210                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
211                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(-ri(i)*scf)
212                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
213                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
214                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
215                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
216               ELSE               ELSE
217                  ! situation stable                  ! situation stable
218                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1.+cd*ri(i))
219                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
220                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
221               ENDIF               ENDIF
222            ELSE            ELSE
223               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)) &
224                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
225               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
226               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k)= l2(i) * coefm(i, k)
227               pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
228            ENDIF            ENDIF
229         ENDDO loop_horiz         ENDDO loop_horiz
230      ENDDO loop_vertical      ENDDO loop_vertical
231    
232      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
233        forall (i = 1: knon)
234      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
235         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
236            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
237    
238    END SUBROUTINE coefkz    END SUBROUTINE coefkz
239    
Legend:
Removed from v.57  
changed lines
  Added in v.233
  ViewVC Help
Powered by ViewVC 1.1.21