/[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.f revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/Sources/phylmd/coefkz.f revision 250 by guez, Fri Jan 5 18:18:53 2018 UTC
# Line 1  Line 1 
1        SUBROUTINE coefkz(nsrf, knon, paprs, pplay,  module coefkz_m
2  cIM 261103  
3       .                  ksta, ksta_ter,    IMPLICIT none
4  cIM 261103  
5       .                  ts, rugos,  contains
6       .                  u,v,t,q,  
7       .                  qsurf,    SUBROUTINE coefkz(nsrf, paprs, pplay, ts, u, v, t, q, zgeop, coefm, coefh)
8       .                  pcfm, pcfh)  
9        use dimens_m      ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
10        use indicesol      ! Date: September 22nd, 1993
11        use dimphy  
12        use iniprint      ! Objet : calculer les coefficients d'échange turbulent dans
13        use SUPHEC_M      ! l'atmosphère.
14        use yoethf_m  
15        use fcttre      USE clesphys, ONLY: ksta, ksta_ter
16        use conf_phys_m      USE conf_phys_m, ONLY: iflag_pbl
17        IMPLICIT none      USE dimphy, ONLY: klev
18  c======================================================================      USE fcttre, ONLY: foede, foeew
19  c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922      USE indicesol, ONLY: is_oce
20  c           (une version strictement identique a l'ancien modele)      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
21  c Objet: calculer le coefficient du frottement du sol (Cdrag) et les      USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
22  c        coefficients d'echange turbulent dans l'atmosphere.  
23  c Arguments:      integer, intent(in):: nsrf ! indicateur de la nature du sol
24  c nsrf-----input-I- indicateur de la nature du sol  
25  c knon-----input-I- nombre de points a traiter      REAL, intent(in):: paprs(:, :) ! (knon, klev + 1)
26  c paprs----input-R- pression a chaque intercouche (en Pa)      ! pression a chaque intercouche (en Pa)
27  c pplay----input-R- pression au milieu de chaque couche (en Pa)  
28  c ts-------input-R- temperature du sol (en Kelvin)      real, intent(in):: pplay(:, :) ! (knon, klev)
29  c rugos----input-R- longeur de rugosite (en m)      ! pression au milieu de chaque couche (en Pa)
30  c u--------input-R- vitesse u  
31  c v--------input-R- vitesse v      REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
32  c t--------input-R- temperature (K)      REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev) wind
33  c q--------input-R- vapeur d'eau (kg/kg)      REAL, intent(in):: t(:, :) ! (knon, klev) temperature (K)
34  c      real, intent(in):: q(:, :) ! (knon, klev) vapeur d'eau (kg/kg)
35  c itop-----output-I- numero de couche du sommet de la couche limite      REAL, intent(in):: zgeop(:, :) ! (knon, klev)
36  c pcfm-----output-R- coefficients a calculer (vitesse)      REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse
37  c pcfh-----output-R- coefficients a calculer (chaleur et humidite)  
38  c======================================================================      real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
39  c      ! coefficient, chaleur et humidité
40  c Arguments:  
41  c      ! Local:
42        INTEGER knon, nsrf  
43        REAL ts(klon)      INTEGER knon ! nombre de points a traiter
44        REAL paprs(klon,klev+1), pplay(klon,klev)  
45        REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)      INTEGER itop(size(ts)) ! (knon)
46        REAL rugos(klon)      ! numero de couche du sommet de la couche limite
47  c  
48        REAL pcfm(klon,klev), pcfh(klon,klev)      ! Quelques constantes et options:
49        INTEGER itop(klon)  
50  c      REAL, PARAMETER:: cepdu2 =0.1**2
51  c Quelques constantes et options:      REAL, PARAMETER:: CKAP = 0.4
52  c      REAL, PARAMETER:: cb = 5.
53        REAL cepdu2, ckap, cb, cc, cd, clam      REAL, PARAMETER:: cc = 5.
54        PARAMETER (cepdu2 =(0.1)**2)      REAL, PARAMETER:: cd = 5.
55        PARAMETER (CKAP=0.4)      REAL, PARAMETER:: clam = 160.
56        PARAMETER (cb=5.0)      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
57        PARAMETER (cc=5.0)  
58        PARAMETER (cd=5.0)      LOGICAL, PARAMETER:: richum = .TRUE.
59        PARAMETER (clam=160.0)      ! utilise le nombre de Richardson humide
60        REAL ratqs ! largeur de distribution de vapeur d'eau  
61        PARAMETER (ratqs=0.05)      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
62        LOGICAL richum ! utilise le nombre de Richardson humide      REAL, PARAMETER:: prandtl = 0.4
63        PARAMETER (richum=.TRUE.)  
64        REAL ric ! nombre de Richardson critique      REAL kstable ! diffusion minimale (situation stable)
65        PARAMETER(ric=0.4)      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
66        REAL prandtl      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
67        PARAMETER (prandtl=0.4)  
68        REAL kstable ! diffusion minimale (situation stable)      LOGICAL, PARAMETER:: tvirtu = .TRUE.
69        ! GKtest      ! calculer Ri d'une maniere plus performante
70        ! PARAMETER (kstable=1.0e-10)  
71        REAL ksta, ksta_ter      LOGICAL, PARAMETER:: opt_ec = .FALSE.
72  cIM: 261103     REAL kstable_ter, kstable_sinon      ! formule du Centre Europeen dans l'atmosphere
73  cIM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)  
74  cIM: 261103     PARAMETER (kstable_ter = 1.0e-8)      INTEGER i, k
75  cIM: 261103   PARAMETER (kstable_ter = 1.0e-10)      REAL zmgeom(size(ts))
76  cIM: 261103   PARAMETER (kstable_sinon = 1.0e-10)      REAL ri(size(ts))
77        ! fin GKtest      REAL l2(size(ts))
78        REAL mixlen ! constante controlant longueur de melange      REAL zdphi, zdu2, ztvd, ztvu, cdn
79        PARAMETER (mixlen=35.0)      REAL scf
80        INTEGER isommet ! le sommet de la couche limite      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
81        PARAMETER (isommet=klev)      logical zdelta
82        LOGICAL tvirtu ! calculer Ri d'une maniere plus performante      REAL z2geomf, zalh2, alm2, zscfh, scfm
83        PARAMETER (tvirtu=.TRUE.)      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
84        LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere  
85        PARAMETER (opt_ec=.FALSE.)      !--------------------------------------------------------------------
86    
87  c      knon = size(ts)
88  c Variables locales:  
89  c      ! Prescrire la valeur de contre-gradient
90        INTEGER i, k, kk !IM 120704      if (iflag_pbl.eq.1) then
91        REAL zgeop(klon,klev)         DO k = 3, klev
92        REAL zmgeom(klon)            gamt(k) = - 1E-3
93        REAL zri(klon)         ENDDO
94        REAL zl2(klon)         gamt(2) = - 2.5E-3
95        else
96        REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)         DO k = 2, klev
97        REAL pcfm1(klon), pcfh1(klon)            gamt(k) = 0.0
98  c         ENDDO
99        REAL zdphi, zdu2, ztvd, ztvu, zcdn      ENDIF
100        REAL zscf  
101        REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs      IF (nsrf .NE. is_oce ) THEN
102        REAL z2geomf, zalh2, zalm2, zscfh, zscfm         kstable = ksta_ter
103        REAL t_coup      ELSE
104        PARAMETER (t_coup=273.15)         kstable = ksta
105  cIM      ENDIF
106        LOGICAL check  
107        PARAMETER (check=.false.)      ! Calculer les coefficients turbulents dans l'atmosphere
108  c  
109  c contre-gradient pour la chaleur sensible: Kelvin/metre      itop = isommet
110        REAL gamt(2:klev)  
111        real qsurf(klon)      loop_vertical: DO k = 2, isommet
112  c         loop_horiz: DO i = 1, knon
113        LOGICAL appel1er            zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
114        SAVE appel1er                 + (v(i, k) - v(i, k - 1))**2)
115  c            zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
116  c Fonctions thermodynamiques et fonctions d'instabilite            zdphi = zmgeom(i) / 2.0
117        REAL fsta, fins, x            zt = (t(i, k) + t(i, k - 1)) * 0.5
118        LOGICAL zxli ! utiliser un jeu de fonctions simples            zq = (q(i, k) + q(i, k - 1)) * 0.5
119        PARAMETER (zxli=.FALSE.)  
120  c            ! calculer Qs et dQs/dT:
121        fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))  
122        fins(x) = SQRT(1.0-18.0*x)            zdelta = RTT >=zt
123  c            zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
124        DATA appel1er /.TRUE./                 / (1. + RVTMP2 * zq)
125  c            zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
126        IF (appel1er) THEN            zqs = MIN(0.5, zqs)
127          if (prt_level > 9) THEN            zcor = 1./(1. - RETV * zqs)
128            print *,'coefkz, opt_ec:', opt_ec            zqs = zqs * zcor
129            print *,'coefkz, richum:', richum            zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
130            IF (richum) print *,'coefkz, ratqs:', ratqs  
131            print *,'coefkz, isommet:', isommet            ! calculer la fraction nuageuse (processus humide):
132            print *,'coefkz, tvirtu:', tvirtu  
133            appel1er = .FALSE.            zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
134          endif            zfr = MAX(0.0, MIN(1.0, zfr))
135        ENDIF            IF (.NOT.richum) zfr = 0.0
136  c  
137  c Initialiser les sorties            !  calculer le nombre de Richardson:
138  c  
139        DO k = 1, klev            IF (tvirtu) THEN
140        DO i = 1, knon               ztvd = (t(i, k) &
141           pcfm(i,k) = 0.0                    + zdphi/RCPD/(1. + RVTMP2 * zq) &
142           pcfh(i,k) = 0.0                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs) ) &
143        ENDDO                    ) * (1. + RETV * q(i, k))
144        ENDDO               ztvu = (t(i, k - 1) &
145        DO i = 1, knon                    - zdphi/RCPD/(1. + RVTMP2 * zq) &
146           itop(i) = 0                    * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs) ) &
147        ENDDO                    ) * (1. + RETV * q(i, k - 1))
148                 ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
149  c               ri(i) = ri(i) &
150  c Prescrire la valeur de contre-gradient                    + zmgeom(i) * zmgeom(i)/RG * gamt(k) &
151  c                    * (paprs(i, k)/101325.0)**RKAPPA &
152        if (iflag_pbl.eq.1) then                    /(zdu2 * 0.5 * (ztvd + ztvu))
153           DO k = 3, klev            ELSE
154              gamt(k) = -1.0E-03               ! calcul de Ridchardson compatible LMD5
155           ENDDO               ri(i) = (RCPD * (t(i, k) - t(i, k - 1)) &
156           gamt(2) = -2.5E-03                    - RD * 0.5 * (t(i, k) + t(i, k - 1))/paprs(i, k) &
157        else                    * (pplay(i, k) - pplay(i, k - 1)) &
158           DO k = 2, klev                    ) * zmgeom(i)/(zdu2 * 0.5 * RCPD * (t(i, k - 1) + t(i, k)))
159              gamt(k) = 0.0               ri(i) = ri(i) + &
160           ENDDO                    zmgeom(i) * zmgeom(i) * gamt(k)/RG &
161        ENDIF                    * (paprs(i, k)/101325.0)**RKAPPA &
162  cIM cf JLD/ GKtest                    /(zdu2 * 0.5 * (t(i, k - 1) + t(i, k)))
163        IF ( nsrf .NE. is_oce ) THEN            ENDIF
164  cIM 261103     kstable = kstable_ter  
165          kstable = ksta_ter            ! finalement, les coefficients d'echange sont obtenus:
166        ELSE  
167  cIM 261103     kstable = kstable_sinon            cdn = SQRT(zdu2) / zmgeom(i) * RG
168          kstable = ksta  
       ENDIF  
 cIM cf JLD/ GKtest fin  
 c  
 c Calculer les geopotentiels de chaque couche  
 c  
       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  
 c  
 c Calculer le frottement au sol (Cdrag)  
 c  
       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  
 c  
       CALL clcdrag(klon, knon, nsrf, zxli,  
      $             u1, v1, t1, q1, z1,  
      $             ts, qsurf, rugos,  
      $             pcfm1, pcfh1)  
 cIM  $             ts, qsurf, rugos,  
 C  
       DO i = 1, knon  
        pcfm(i,1)=pcfm1(i)  
        pcfh(i,1)=pcfh1(i)  
       ENDDO  
 c  
 c Calculer les coefficients turbulents dans l'atmosphere  
 c  
       DO i = 1, knon  
          itop(i) = isommet  
       ENDDO  
   
   
       DO k = 2, isommet  
       DO i = 1, knon  
             zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2  
      .                     +(v(i,k)-v(i,k-1))**2)  
             zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)  
             zdphi =zmgeom(i) / 2.0  
             zt = (t(i,k)+t(i,k-1)) * 0.5  
             zq = (q(i,k)+q(i,k-1)) * 0.5  
 c  
 c           calculer Qs et dQs/dT:  
 c  
             IF (thermcep) THEN  
               zdelta = MAX(0.,SIGN(1.,RTT-zt))  
               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  
      .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta  
               zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)  
               zqs = MIN(0.5,zqs)  
               zcor = 1./(1.-RETV*zqs)  
               zqs = 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  
 c  
 c           calculer la fraction nuageuse (processus humide):  
 c  
             zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)  
             zfr = MAX(0.0,MIN(1.0,zfr))  
             IF (.NOT.richum) zfr = 0.0  
 c  
 c           calculer le nombre de Richardson:  
 c  
             IF (tvirtu) THEN  
             ztvd =( t(i,k)  
      .             + zdphi/RCPD/(1.+RVTMP2*zq)  
      .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )  
      .            )*(1.+RETV*q(i,k))  
             ztvu =( t(i,k-1)  
      .             - zdphi/RCPD/(1.+RVTMP2*zq)  
      .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )  
      .            )*(1.+RETV*q(i,k-1))  
             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))  
             zri(i) = zri(i)  
      .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)  
      .               *(paprs(i,k)/101325.0)**RKAPPA  
      .               /(zdu2*0.5*(ztvd+ztvu))  
 c  
             ELSE ! calcul de Ridchardson compatible LMD5  
 c  
             zri(i) =(RCPD*(t(i,k)-t(i,k-1))  
      .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)  
      .               *(pplay(i,k)-pplay(i,k-1))  
      .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))  
             zri(i) = zri(i) +  
      .             zmgeom(i)*zmgeom(i)*gamt(k)/RG  
 cSB     .             /(paprs(i,k)/101325.0)**RKAPPA  
      .             *(paprs(i,k)/101325.0)**RKAPPA  
      .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))  
             ENDIF  
 c  
 c           finalement, les coefficients d'echange sont obtenus:  
 c  
             zcdn=SQRT(zdu2) / zmgeom(i) * RG  
 c  
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  ! situation instable               IF (ri(i) < 0.) THEN
176                 zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3                  ! situation instable
177       .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)                  scf = ((zgeop(i, k)/zgeop(i, k - 1))**(1./3.) - 1.)**3 &
178                 zscf = SQRT(-zri(i)*zscf)                       / (zmgeom(i)/RG)**3 / (zgeop(i, k - 1)/RG)
179                 zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)                  scf = SQRT(- ri(i) * scf)
180                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)                  scfm = 1.0 / (1.0 + 3.0 * cb * cc * alm2 * scf)
181                 pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)                  zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * scf)
182                 pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)                  coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
183              ELSE ! situation stable                  coefh(i, k) = cdn * zalh2 * (1. - 3.0 * cb * ri(i) * zscfh)
184                 zscf=SQRT(1.+cd*zri(i))               ELSE
185                 pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)                  ! situation stable
186                 pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)                  scf = SQRT(1. + cd * ri(i))
187              ENDIF                  coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
188                    coefh(i, k) = cdn * zalh2/(1. + 3.0 * cb * ri(i) * scf)
189                 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         ENDDO loop_horiz
198        ENDDO      ENDDO loop_vertical
199  c  
200  c Au-dela du sommet, pas de diffusion turbulente:      ! Au-delà du sommet, pas de diffusion turbulente :
201  c      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
205                 pcfh(i,k) = 0.0  
206                 pcfm(i,k) = 0.0    END SUBROUTINE coefkz
207              ENDDO  
208           ENDIF  end module coefkz_m
       ENDDO  
 c  
       RETURN  
       END  

Legend:
Removed from v.38  
changed lines
  Added in v.250

  ViewVC Help
Powered by ViewVC 1.1.21