/[lmdze]/trunk/phylmd/pbl_surface.f

Diff of /trunk/phylmd/pbl_surface.f

Parent Directory | Revision Log | View Patch Patch

-revision 226 by guez,
Mon Oct 16 13:04:05 2017 UTC
+revision 249 by guez,
Fri Jan  5 17:15:05 2018 UTC
 Line 8 
 contains
         cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &
         qsurf, evap, falbe, fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, &
         agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, &
-        flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, t2m, q2m, &
+        flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, coefh, t2m, q2m, &
         u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, therm, trmb1, &
         trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
 Line 21 
 contains
      ! ne tient pas compte de la diff\'erentiation des sous-fractions
      ! de sol.
+     use clcdrag_m, only: clcdrag
      use clqh_m, only: clqh
      use clvent_m, only: clvent
      use coefkz_m, only: coefkz
      use coefkzmin_m, only: coefkzmin
-     USE conf_gcm_m, ONLY: prt_level, lmt_pas
+     use coefkz2_m, only: coefkz2
+     USE conf_gcm_m, ONLY: lmt_pas
      USE conf_phys_m, ONLY: iflag_pbl
      USE dimphy, ONLY: klev, klon, zmasq
      USE dimsoil, ONLY: nsoilmx
-Line 36 
 contains
+Line 38 
 contains
      USE suphec_m, ONLY: rd, rg, rkappa
      use time_phylmdz, only: itap
      use ustarhb_m, only: ustarhb
-     use vdif_kcay_m, only: vdif_kcay
      use yamada4_m, only: yamada4
      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
-Line 96 
 contains
+Line 97 
 contains
      ! flux de vapeur d'eau (kg / m2 / s) à la surface
      REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
-     ! tension du vent à la surface, en Pa
+     ! tension du vent (flux turbulent de vent) à la surface, en Pa
      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
      real q2(klon, klev + 1, nbsrf)
-Line 106 
 contains
+Line 107 
 contains
      ! dflux_q derive du flux latent
      ! IM "slab" ocean
-     REAL, intent(out):: ycoefh(klon, klev)
+     REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
-     ! "ycoefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
+     ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
      ! ce champ sur les quatre sous-surfaces du mod\`ele.
      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
-Line 150 
 contains
+Line 151 
 contains
      real y_run_off_lic_0(klon)
      REAL rugmer(klon)
      REAL ytsoil(klon, nsoilmx)
-     REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
+     REAL yts(klon), ypct(klon), yz0_new(klon)
+     real yrugos(klon) ! longeur de rugosite (en m)
      REAL yalb(klon)
-     REAL u1lay(klon), v1lay(klon) ! vent dans la premi\`ere couche, pour
-                               ! une sous-surface donnée
      REAL snow(klon), yqsurf(klon), yagesno(klon)
      real yqsol(klon) ! column-density of water in soil, in kg m-2
      REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
-Line 168 
 contains
+Line 166 
 contains
      REAL y_flux_t(klon), y_flux_q(klon)
      REAL y_flux_u(klon), y_flux_v(klon)
      REAL y_dflux_t(klon), y_dflux_q(klon)
-     REAL coefh(klon, klev), coefm(klon, klev)
+     REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
+     real ycdragh(klon), ycdragm(klon)
      REAL yu(klon, klev), yv(klon, klev)
      REAL yt(klon, klev), yq(klon, klev)
      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
+     REAL ycoefm0(klon, 2:klev), ycoefh0(klon, 2:klev)
-     REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
+     REAL yzlay(klon, klev), zlev(klon, klev + 1), yteta(klon, klev)
-     REAL yzlay(klon, klev), yzlev(klon, klev + 1), yteta(klon, klev)
-     REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)
-     REAL ykmq(klon, klev + 1)
      REAL yq2(klon, klev + 1)
-     REAL q2diag(klon, klev + 1)
      REAL delp(klon, klev)
      INTEGER i, k, nsrf
      INTEGER ni(klon), knon, j
      REAL pctsrf_pot(klon, nbsrf)
      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
      ! apparitions ou disparitions de la glace de mer
-     REAL yt2m(klon), yq2m(klon), yu10m(klon)
+     REAL yt2m(klon), yq2m(klon), wind10m(klon)
-     REAL yustar(klon)
+     REAL ustar(klon)
      REAL yt10m(klon), yq10m(klon)
      REAL ypblh(klon)
-Line 204 
 contains
+Line 196 
 contains
      REAL ytrmb1(klon)
      REAL ytrmb2(klon)
      REAL ytrmb3(klon)
-     REAL uzon(klon), vmer(klon)
+     REAL u1(klon), v1(klon)
      REAL tair1(klon), qair1(klon), tairsol(klon)
      REAL psfce(klon), patm(klon)
      REAL qairsol(klon), zgeo1(klon)
      REAL rugo1(klon)
+     REAL zgeop(klon, klev)
-     ! utiliser un jeu de fonctions simples
-     LOGICAL zxli
-     PARAMETER (zxli=.FALSE.)
      !------------------------------------------------------------
-Line 256 
 contains
+Line 245 
 contains
      d_q = 0.
      d_u = 0.
      d_v = 0.
-     ycoefh = 0.
+     coefh = 0.
      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
-Line 300 
 contains
+Line 289 
 contains
               yagesno(j) = agesno(i, nsrf)
               yrugos(j) = frugs(i, nsrf)
               yrugoro(j) = rugoro(i)
-              u1lay(j) = u(i, 1)
-              v1lay(j) = v(i, 1)
               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
               ypaprs(j, klev + 1) = paprs(i, klev + 1)
               y_run_off_lic_0(j) = run_off_lic_0(i)
-Line 325 
 contains
+Line 312 
 contains
               END DO
            END DO
-           ! calculer Cdrag et les coefficients d'echange
+           ! Calculer les géopotentiels de chaque couche:
-           CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
-                yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &
+           zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
-                coefh(:knon, :))
+                + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
+           DO k = 2, klev
+              zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
+                   * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
+                   * (ypplay(:knon, k - 1) - ypplay(:knon, k))
+           ENDDO
+           CALL clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
+                yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
+                yrugos(:knon), ycdragm(:knon), ycdragh(:knon))
            IF (iflag_pbl == 1) THEN
-              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
+              ycdragm(:knon) = max(ycdragm(:knon), 0.)
-              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
+              ycdragh(:knon) = max(ycdragh(:knon), 0.)
-              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
+           end IF
+           ! on met un seuil pour ycdragm et ycdragh
+           IF (nsrf == is_oce) THEN
+              ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
+              ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
            END IF
-           ! on met un seuil pour coefm et coefh
+           CALL coefkz(nsrf, ypaprs(:knon, :), ypplay(:knon, :), ksta, &
-           IF (nsrf == is_oce) THEN
+                ksta_ter, yts(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
-              coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
+                yq(:knon, :), zgeop(:knon, :), ycoefm(:knon, :), &
-              coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
+                ycoefh(:knon, :))
+           IF (iflag_pbl == 1) THEN
+              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, :), &
+                   ycoefh0(:knon, :))
+              ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
+              ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
            END IF
            IF (ok_kzmin) THEN
               ! Calcul d'une diffusion minimale pour les conditions tres stables
               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
-                   coefm(:knon, 1), ycoefm0, ycoefh0)
+                   ycdragm(:knon), ycoefh0(:knon, :))
-              coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
+              ycoefm0(:knon, :) = ycoefh0(:knon, :)
-              coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
+              ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
+              ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
            END IF
-           IF (iflag_pbl >= 3) THEN
+           IF (iflag_pbl >= 6) THEN
               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
               ! Fr\'ed\'eric Hourdin
               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
                    + ypplay(:knon, 1))) &
                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
               DO k = 2, klev
-                 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
+                 yzlay(:knon, k) = yzlay(:knon, k-1) &
                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
                       / ypaprs(1:knon, k) &
                       * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
               END DO
               DO k = 1, klev
                  yteta(1:knon, k) = yt(1:knon, k) * (ypaprs(1:knon, 1) &
                       / ypplay(1:knon, k))**rkappa * (1. + 0.61 * yq(1:knon, k))
               END DO
-              yzlev(1:knon, 1) = 0.
-              yzlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
+              zlev(:knon, 1) = 0.
+              zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
                    - yzlay(:knon, klev - 1)
               DO k = 2, klev
-                 yzlev(1:knon, k) = 0.5 * (yzlay(1:knon, k) + yzlay(1:knon, k-1))
+                 zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(:knon, k-1))
               END DO
               DO k = 1, klev + 1
                  DO j = 1, knon
                     i = ni(j)
-Line 378 
 contains
+Line 393 
 contains
                  END DO
               END DO
-              CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
+              ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))
-              IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
+              CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
+                   yu(:knon, :), yv(:knon, :), yteta(:knon, :), yq2(:knon, :), &
-              ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
+                   ycoefm(:knon, :), ycoefh(:knon, :), ustar(:knon))
-              IF (iflag_pbl >= 11) THEN
-                 CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
-                      yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
-                      iflag_pbl)
-              ELSE
-                 CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
-                      coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
-              END IF
-              coefm(:knon, 2:) = ykmm(:knon, 2:klev)
-              coefh(:knon, 2:) = ykmn(:knon, 2:klev)
            END IF
-           ! calculer la diffusion des vitesses "u" et "v"
+           CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
-           CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &
+                ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
-                coefm(:knon, :), yt, yu, ypaprs, ypplay, ydelp, y_d_u, &
+                ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
                 y_flux_u(:knon))
-           CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &
+           CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
-                coefm(:knon, :), yt, yv, ypaprs, ypplay, ydelp, y_d_v, &
+                ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
+                ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
                 y_flux_v(:knon))
            ! calculer la diffusion de "q" et de "h"
            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
                 ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
-                u1lay(:knon), v1lay(:knon), coefh(:knon, :), yt, yq, &
+                yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
-                yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), yalb(:knon), &
+                yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
-                snow(:knon), yqsurf, yrain_f, ysnow_f, yfluxlat(:knon), &
+                yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
-                pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &
+                yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
-                yz0_new, y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &
+                y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
-                y_dflux_q(:knon), y_fqcalving, y_ffonte, y_run_off_lic_0)
+                y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
+                y_run_off_lic_0)
            ! calculer la longueur de rugosite sur ocean
            yrugm = 0.
            IF (nsrf == is_oce) THEN
               DO j = 1, knon
-                 yrugm(j) = 0.018 * coefm(j, 1) * (u1lay(j)**2 + v1lay(j)**2) &
+                 yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
                       / rg + 0.11 * 14E-6 &
-                      / sqrt(coefm(j, 1) * (u1lay(j)**2 + v1lay(j)**2))
+                      / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
                  yrugm(j) = max(1.5E-05, yrugm(j))
               END DO
            END IF
-Line 432 
 contains
+Line 437 
 contains
            DO k = 1, klev
               DO j = 1, knon
                  i = ni(j)
-                 coefh(j, k) = coefh(j, k) * ypct(j)
-                 coefm(j, k) = coefm(j, k) * ypct(j)
                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
-Line 467 
 contains
+Line 470 
 contains
               agesno(i, nsrf) = yagesno(j)
               fqcalving(i, nsrf) = y_fqcalving(j)
               ffonte(i, nsrf) = y_ffonte(j)
-              cdragh(i) = cdragh(i) + coefh(j, 1)
+              cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
-              cdragm(i) = cdragm(i) + coefm(j, 1)
+              cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
               dflux_t(i) = dflux_t(i) + y_dflux_t(j)
               dflux_q(i) = dflux_q(i) + y_dflux_q(j)
            END DO
-Line 491 
 contains
+Line 494 
 contains
                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)
                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
-                 ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
               END DO
            END DO
+           forall (k = 2:klev) coefh(ni(:knon), k) &
+                = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
            ! diagnostic t, q a 2m et u, v a 10m
            DO j = 1, knon
               i = ni(j)
-              uzon(j) = yu(j, 1) + y_d_u(j, 1)
+              u1(j) = yu(j, 1) + y_d_u(j, 1)
-              vmer(j) = yv(j, 1) + y_d_v(j, 1)
+              v1(j) = yv(j, 1) + y_d_v(j, 1)
               tair1(j) = yt(j, 1) + y_d_t(j, 1)
               qair1(j) = yq(j, 1) + y_d_q(j, 1)
               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
-Line 516 
 contains
+Line 521 
 contains
               qairsol(j) = yqsurf(j)
            END DO
-           CALL stdlevvar(klon, knon, nsrf, zxli, uzon(:knon), vmer(:knon), &
+           CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
-                tair1, qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, &
+                qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
-                yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+                yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
            DO j = 1, knon
               i = ni(j)
               t2m(i, nsrf) = yt2m(j)
               q2m(i, nsrf) = yq2m(j)
-              u10m_srf(i, nsrf) = (yu10m(j) * uzon(j)) &
+              u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
-                   / sqrt(uzon(j)**2 + vmer(j)**2)
+                   / sqrt(u1(j)**2 + v1(j)**2)
-              v10m_srf(i, nsrf) = (yu10m(j) * vmer(j)) &
+              v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
-                   / sqrt(uzon(j)**2 + vmer(j)**2)
+                   / sqrt(u1(j)**2 + v1(j)**2)
            END DO
-           CALL hbtm(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &
+           CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)

 Legend:



Removed from v.226
 


changed lines


 
Added in v.249
 Legend:



Removed from v.226
 


changed lines


 
Added in v.249
-Removed from v.226
+Added in v.249

	ViewVC Help
Powered by ViewVC 1.1.21