/[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 229 by guez, Mon Nov 6 17:20:45 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)
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(:, :) ! (knon, klev) coefficient, vitesse
39      real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité  
40        real, intent(out):: coefh(:, :) ! (knon, klev)
41        ! coefficient, chaleur et humidité
42    
43      ! Local:      ! Local:
44    
45      INTEGER itop(klon) ! numero de couche du sommet de la couche limite      INTEGER knon ! nombre de points a traiter
46        INTEGER itop(size(coefm, 1)) ! (knon) numero de couche du sommet
47                                     ! 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 zgeop(klon, klev)
77      REAL zmgeom(klon)      REAL zmgeom(klon)
78      REAL zri(klon)      REAL ri(klon)
79      REAL zl2(klon)      REAL l2(klon)
80    
81      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)  
82    
83      REAL zdphi, zdu2, ztvd, ztvu, zcdn      REAL zdphi, zdu2, ztvd, ztvu, cdn
84      REAL zscf      REAL scf
85      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
86      REAL z2geomf, zalh2, zalm2, zscfh, zscfm      logical zdelta
87      REAL, PARAMETER:: t_coup = 273.15      REAL z2geomf, zalh2, alm2, zscfh, scfm
88      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
89    
90      !--------------------------------------------------------------------      !--------------------------------------------------------------------
91    
92      ! 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  
93    
94      ! Prescrire la valeur de contre-gradient      ! Prescrire la valeur de contre-gradient
95      if (iflag_pbl.eq.1) then      if (iflag_pbl.eq.1) then
# Line 138  contains Line 132  contains
132         z1(i) = zgeop(i, 1)         z1(i) = zgeop(i, 1)
133      ENDDO      ENDDO
134    
135      CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &      CALL clcdrag(nsrf, u1, v1, t1, q1, z1, ts, qsurf, rugos, coefm(:, 1), &
136           rugos, pcfm1, pcfh1)           coefh(:, 1))
   
     DO i = 1, knon  
        pcfm(i, 1) = pcfm1(i)  
        pcfh(i, 1) = pcfh1(i)  
     ENDDO  
137    
138      ! Calculer les coefficients turbulents dans l'atmosphere      ! Calculer les coefficients turbulents dans l'atmosphere
139    
140      DO i = 1, knon      itop = isommet
        itop(i) = isommet  
     ENDDO  
141    
142      loop_vertical: DO k = 2, isommet      loop_vertical: DO k = 2, isommet
143         loop_horiz: DO i = 1, knon         loop_horiz: DO i = 1, knon
# Line 163  contains Line 150  contains
150    
151            ! calculer Qs et dQs/dT:            ! calculer Qs et dQs/dT:
152    
153            IF (thermcep) THEN            zdelta = RTT >=zt
154               zdelta = MAX(0., SIGN(1., RTT-zt))            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
155               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &                 / (1. + RVTMP2*zq)
156                    + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
157               zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)            zqs = MIN(0.5, zqs)
158               zqs = MIN(0.5, zqs)            zcor = 1./(1.-RETV*zqs)
159               zcor = 1./(1.-RETV*zqs)            zqs = zqs*zcor
160               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  
161    
162            ! calculer la fraction nuageuse (processus humide):            ! calculer la fraction nuageuse (processus humide):
163    
# Line 199  contains Line 176  contains
176                    - zdphi/RCPD/(1.+RVTMP2*zq) &                    - zdphi/RCPD/(1.+RVTMP2*zq) &
177                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &                    *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
178                    )*(1.+RETV*q(i, k-1))                    )*(1.+RETV*q(i, k-1))
179               zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))               ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
180               zri(i) = zri(i) &               ri(i) = ri(i) &
181                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &                    + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
182                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
183                    /(zdu2*0.5*(ztvd+ztvu))                    /(zdu2*0.5*(ztvd+ztvu))
184            ELSE            ELSE
185               ! calcul de Ridchardson compatible LMD5               ! calcul de Ridchardson compatible LMD5
186               zri(i) =(RCPD*(t(i, k)-t(i, k-1)) &               ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
187                    -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) &
188                    *(pplay(i, k)-pplay(i, k-1)) &                    *(pplay(i, k)-pplay(i, k-1)) &
189                    )*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)))
190               zri(i) = zri(i) + &               ri(i) = ri(i) + &
191                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &                    zmgeom(i)*zmgeom(i)*gamt(k)/RG &
192                    *(paprs(i, k)/101325.0)**RKAPPA &                    *(paprs(i, k)/101325.0)**RKAPPA &
193                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))                    /(zdu2*0.5*(t(i, k-1)+t(i, k)))
# Line 218  contains Line 195  contains
195    
196            ! finalement, les coefficients d'echange sont obtenus:            ! finalement, les coefficients d'echange sont obtenus:
197    
198            zcdn = SQRT(zdu2) / zmgeom(i) * RG            cdn = SQRT(zdu2) / zmgeom(i) * RG
199    
200            IF (opt_ec) THEN            IF (opt_ec) THEN
201               z2geomf = zgeop(i, k-1)+zgeop(i, k)               z2geomf = zgeop(i, k-1)+zgeop(i, k)
202               zalm2 = (0.5*ckap/RG*z2geomf &               alm2 = (0.5*ckap/RG*z2geomf &
203                    /(1.+0.5*ckap/rg/clam*z2geomf))**2                    /(1.+0.5*ckap/rg/clam*z2geomf))**2
204               zalh2 = (0.5*ckap/rg*z2geomf &               zalh2 = (0.5*ckap/rg*z2geomf &
205                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2                    /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
206               IF (zri(i).LT.0.0) THEN               IF (ri(i) < 0.) THEN
207                  ! situation instable                  ! situation instable
208                  zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &                  scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
209                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
210                  zscf = SQRT(-zri(i)*zscf)                  scf = SQRT(-ri(i)*scf)
211                  zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
212                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
213                  pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
214                  pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
215               ELSE               ELSE
216                  ! situation stable                  ! situation stable
217                  zscf = SQRT(1.+cd*zri(i))                  scf = SQRT(1.+cd*ri(i))
218                  pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
219                  pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
220               ENDIF               ENDIF
221            ELSE            ELSE
222               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)) &
223                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2                    /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
224               pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))               coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
225               pcfm(i, k)= zl2(i)* pcfm(i, k)               coefm(i, k)= l2(i) * coefm(i, k)
226               pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different               coefh(i, k) = coefm(i, k) / prandtl ! h et m different
227            ENDIF            ENDIF
228         ENDDO loop_horiz         ENDDO loop_horiz
229      ENDDO loop_vertical      ENDDO loop_vertical
230    
231      ! Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
232        forall (i = 1: knon)
233      DO i = 1, knon         coefh(i, itop(i) + 1:) = 0.
234         IF (itop(i)+1 .LE. klev) THEN         coefm(i, itop(i) + 1:) = 0.
235            DO k = itop(i)+1, klev      END forall
              pcfh(i, k) = 0.  
              pcfm(i, k) = 0.  
           ENDDO  
        ENDIF  
     ENDDO  
236    
237    END SUBROUTINE coefkz    END SUBROUTINE coefkz
238    
Legend:
Removed from v.57  
changed lines
  Added in v.229
  ViewVC Help
Powered by ViewVC 1.1.21