Changeset 4593 for palm


Ignore:
Timestamp:
Jul 9, 2020 12:48:18 PM (4 years ago)
Author:
suehring
Message:

Surface flux calculation: Pre-calculate ln(z/z0) at each timestep in order to reduce the number of log-calculations; Bugfix - add missing density to fluxes of passive-scalars, chemistry and cloud-phyiscal quantities at upward-facing surfaces; Move if-statement out of inner loop; Remove unnecessary index referencing

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r4562 r4593  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! - Pre-calculate ln(z/z0) at each timestep in order to reduce the number of log-calculations
     28! - Bugfix - add missing density to fluxes of passive-scalars, chemistry and cloud-phyiscal
     29!   quantities at upward-facing surfaces
     30! - Move if-statement out of inner loop
     31! - Remove unnecessary index referencing
     32!
     33! 4562 2020-06-12 08:38:47Z raasch
    2734! File re-formatted to follow the PALM coding standard
    2835!
     
    268275    downward      = .FALSE. !< flag indicating downward-facing surface elements
    269276!
     277!-- First, precalculate ln(z/z0) for all surfaces. This is done each timestep, in order to account
     278!-- for time-dependent roughness or user-modifications.
     279    DO  l = 0, 1
     280       IF ( surf_def_h(l)%ns >= 1 )  THEN
     281          surf => surf_def_h(l)
     282          CALL calc_ln
     283       ENDIF
     284    ENDDO
     285    IF ( surf_lsm_h%ns >= 1 )  THEN
     286       surf => surf_lsm_h
     287       CALL calc_ln
     288    ENDIF
     289    IF ( surf_usm_h%ns >= 1 )  THEN
     290       surf => surf_usm_h
     291       CALL calc_ln
     292    ENDIF
     293
     294    DO  l = 0, 3
     295       IF ( surf_def_v(l)%ns >= 1 )  THEN
     296          surf => surf_def_v(l)
     297          CALL calc_ln
     298       ENDIF
     299       IF ( surf_lsm_v(l)%ns >= 1 )  THEN
     300          surf => surf_lsm_v(l)
     301          CALL calc_ln
     302       ENDIF
     303       IF ( surf_usm_v(l)%ns >= 1 )  THEN
     304          surf => surf_usm_v(l)
     305          CALL calc_ln
     306       ENDIF
     307    ENDDO
     308!
    270309!-- Derive potential temperature and specific humidity at first grid level from the fields pt and q
    271 !
    272310!-- First call for horizontal default-type surfaces (l=0 - upward facing, l=1 - downward facing)
    273311    DO  l = 0, 1
     
    314352
    315353!
    316 !-- First, calculate the new Obukhov length, then new friction velocity, followed by the new scaling
     354!-- First, calculate the new Obukhov length from precalculated values of log(z/z0) and wind speeds.
     355!-- As a second step, then calculate new friction velocity, followed by the new scaling
    317356!-- parameters (th*, q*, etc.), and the new surface fluxes, if required. Note, each routine is called
    318357!-- for different surface types. First call for default-type horizontal surfaces, for natural- and
    319358!-- urban-type surfaces. Note, here only upward-facing horizontal surfaces are treated.
    320 
    321 !
    322 !-- Default-type upward-facing horizontal surfaces
     359!-- Note, calculation of log(z/z0) is redone each timestep, in order to account for time-dependent
     360!-- values.
     361!-- Start with default-type upward-facing horizontal surfaces
    323362    IF ( surf_def_h(0)%ns >= 1 )  THEN
    324363       surf => surf_def_h(0)
     
    645684! Description:
    646685! ------------
     686!> Compute ln(z/z0).
     687!--------------------------------------------------------------------------------------------------!
     688 SUBROUTINE calc_ln
     689
     690    INTEGER(iwp) ::  m     !< running index surface elements
     691
     692!
     693!-- Note, ln(z/z0h) and ln(z/z0q) is also calculated even if neural simulations are applied.
     694!-- This is because the scalar coefficients are also used for other scalars such as passive scalars,
     695!-- chemistry and aerosols.
     696    !$OMP PARALLEL  DO PRIVATE( z_mo )
     697    !$ACC PARALLEL LOOP PRIVATE(z_mo) &
     698    !$ACC PRESENT(surf)
     699    DO  m = 1, surf%ns
     700       z_mo = surf%z_mo(m)
     701       surf%ln_z_z0(m)  = LOG( z_mo / surf%z0(m) )
     702       surf%ln_z_z0h(m) = LOG( z_mo / surf%z0h(m) )
     703       surf%ln_z_z0q(m) = LOG( z_mo / surf%z0q(m) )
     704    ENDDO
     705
     706 END SUBROUTINE calc_ln
     707 
     708!--------------------------------------------------------------------------------------------------!
     709! Description:
     710! ------------
    647711!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
    648712!> surface elements. This is required by all methods.
     
    663727    ibit = MERGE( 1, 0, .NOT. downward )
    664728
    665     !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
    666     !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
    667     !$ACC PRESENT(surf, u, v)
    668     DO  m = 1, surf%ns
    669        i   = surf%i(m)
    670        j   = surf%j(m)
    671        k   = surf%k(m)
    672 
    673 !
    674 !--    Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This
    675 !--    will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is
    676 !--    used to set w_lfc to zero under stable conditions.
    677        IF ( use_free_convection_scaling )  THEN
     729    IF ( use_free_convection_scaling )  THEN
     730       !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
     731       !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
     732       !$ACC PRESENT(surf, u, v)
     733       DO  m = 1, surf%ns
     734          i   = surf%i(m)
     735          j   = surf%j(m)
     736          k   = surf%k(m)
     737
     738!
     739!--       Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This
     740!--       will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is
     741!--       used to set w_lfc to zero under stable conditions.
    678742          w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
    679743          w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
    680        ELSE
    681           w_lfc = 0.0_wp
    682        ENDIF
    683 
    684 !
    685 !--    Compute the absolute value of the horizontal velocity. (relative to the surface in case the
    686 !--    lower surface is the ocean). Please note, in new surface modelling concept the index values
    687 !--    changed, i.e. the reference grid point is not the surface-grid point itself but the first
    688 !--    grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
    689 !--    relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
    690 !--    to calculate the absolute velocity. However, this does not apply for downward-facing walls.
    691 !--    To mask this, use ibit, which checks for upward/downward-facing surfaces.
    692        surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) )  &
    693                                             * ibit ) )**2                                          &
    694                                + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
    695                                                                         ) * ibit )                 &
     744!
     745!--       Compute the absolute value of the horizontal velocity. (relative to the surface in case the
     746!--       lower surface is the ocean). Please note, in new surface modelling concept the index values
     747!--       changed, i.e. the reference grid point is not the surface-grid point itself but the first
     748!--       grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
     749!--       relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
     750!--       to calculate the absolute velocity. However, this does not apply for downward-facing walls.
     751!--       To mask this, use ibit, which checks for upward/downward-facing surfaces.
     752          surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1)  &
     753                                                                       ) * ibit )                   &
     754                                  )**2                                                              &
     755                                + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
     756                                                                       ) * ibit )                   &
    696757                                  )**2  + w_lfc**2 )
    697     ENDDO
     758       ENDDO
     759    ELSE
     760       !$OMP PARALLEL DO PRIVATE(i, j, k)
     761       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
     762       !$ACC PRESENT(surf, u, v)
     763       DO  m = 1, surf%ns
     764          i   = surf%i(m)
     765          j   = surf%j(m)
     766          k   = surf%k(m)
     767!
     768!--       Compute the absolute value of the horizontal velocity. (relative to the surface in case the
     769!--       lower surface is the ocean). Please note, in new surface modelling concept the index values
     770!--       changed, i.e. the reference grid point is not the surface-grid point itself but the first
     771!--       grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
     772!--       relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
     773!--       to calculate the absolute velocity. However, this does not apply for downward-facing walls.
     774!--       To mask this, use ibit, which checks for upward/downward-facing surfaces.
     775          surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1)  &
     776                                                                       ) * ibit )                   &
     777                                  )**2                                                              &
     778                                + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i)  &
     779                                                                       ) * ibit )                   &
     780                                  )**2 )
     781       ENDDO
     782    ENDIF
    698783
    699784 END SUBROUTINE calc_uvw_abs
     
    9801065!
    9811066!--             Calculate f = Ri - [...]/[...]^2 = 0
    982                 f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) )                    &
     1067                f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_z0h(m)                             &
    9831068                                                      - psi_h( z_mo / ol_m )                       &
    9841069                                                      + psi_h( surf%z0h(m) / ol_m ) ) /            &
    985                     ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m )                              &
    986                       + psi_m( surf%z0(m) /  ol_m ) )**2
     1070                    ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m )                                       &
     1071                                      + psi_m( surf%z0(m) /  ol_m ) )**2
    9871072
    9881073!
    9891074!--             Calculate df/dL
    990                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) )                         &
     1075                f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_z0h(m)                                  &
    9911076                                                 - psi_h( z_mo / ol_u )                            &
    9921077                                                 + psi_h( surf%z0h(m) / ol_u ) ) /                 &
    993                            ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_u )                       &
    994                              + psi_m( surf%z0(m) / ol_u ) )**2                                     &
    995                            + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / ol_l )  &
    996                                                  + psi_h( surf%z0h(m) / ol_l ) ) /                 &
    997                            ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_l )                       &
    998                              + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )
     1078                           ( surf%ln_z_z0(m) - psi_m( z_mo / ol_u )                                &
     1079                                             + psi_m( surf%z0(m) / ol_u ) )**2                     &
     1080                           + ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l )           &
     1081                                                                  + psi_h( surf%z0h(m) / ol_l ) ) /&
     1082                           ( surf%ln_z_z0(m) - psi_m( z_mo / ol_l )                                &
     1083                                             + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )
    9991084             ELSE
    10001085!
    10011086!--             Calculate f = Ri - 1 /[...]^3 = 0
    10021087                f = surf%rib(m) - ( z_mo / ol_m ) /                                                &
    1003                     ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) &
    1004                     )**3
     1088                    ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) )**3
    10051089
    10061090!
    10071091!--             Calculate df/dL
    1008                 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )                          &
     1092                f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_z0(m)                                   &
    10091093                                                 - psi_m( z_mo / ol_u )                            &
    10101094                                                 + psi_m( surf%z0(m) / ol_u )                      &
    10111095                                               )**3                                                &
    1012                            + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )                          &
     1096                           + ( z_mo / ol_l ) / ( surf%ln_z_z0(m)                                   &
    10131097                                                 - psi_m( z_mo / ol_l )                            &
    10141098                                                 + psi_m( surf%z0(m) / ol_l )                      &
     
    10951179!
    10961180!--             Calculate f = Ri - [...]/[...]^2 = 0
    1097                 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )      &
     1181                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( surf%ln_z_z0(m)                       &
    10981182                                                           - psi_h( z_mo_vec(m) / ol_m )           &
    10991183                                                           + psi_h( surf%z0h(m) / ol_m )           &
    11001184                                                           ) /                                     &
    1101                                                            ( LOG( z_mo_vec(m) / surf%z0(m) )       &
     1185                                                           ( surf%ln_z_z0(m)                       &
    11021186                                                          - psi_m( z_mo_vec(m) / ol_m )            &
    11031187                                                          + psi_m( surf%z0(m) /  ol_m )            &
     
    11061190!
    11071191!--             Calculate df/dL
    1108                 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )           &
     1192                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( surf%ln_z_z0h(m)                           &
    11091193                                                        - psi_h( z_mo_vec(m) / ol_u )              &
    11101194                                                        + psi_h( surf%z0h(m) / ol_u )              &
    11111195                                                      ) /                                          &
    1112                                                       ( LOG( z_mo_vec(m) / surf%z0(m) )            &
     1196                                                      ( surf%ln_z_z0(m)                            &
    11131197                                                        - psi_m( z_mo_vec(m) / ol_u )              &
    11141198                                                        + psi_m( surf%z0(m) / ol_u )               &
    11151199                                                      )**2                                         &
    1116                            + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )           &
     1200                           + ( z_mo_vec(m) / ol_l ) * ( surf%ln_z_z0h(m)                           &
    11171201                                                        - psi_h( z_mo_vec(m) / ol_l )              &
    11181202                                                        + psi_h( surf%z0h(m) / ol_l )              &
    11191203                                                      ) /                                          &
    1120                                                       ( LOG( z_mo_vec(m) / surf%z0(m) )            &
     1204                                                      ( surf%ln_z_z0(m)                            &
    11211205                                                        - psi_m( z_mo_vec(m) / ol_l )              &
    11221206                                                        + psi_m( surf%z0(m) / ol_l )               &
     
    11261210!
    11271211!--             Calculate f = Ri - 1 /[...]^3 = 0
    1128                 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) )       &
     1212                f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( surf%ln_z_z0(m)                       &
    11291213                                                             - psi_m( z_mo_vec(m) / ol_m )         &
    11301214                                                             + psi_m( surf%z0(m) / ol_m )          &
     
    11331217!
    11341218!--             Calculate df/dL
    1135                 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) )            &
     1219                f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( surf%ln_z_z0(m)                            &
    11361220                                                        - psi_m( z_mo_vec(m) / ol_u )              &
    11371221                                                        + psi_m( surf%z0(m) / ol_u )               &
    11381222                                                      )**3                                         &
    1139                            + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) )            &
     1223                           + ( z_mo_vec(m) / ol_l ) / ( surf%ln_z_z0(m)                            &
    11401224                                                        - psi_m( z_mo_vec(m) / ol_l )              &
    11411225                                                        + psi_m( surf%z0(m) / ol_l )               &
     
    12101294!
    12111295!--          Compute u* at the scalars' grid points
    1212              surf%us(m) = kappa * surf%uvw_abs(m) / ( LOG( z_mo / surf%z0(m) )                     &
     1296             surf%us(m) = kappa * surf%uvw_abs(m) / ( surf%ln_z_z0(m)                              &
    12131297                                                      - psi_m( z_mo / surf%ol(m) )                 &
    12141298                                                      + psi_m( surf%z0(m) / surf%ol(m) ) )
     
    12171301!--    Compute u* at downward-facing surfaces. This case, do not consider any stability.
    12181302       ELSE
    1219           !$OMP PARALLEL  DO PRIVATE( z_mo )
    1220           !$ACC PARALLEL LOOP PRIVATE(z_mo) &
     1303          !$OMP PARALLEL  DO
     1304          !$ACC PARALLEL LOOP &
    12211305          !$ACC PRESENT(surf)
    12221306          DO  m = 1, surf%ns
    1223              z_mo = surf%z_mo(m)
    12241307!
    12251308!--          Compute u* at the scalars' grid points
    1226              surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
     1309             surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
    12271310          ENDDO
    12281311       ENDIF
     
    12311314!-- No stability is considered in this case.
    12321315    ELSE
    1233        !$OMP PARALLEL DO PRIVATE( z_mo )
    1234        !$ACC PARALLEL LOOP PRIVATE(z_mo) &
     1316       !$OMP PARALLEL DO
     1317       !$ACC PARALLEL LOOP &
    12351318       !$ACC PRESENT(surf)
    12361319       DO  m = 1, surf%ns
    1237           z_mo = surf%z_mo(m)
    1238           surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
     1320          surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
    12391321       ENDDO
    12401322    ENDIF
     
    14111493          z_mo = surf%z_mo(m)
    14121494          surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )                                &
    1413                        / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )                  &
    1414                                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1495                       / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) )                           &
     1496                                            + psi_h( surf%z0h(m) / surf%ol(m) ) )
    14151497       ENDDO
    14161498
     
    14821564                z_mo = surf%z_mo(m)
    14831565                surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) )                           &
    1484                              / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )            &
    1485                                                            + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1566                             / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                     &
     1567                                                  + psi_h( surf%z0q(m) / surf%ol(m) ) )
    14861568             ENDDO
    14871569          ELSE
     
    14931575                z_mo = surf%z_mo(m)
    14941576                surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )                                     &
    1495                              / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )            &
    1496                                                            + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1577                             / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                     &
     1578                                                  + psi_h( surf%z0q(m) / surf%ol(m) ) )
    14971579             ENDDO
    14981580          ENDIF
     
    15361618
    15371619             surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) )                                        &
    1538                           / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )               &
    1539                                                         + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1620                          / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) )                        &
     1621                                               + psi_h( surf%z0h(m) / surf%ol(m) ) )
    15401622          ENDDO
    15411623       ENDIF
     
    15941676
    15951677          surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                                        &
    1596                         / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
    1597                                                       + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1678                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
     1679                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
    15981680
    15991681          surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                                        &
    1600                         / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
    1601                                                       + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1682                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
     1683                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
    16021684       ENDDO
    16031685
     
    16161698
    16171699          surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                                        &
    1618                         / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
    1619                                                       + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1700                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
     1701                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
    16201702
    16211703          surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                                        &
    1622                         / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )                 &
    1623                                                       + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1704                        / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) )                          &
     1705                                             + psi_h( surf%z0q(m) / surf%ol(m) ) )
    16241706       ENDDO
    16251707
     
    16711753
    16721754             surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )                                      &
    1673                             / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )              &
    1674                                                          + psi_m( surf%z0(m) / surf%ol(m) ) )
     1755                            / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) )                       &
     1756                                                + psi_m( surf%z0(m) / surf%ol(m) ) )
    16751757!
    16761758!--          Please note, the computation of usws is not fully accurate. Actually a further
     
    16841766!--    At downward-facing surfaces
    16851767       ELSE
    1686           !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1768          !$OMP PARALLEL DO PRIVATE( i, j, k )
    16871769          DO  m = 1, surf%ns
    16881770             i = surf%i(m)
     
    16901772             k = surf%k(m)
    16911773
    1692              z_mo = surf%z_mo(m)
    1693 
    1694              surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
     1774             surf%usws(m) = kappa * u(k,j,i) / surf%ln_z_z0(m)
    16951775             surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
    16961776          ENDDO
     
    17131793
    17141794             surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )                                      &
    1715                             / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )              &
    1716                                                          + psi_m( surf%z0(m) / surf%ol(m) ) )
     1795                            / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) )                       &
     1796                                                + psi_m( surf%z0(m) / surf%ol(m) ) )
    17171797!
    17181798!--          Please note, the computation of vsws is not fully accurate. Actually a further
     
    17261806!--    Downward-facing surfaces
    17271807       ELSE
    1728           !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1808          !$OMP PARALLEL DO PRIVATE( i, j, k )
    17291809          DO  m = 1, surf%ns
    17301810             i = surf%i(m)
     
    17321812             k = surf%k(m)
    17331813
    1734              z_mo = surf%z_mo(m)
    1735 
    1736              surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
     1814             surf%vsws(m) = kappa * v(k,j,i) / surf%ln_z_z0(m)
    17371815             surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
    17381816          ENDDO
    17391817       ENDIF
    17401818!
    1741 !--    Compute the vertical kinematic heat flux
     1819!--    Compute the vertical kinematic heat flux. Note, only upward-facing surfaces are considered,
     1820!--    at downward-facing surfaces the flux is not parametrized with a scaling parameter.
    17421821       IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point <=  skip_time_do_lsm  &
    17431822             .AND.  simulated_time > 0.0_wp )  .OR.  .NOT.  land_surface )  .AND.                  &
    17441823             .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
    1745           !$OMP PARALLEL DO PRIVATE( i, j, k )
    1746           DO  m = 1, surf%ns
    1747              i = surf%i(m)
    1748              j = surf%j(m)
     1824          !$OMP PARALLEL DO PRIVATE( k )
     1825          DO  m = 1, surf%ns
    17491826             k = surf%k(m)
    17501827             surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
     
    17571834            .OR.  .NOT.  land_surface  )  .AND.  .NOT.  urban_surface  .AND.  .NOT.  downward )    &
    17581835            THEN
    1759           !$OMP PARALLEL DO PRIVATE( i, j, k )
    1760           DO  m = 1, surf%ns
    1761              i    = surf%i(m)
    1762              j    = surf%j(m)
     1836          !$OMP PARALLEL DO PRIVATE( k )
     1837          DO  m = 1, surf%ns
    17631838             k    = surf%k(m)
    17641839             surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
     
    17681843!--    Compute the vertical scalar flux
    17691844       IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.  .NOT.  downward )  THEN
    1770           !$OMP PARALLEL DO PRIVATE( i, j )
    1771           DO  m = 1, surf%ns
    1772              i = surf%i(m)
    1773              j = surf%j(m)
    1774              surf%ssws(m) = -surf%ss(m) * surf%us(m)
     1845          !$OMP PARALLEL DO
     1846          DO  m = 1, surf%ns
     1847             surf%ssws(m) = -surf%ss(m) * surf%us(m) * rho_air_zw(k-1)
    17751848          ENDDO
    17761849       ENDIF
     
    17791852       DO  lsp = 1, nvar
    17801853          IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.  .NOT.  downward )  THEN
    1781              !$OMP PARALLEL DO PRIVATE( i, j )
     1854             !$OMP PARALLEL DO
    17821855             DO  m = 1, surf%ns
    1783                 i = surf%i(m)
    1784                 j = surf%j(m)
    1785                 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
     1856                surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) * rho_air_zw(k-1)
    17861857             ENDDO
    17871858          ENDIF
     
    17911862!--    Compute (turbulent) fluxes of cloud water content and cloud drop conc.
    17921863       IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.  .NOT. downward)  THEN
    1793           !$OMP PARALLEL DO PRIVATE( i, j )
    1794           DO  m = 1, surf%ns
    1795              i = surf%i(m)
    1796              j = surf%j(m)
    1797 
    1798              surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
    1799              surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
     1864          !$OMP PARALLEL DO
     1865          DO  m = 1, surf%ns
     1866             surf%qcsws(m) = -surf%qcs(m) * surf%us(m) * rho_air_zw(k-1)
     1867             surf%ncsws(m) = -surf%ncs(m) * surf%us(m) * rho_air_zw(k-1)
    18001868          ENDDO
    18011869       ENDIF
     
    18031871!--    Compute (turbulent) fluxes of rain water content and rain drop conc.
    18041872       IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.  .NOT. downward)  THEN
    1805           !$OMP PARALLEL DO PRIVATE( i, j )
    1806           DO  m = 1, surf%ns
    1807              i = surf%i(m)
    1808              j = surf%j(m)
    1809 
    1810              surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
    1811              surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
     1873          !$OMP PARALLEL DO
     1874          DO  m = 1, surf%ns
     1875             surf%qrsws(m) = -surf%qrs(m) * surf%us(m) * rho_air_zw(k-1)
     1876             surf%nrsws(m) = -surf%nrs(m) * surf%us(m) * rho_air_zw(k-1)
    18121877          ENDDO
    18131878       ENDIF
     
    18321897!
    18331898!-- Calcuate surface fluxes at vertical surfaces. No stability is considered.
     1899!-- Further, no density needs to be considered here.
    18341900    ELSE
    18351901!
     
    18411907          flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )
    18421908          flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )
    1843           !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
     1909          !$OMP PARALLEL  DO PRIVATE( i, j, k )
    18441910          DO  m = 1, surf%ns
    18451911             i = surf%i(m)
     
    18471913             k = surf%k(m)
    18481914
    1849              z_mo = surf%z_mo(m)
    1850 
    18511915             surf%mom_flux_uv(m) = kappa * ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /            &
    1852                                    LOG( z_mo / surf%z0(m) )
     1916                                   surf%ln_z_z0(m)
    18531917
    18541918            surf%mom_flux_uv(m) = - surf%mom_flux_uv(m) * surf%us(m)
     
    18581922!--    Compute wsus l={0,1} and wsvs l={2,3}
    18591923       IF ( mom_w )  THEN
    1860           !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
     1924          !$OMP PARALLEL  DO PRIVATE( i, j, k )
    18611925          DO  m = 1, surf%ns
    18621926             i = surf%i(m)
     
    18641928             k = surf%k(m)
    18651929
    1866              z_mo = surf%z_mo(m)
    1867 
    1868              surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
     1930             surf%mom_flux_w(m) = kappa * w(k,j,i) / surf%ln_z_z0(m)
    18691931
    18701932             surf%mom_flux_w(m) = - surf%mom_flux_w(m) * surf%us(m)
     
    19071969          ENDIF
    19081970
    1909           !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
     1971          !$OMP PARALLEL DO PRIVATE( i, j, dum )
    19101972          DO  m = 1, surf%ns
    19111973             i = surf%i(m)
    19121974             j = surf%j(m)
    19131975
    1914              z_mo = surf%z_mo(m)
    1915 
    1916              dum = kappa / LOG( z_mo / surf%z0(m) )
     1976             dum = kappa / surf%ln_z_z0(m)
    19171977!
    19181978!--          usvs (l=0,1) and vsus (l=2,3)
     
    19462006    CHARACTER(LEN = *), INTENT(IN) ::  z_char  !< string identifier to identify z level
    19472007
    1948     INTEGER(iwp) ::  i  !< grid index x-dimension
    1949     INTEGER(iwp) ::  j  !< grid index y-dimension
    1950     INTEGER(iwp) ::  k  !< grid index z-dimension
    19512008    INTEGER(iwp) ::  m  !< running index for surface elements
    19522009
    1953 
    19542010    SELECT CASE ( z_char)
    19552011
     
    19572013
    19582014          DO  m = 1, surf%ns
    1959 
    1960              i = surf%i(m)
    1961              j = surf%j(m)
    1962              k = surf%k(m)
    1963 
    19642015             surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa                             &
    19652016                               * ( LOG( 0.1_wp /  surf%z0h(m) ) - psi_h( 0.1_wp / surf%ol(m) )     &
  • palm/trunk/SOURCE/surface_mod.f90

    r4586 r4593  
    2525! -----------------
    2626! $Id$
     27! Add arrays for pre-calculated ln(z/z0)
     28!
     29! 4586 2020-07-01 16:16:43Z gronemeier
    2730! renamed Richardson flux number into gradient Richardson number (1D model)
    2831!
     
    248251
    249252       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css    !< scaling parameter chemical species
     253!
     254!--    Pre-defined arrays for ln(z/z0)
     255       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0  !< ln(z/z0)
     256       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0h !< ln(z/z0h)
     257       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0q !< ln(z/z0q)
    250258!
    251259!--    Define arrays for surface fluxes
     
    12811289    DEALLOCATE ( surfaces%uvw_abs )
    12821290!
     1291!-- Pre-calculated ln(z/z0)
     1292    DEALLOCATE ( surfaces%ln_z_z0  )
     1293    DEALLOCATE ( surfaces%ln_z_z0h )
     1294    DEALLOCATE ( surfaces%ln_z_z0q )
     1295!
    12831296!-- Roughness
    12841297    DEALLOCATE ( surfaces%z0 )
     
    14111424    ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
    14121425!
     1426!-- Precalculated ln(z/z0)
     1427    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns)  )
     1428    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
     1429    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
     1430!
    14131431!-- Roughness
    14141432    ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
     
    15201538    !$ACC DELETE(surfaces%z_mo(1:surfaces%ns)) &
    15211539    !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
     1540    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns)) &
     1541    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns)) &
     1542    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns)) &
    15221543    !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
    15231544    !$ACC COPYOUT(surfaces%us(1:surfaces%ns)) &
     
    15611582    !$ACC COPYIN(surfaces%z_mo(1:surfaces%ns)) &
    15621583    !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
     1584    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns)) &
     1585    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns)) &
     1586    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns)) &
    15631587    !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
    15641588    !$ACC COPYIN(surfaces%us(1:surfaces%ns)) &
     
    18301854    DEALLOCATE ( surfaces%uvw_abs )
    18311855!
     1856!-- Precalculated ln(z/z0)
     1857    DEALLOCATE ( surfaces%ln_z_z0  )
     1858    DEALLOCATE ( surfaces%ln_z_z0h )
     1859    DEALLOCATE ( surfaces%ln_z_z0q )
     1860!
    18321861!-- Roughness
    18331862    DEALLOCATE ( surfaces%z0 )
    18341863    DEALLOCATE ( surfaces%z0h )
    18351864    DEALLOCATE ( surfaces%z0q )
    1836 
    18371865!
    18381866!-- Friction velocity
     
    18791907!-- Scaling parameter (cs*) and surface flux of chemical species
    18801908    IF ( air_chemistry )  THEN
    1881           DEALLOCATE ( surfaces%css )
    1882           DEALLOCATE ( surfaces%cssws )
     1909       DEALLOCATE ( surfaces%css )
     1910       DEALLOCATE ( surfaces%cssws )
    18831911    ENDIF
    18841912!
     
    19581986    ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
    19591987!
     1988!-- Precalculated ln(z/z0)
     1989    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns) )
     1990    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
     1991    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
     1992!
    19601993!-- Roughness
    19611994    ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
     
    20072040!-- Scaling parameter (cs*) and surface flux of chemical species
    20082041    IF ( air_chemistry )  THEN
    2009           ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
    2010           ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
     2042       ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
     2043       ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
    20112044    ENDIF
    20122045!
     
    20622095    !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
    20632096    !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
     2097    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns) ) &
     2098    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns) ) &
     2099    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns) ) &
    20642100    !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
    20652101    !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
     
    20942130    !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
    20952131    !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
     2132    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns) ) &
     2133    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns) ) &
     2134    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns) ) &
    20962135    !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
    20972136    !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
     
    25692608    surf%rib(num_h)     = 0.0_wp
    25702609    surf%uvw_abs(num_h) = 0.0_wp
     2610!
     2611!-- Initialize ln(z/z0)
     2612    surf%ln_z_z0(num_h)  = LOG( surf%z_mo(num_h) / surf%z0(num_h) )
     2613    surf%ln_z_z0h(num_h) = LOG( surf%z_mo(num_h) / surf%z0h(num_h) )
     2614    surf%ln_z_z0q(num_h) = LOG( surf%z_mo(num_h) / surf%z0q(num_h) )
    25712615
    25722616    IF ( .NOT. constant_diffusion )  THEN
     
    28712915
    28722916    surf%us(num_v)  = 0.0_wp
     2917!
     2918!-- Initialize ln(z/z0)
     2919    surf%ln_z_z0(num_v)  = LOG( surf%z_mo(num_v) / surf%z0(num_v)  )
     2920    surf%ln_z_z0h(num_v) = LOG( surf%z_mo(num_v) / surf%z0h(num_v) )
     2921    surf%ln_z_z0q(num_v) = LOG( surf%z_mo(num_v) / surf%z0q(num_v) )
    28732922!
    28742923!-- If required, initialize Obukhov length
     
    51065155    ENDIF
    51075156    IF ( INDEX( restart_string(1:length), '%cssws' ) /= 0 )  THEN
    5108        IF ( ALLOCATED( surf_target%cssws )  .AND.  ALLOCATED( surf_file%cssws   ) )  THEN
     5157       IF ( ALLOCATED( surf_target%cssws )  .AND.  ALLOCATED( surf_file%cssws ) )  THEN
    51095158          DO  lsp = 1, nvar
    51105159             surf_target%cssws(lsp,m_target) = surf_file%cssws(lsp,m_file)
     
    55075556!
    55085557!-- Redistribute surface elements on its respective type.
    5509 
    55105558    DO  l = 0, 2
    55115559       CALL restore_surface_elements( surf_def_h(l), surf_h(l) )
    5512        CALL restore_surface_elements( surf_lsm_h, surf_h(l) )
    5513        CALL restore_surface_elements( surf_usm_h, surf_h(l) )
    55145560    ENDDO
     5561    CALL restore_surface_elements( surf_lsm_h, surf_h(0) )
     5562    CALL restore_surface_elements( surf_usm_h, surf_h(0) )
    55155563
    55165564    DO  l = 0, 3
Note: See TracChangeset for help on using the changeset viewer.