/[lmdze]/trunk/phylmd/Interface_surf/coefkz.f

Diff of /trunk/phylmd/Interface_surf/coefkz.f

Parent Directory | Revision Log | View Patch Patch

-trunk/phylmd/coefkz.f
revision 103 by guez,
Fri Aug 29 13:00:05 2014 UTC
+trunk/phylmd/Interface_surf/coefkz.f
revision 288 by guez,
Tue Jul 24 16:27:12 2018 UTC
 Line 4 
 module coefkz_m
  contains
-   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, coefm, coefh)
      ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
-     ! 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.
+     ! Objet : calculer les coefficients d'échange turbulent dans
+     ! l'atmosphère.
+     USE clesphys, ONLY: ksta, ksta_ter
+     USE conf_phys_m, ONLY: iflag_pbl
+     USE dimphy, ONLY: klev
+     USE fcttre, ONLY: foede, foeew
      USE indicesol, ONLY: is_oce
-     USE dimphy, ONLY: klev, klon
-     USE conf_gcm_m, 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
-     use clcdrag_m, only: clcdrag
-     ! Arguments:
      integer, intent(in):: nsrf ! indicateur de la nature du sol
-     INTEGER, intent(in):: knon ! nombre de points a traiter
-     REAL, intent(in):: paprs(klon, klev+1)
+     REAL, intent(in):: paprs(:, :) ! (knon, klev + 1)
      ! pression a chaque intercouche (en Pa)
-     real, intent(in):: pplay(klon, klev)
+     real, intent(in):: pplay(:, :) ! (knon, klev)
      ! pression au milieu de chaque couche (en Pa)
-     REAL, intent(in):: ksta, ksta_ter
+     REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
-     REAL, intent(in):: ts(klon) ! temperature du sol (en Kelvin)
+     REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev) wind
-     REAL, intent(in):: rugos(klon) ! longeur de rugosite (en m)
+     REAL, intent(in):: t(:, :) ! (knon, klev) temperature (K)
-     REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind
+     real, intent(in):: q(:, :) ! (knon, klev) vapeur d'eau (kg/kg)
-     REAL, intent(in):: t(klon, klev) ! temperature (K)
+     REAL, intent(in):: zgeop(:, :) ! (knon, klev)
-     real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
+     REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev) coefficient, vitesse
-     real, intent(in):: qsurf(klon)
-     REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
-     real, intent(out):: coefh(:, :) ! (knon, klev)
+     real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
      ! coefficient, chaleur et humidité
      ! Local:
-     INTEGER itop(knon) ! numero de couche du sommet de la couche limite
+     INTEGER knon ! nombre de points a traiter
+     INTEGER itop(size(ts)) ! (knon)
+     ! numero de couche du sommet de la couche limite
      ! Quelques constantes et options:
      REAL, PARAMETER:: cepdu2 =0.1**2
-     REAL, PARAMETER:: CKAP = 0.4
-     REAL, PARAMETER:: cb = 5.
-     REAL, PARAMETER:: cc = 5.
-     REAL, PARAMETER:: cd = 5.
-     REAL, PARAMETER:: clam = 160.
      REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
-     LOGICAL, PARAMETER:: richum = .TRUE.
-     ! utilise le nombre de Richardson humide
      REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
      REAL, PARAMETER:: prandtl = 0.4
      REAL kstable ! diffusion minimale (situation stable)
      REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
      INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
+     INTEGER i, k
-     LOGICAL, PARAMETER:: tvirtu = .TRUE.
+     REAL zmgeom(size(ts))
-     ! calculer Ri d'une maniere plus performante
+     REAL ri(size(ts))
+     REAL l2(size(ts))
-     LOGICAL, PARAMETER:: opt_ec = .FALSE.
-     ! formule du Centre Europeen dans l'atmosphere
-     INTEGER i, k, kk
-     REAL zgeop(klon, klev)
-     REAL zmgeom(klon)
-     REAL ri(klon)
-     REAL l2(klon)
-     REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
      REAL zdphi, zdu2, ztvd, ztvu, cdn
-     REAL scf
      REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
      logical zdelta
-     REAL z2geomf, zalh2, alm2, zscfh, scfm
-     REAL, PARAMETER:: t_coup = 273.15
      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
      !--------------------------------------------------------------------
+     knon = size(ts)
      ! Prescrire la valeur de contre-gradient
-     if (iflag_pbl.eq.1) then
+     if (iflag_pbl == 1) then
         DO k = 3, klev
-           gamt(k) = -1.0E-03
+           gamt(k) = - 1E-3
         ENDDO
-        gamt(2) = -2.5E-03
+        gamt(2) = - 2.5E-3
      else
         DO k = 2, klev
            gamt(k) = 0.0
         ENDDO
      ENDIF
-     IF ( nsrf .NE. is_oce ) THEN
+     kstable = merge(ksta, ksta_ter, nsrf == is_oce)
-        kstable = ksta_ter
-     ELSE
-        kstable = ksta
-     ENDIF
-     ! 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, coefm(:, 1), coefh(:, 1))
      ! Calculer les coefficients turbulents dans l'atmosphere
      itop = isommet
-     loop_vertical: DO k = 2, isommet
+     DO k = 2, isommet
-        loop_horiz: DO i = 1, knon
+        DO i = 1, knon
-           zdu2 = MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &
+           zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
-                +(v(i, k)-v(i, k-1))**2)
+                + (v(i, k) - v(i, k - 1))**2)
-           zmgeom(i) = zgeop(i, k)-zgeop(i, k-1)
+           zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
-           zdphi =zmgeom(i) / 2.0
+           zdphi = zmgeom(i) / 2.0
-           zt = (t(i, k)+t(i, k-1)) * 0.5
+           zt = (t(i, k) + t(i, k - 1)) * 0.5
-           zq = (q(i, k)+q(i, k-1)) * 0.5
+           zq = (q(i, k) + q(i, k - 1)) * 0.5
            ! calculer Qs et dQs/dT:
+           zdelta = RTT >=zt
-           IF (thermcep) THEN
+           zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
-              zdelta = RTT >=zt
+                / (1. + RVTMP2 * zq)
-              zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
+           zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
-                   / (1. + RVTMP2*zq)
+           zqs = MIN(0.5, zqs)
-              zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
+           zcor = 1./(1. - RETV * zqs)
-              zqs = MIN(0.5, zqs)
+           zqs = zqs * zcor
-              zcor = 1./(1.-RETV*zqs)
+           zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
-              zqs = zqs*zcor
-              zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
-           ELSE
-              IF (zt  <  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
            ! calculer la fraction nuageuse (processus humide):
+           zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
-           zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
            zfr = MAX(0.0, MIN(1.0, zfr))
-           IF (.NOT.richum) zfr = 0.0
-           !  calculer le nombre de Richardson:
+           ! calculer le nombre de Richardson:
+           ztvd = (t(i, k) &
-           IF (tvirtu) THEN
+                + zdphi/RCPD/(1. + RVTMP2 * zq) &
-              ztvd =( t(i, k) &
+                * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
-                   + zdphi/RCPD/(1.+RVTMP2*zq) &
+                ) * (1. + RETV * q(i, k))
-                   *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
+           ztvu = (t(i, k - 1) &
-                   )*(1.+RETV*q(i, k))
+                - zdphi/RCPD/(1. + RVTMP2 * zq) &
-              ztvu =( t(i, k-1) &
+                * ((1. - zfr) + zfr * (1. + RLVTT * zqs/RD/zt)/(1. + zdqs)) &
-                   - zdphi/RCPD/(1.+RVTMP2*zq) &
+                ) * (1. + RETV * q(i, k - 1))
-                   *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
+           ri(i) = zmgeom(i) * (ztvd - ztvu)/(zdu2 * 0.5 * (ztvd + ztvu))
-                   )*(1.+RETV*q(i, k-1))
+           ri(i) = ri(i) &
-              ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
+                + zmgeom(i) * zmgeom(i)/RG * gamt(k) &
-              ri(i) = ri(i) &
+                * (paprs(i, k)/101325.0)**RKAPPA &
-                   + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
+                /(zdu2 * 0.5 * (ztvd + ztvu))
-                   *(paprs(i, k)/101325.0)**RKAPPA &
-                   /(zdu2*0.5*(ztvd+ztvu))
-           ELSE
-              ! calcul de Ridchardson compatible LMD5
-              ri(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)))
-              ri(i) = ri(i) + &
-                   zmgeom(i)*zmgeom(i)*gamt(k)/RG &
-                   *(paprs(i, k)/101325.0)**RKAPPA &
-                   /(zdu2*0.5*(t(i, k-1)+t(i, k)))
-           ENDIF
            ! finalement, les coefficients d'echange sont obtenus:
            cdn = SQRT(zdu2) / zmgeom(i) * RG
-           IF (opt_ec) THEN
+           l2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, itop(i) + 1)) &
-              z2geomf = zgeop(i, k-1)+zgeop(i, k)
+                /(paprs(i, 2) - paprs(i, itop(i) + 1))))**2
-              alm2 = (0.5*ckap/RG*z2geomf &
+           coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
-                   /(1.+0.5*ckap/rg/clam*z2geomf))**2
+           coefm(i, k) = l2(i) * coefm(i, k)
-              zalh2 = (0.5*ckap/rg*z2geomf &
+           coefh(i, k) = coefm(i, k) / prandtl ! h et m different
-                   /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
+        ENDDO
-              IF (ri(i) < 0.) THEN
+     ENDDO
-                 ! situation instable
-                 scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
-                      / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
-                 scf = SQRT(-ri(i)*scf)
-                 scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
-                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
-                 coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
-                 coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
-              ELSE
-                 ! situation stable
-                 scf = SQRT(1.+cd*ri(i))
-                 coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
-                 coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
-              ENDIF
-           ELSE
-              l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
-                   /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
-              coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
-              coefm(i, k)= l2(i) * coefm(i, k)
-              coefh(i, k) = coefm(i, k) / prandtl ! h et m different
-           ENDIF
-        ENDDO loop_horiz
-     ENDDO loop_vertical
      ! Au-delà du sommet, pas de diffusion turbulente :
      forall (i = 1: knon)

 Legend:



Removed from v.103
 


changed lines


 
Added in v.288
 Legend:



Removed from v.103
 


changed lines


 
Added in v.288
-Removed from v.103
+Added in v.288

	ViewVC Help
Powered by ViewVC 1.1.21