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

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

Parent Directory | Revision Log | View Patch Patch

-trunk/libf/phylmd/coefkz.f90
revision 47 by guez,
Fri Jul  1 15:00:48 2011 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, pcfm, pcfh)
-     ! Authors: F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS)
+     ! 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.
-     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
-     ! Arguments:
+     ! 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 suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
+     USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
      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(inout):: pcfm(klon, klev) ! coefficient, vitesse
+     real, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
-     real, intent(inout):: pcfh(klon, klev) ! coefficient, chaleur et humidité
+     ! coefficient, chaleur et humidité
      ! Local:
-     INTEGER itop(klon) ! 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:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
-     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 controlant longueur de melange
+     REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
-     INTEGER, PARAMETER:: isommet = klev ! le sommet de la couche limite
+     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.
+     REAL zdphi, zdu2, ztvd, ztvu, cdn
-     ! formule du Centre Europeen dans l'atmosphere
+     REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
+     logical zdelta
-     INTEGER i, k, kk
-     REAL zgeop(klon, klev)
-     REAL zmgeom(klon)
-     REAL zri(klon)
-     REAL zl2(klon)
-     REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
-     REAL pcfm1(klon), pcfh1(klon)
-     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
      REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
      !--------------------------------------------------------------------
-     ! 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
      ! 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, pcfm1, pcfh1)
-     DO i = 1, knon
-        pcfm(i, 1) = pcfm1(i)
-        pcfh(i, 1) = pcfh1(i)
-     ENDDO
      ! Calculer les coefficients turbulents dans l'atmosphere
-     DO i = 1, knon
+     itop = isommet
-        itop(i) = isommet
-     ENDDO
-     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 = MAX(0., SIGN(1., RTT-zt))
+                / (1. + RVTMP2 * zq)
-              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)  &
+           zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
-                   + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
+           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 .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
            ! 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) &
-              zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
+                + zmgeom(i) * zmgeom(i)/RG * gamt(k) &
-              zri(i) = zri(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
-              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 &
-                   *(paprs(i, k)/101325.0)**RKAPPA &
-                   /(zdu2*0.5*(t(i, k-1)+t(i, k)))
-           ENDIF
            ! finalement, les coefficients d'echange sont obtenus:
-           zcdn = SQRT(zdu2) / zmgeom(i) * RG
+           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
-              zalm2 = (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 (zri(i).LT.0.0) THEN
-                 ! situation instable
-                 zscf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
-                      / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
-                 zscf = SQRT(-zri(i)*zscf)
-                 zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
-                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
-                 pcfm(i, k) = zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
-                 pcfh(i, k) = zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
-              ELSE
-                 ! situation stable
-                 zscf = SQRT(1.+cd*zri(i))
-                 pcfm(i, k) = zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
-                 pcfh(i, k) = zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
-              ENDIF
-           ELSE
-              zl2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
-                   /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
-              pcfm(i, k) = sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
-              pcfm(i, k)= zl2(i)* pcfm(i, k)
-              pcfh(i, k) = pcfm(i, k) /prandtl ! h et m different
-           ENDIF
-        ENDDO loop_horiz
-     ENDDO loop_vertical
-     ! Au-dela du sommet, pas de diffusion turbulente:
-     DO i = 1, knon
-        IF (itop(i)+1 .LE. klev) THEN
-           DO k = itop(i)+1, klev
-              pcfh(i, k) = 0.
-              pcfm(i, k) = 0.
-           ENDDO
-        ENDIF
      ENDDO
+     ! Au-delà du sommet, pas de diffusion turbulente :
+     forall (i = 1: knon)
+        coefh(i, itop(i) + 1:) = 0.
+        coefm(i, itop(i) + 1:) = 0.
+     END forall
    END SUBROUTINE coefkz
  end module coefkz_m

 Legend:



Removed from v.47
 


changed lines


 
Added in v.288
 Legend:



Removed from v.47
 


changed lines


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

	ViewVC Help
Powered by ViewVC 1.1.21