Ignore:
Timestamp:
Jul 18, 2018 10:36:19 PM (6 years ago)
Author:
maronga
Message:

major bugfix in calculation of Obukhov length and subsequent handling of surface information

File:
1 edited

Legend:

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

    r3130 r3146  
    2121! Current revisions:
    2222! ------------------
    23 !
     23! Major bugfix in calculation of Obukhov length
    2424!
    2525! Former revisions:
     
    329329             surf => surf_def_h(l)
    330330             CALL calc_pt_q
    331              IF ( .NOT. neutral )  CALL calc_pt_surface
     331             IF ( .NOT. neutral )  THEN
     332                CALL calc_pt_surface
     333                IF ( humidity )  THEN
     334                   CALL calc_vpt_surface
     335                ENDIF
     336             ENDIF
    332337          ENDIF
    333338       ENDDO
     
    354359          ENDIF
    355360
    356 !
    357 !--    Call for urban-type vertical surfaces
     361!--       Call for urban-type vertical surfaces
    358362          IF ( surf_usm_v(l)%ns >= 1  )  THEN
    359363             surf => surf_usm_v(l)
     
    11841188          IF ( ibc_pt_b /= 1 )  THEN
    11851189             IF ( humidity )  THEN
    1186                 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1190                !$OMP PARALLEL DO PRIVATE( z_mo )
    11871191                DO  m = 1, surf%ns
    11881192
    1189                    i   = surf%i(m)           
    1190                    j   = surf%j(m)
    1191                    k   = surf%k(m)
    1192 
    11931193                   z_mo = surf%z_mo(m)
    11941194
    1195                    surf%rib(m) = g * z_mo *                                    &
    1196                                          ( vpt(k,j,i) - vpt(k-1,j,i) ) /       &
    1197                       ( surf%uvw_abs(m)**2 * vpt(k,j,i) + 1.0E-20_wp )
     1195                   surf%rib(m) = g * z_mo                                      &
     1196                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
     1197                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
     1198                                 + 1.0E-20_wp )
    11981199                ENDDO
    11991200             ELSE
    1200                 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1201                !$OMP PARALLEL DO PRIVATE( z_mo )
    12011202                DO  m = 1, surf%ns
    12021203
    1203                    i   = surf%i(m)           
    1204                    j   = surf%j(m)
    1205                    k   = surf%k(m)
    1206 
    12071204                   z_mo = surf%z_mo(m)
    12081205
    1209                    surf%rib(m) = g * z_mo *                                    &
    1210                                          ( pt(k,j,i) - pt(k-1,j,i)   ) /       &
    1211                       ( surf%uvw_abs(m)**2 * pt(k,j,i)  + 1.0E-20_wp )
     1206                   surf%rib(m) = g * z_mo                                      &
     1207                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
     1208                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
     1209                                 + 1.0E-20_wp )
    12121210                ENDDO
    12131211             ENDIF
    12141212          ELSE
    12151213             IF ( humidity )  THEN
    1216                 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1214                !$OMP PARALLEL DO PRIVATE( k, z_mo )
    12171215                DO  m = 1, surf%ns
    12181216
    1219                    i   = surf%i(m)           
    1220                    j   = surf%j(m)
    12211217                   k   = surf%k(m)
    12221218
     
    12241220
    12251221                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
    1226                            * q(k,j,i) ) * surf%shf(m) + 0.61_wp                &
    1227                            * pt(k,j,i) * surf%qsws(m) ) *                      &
     1222                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
     1223                           * surf%pt1(m) * surf%qsws(m) ) *                    &
    12281224                             drho_air_zw(k-1)                /                 &
    1229                          ( surf%uvw_abs(m)**3 * vpt(k,j,i) * kappa**2          &
     1225                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
    12301226                           + 1.0E-20_wp )
    12311227                ENDDO
    12321228             ELSE
    1233                 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1229                !$OMP PARALLEL DO PRIVATE( k, z_mo )
    12341230                DO  m = 1, surf%ns
    12351231
    1236                    i   = surf%i(m)           
    1237                    j   = surf%j(m)
    12381232                   k   = surf%k(m)
    12391233
     
    12421236                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
    12431237                                        drho_air_zw(k-1)            /          &
    1244                         ( surf%uvw_abs(m)**3 * pt(k,j,i) * kappa**2            &
     1238                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
    12451239                           + 1.0E-20_wp )
    12461240                ENDDO
     
    14181412
    14191413          IF ( .NOT. humidity )  THEN
    1420              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1414             !$OMP PARALLEL DO PRIVATE( z_mo )
    14211415             DO  m = 1, surf%ns
    14221416
    1423                 i   = surf%i(m)           
    1424                 j   = surf%j(m)
    1425                 k   = surf%k(m)
    1426 
    14271417                z_mo = surf%z_mo(m)
    14281418
    1429                 surf%ol(m) =  ( pt(k,j,i) *  surf%us(m)**2 ) /                 &
     1419                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
    14301420                                       ( kappa * g *                           &
    14311421                                         surf%ts(m) + 1E-30_wp )
     
    14531443
    14541444
    1455                 surf%ol(m) =  ( vpt(k,j,i) * surf%us(m)**2 ) /                 &
     1445                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
    14561446                    ( kappa * g * ( surf%ts(m) +                               &
    14571447                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
     
    14721462          ELSE
    14731463
    1474              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1464             !$OMP PARALLEL DO PRIVATE( z_mo )
    14751465             DO  m = 1, surf%ns
    14761466
    1477                 i   = surf%i(m)           
    1478                 j   = surf%j(m)
    1479                 k   = surf%k(m)
    1480 
    14811467                z_mo = surf%z_mo(m)
    14821468
    1483                 surf%ol(m) =  ( vpt(k,j,i) *  surf%us(m)**2 ) /                &
    1484                     ( kappa * g * ( surf%ts(m) + 0.61_wp * pt(k,j,i) *         &
    1485                                     surf%qs(m) + 0.61_wp * q(k,j,i)  *         &
     1469                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
     1470                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
     1471                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
    14861472                                    surf%ts(m) ) + 1E-30_wp )
    14871473
     
    15601546
    15611547!
    1562 !-- Calculate potential temperature and specific humidity at first grid level
     1548!-- Calculate potential temperature, specific humidity, and virtual potential
     1549!-- temperature at first grid level
    15631550!-- ( only for upward-facing surfs )
    15641551    SUBROUTINE calc_pt_q
     
    15901577          ENDIF
    15911578
     1579          IF ( humidity )  THEN
     1580             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
     1581          ENDIF
     1582         
    15921583       ENDDO
    15931584
     
    15961587
    15971588!
    1598 !-- Calculate potential temperature and specific humidity at first grid level
     1589!-- Calculate potential temperature at first grid level
    15991590!-- ( only for upward-facing surfs )
    16001591    SUBROUTINE calc_pt_surface
     
    16021593       IMPLICIT NONE
    16031594
    1604        INTEGER(iwp) ::  koff    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
     1595       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
    16051596       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
    16061597       
    1607        koff = surf%koff
     1598       k_off = surf%koff
    16081599       !$OMP PARALLEL DO PRIVATE( i, j, k )
    16091600       DO  m = 1, surf%ns
     
    16131604          k   = surf%k(m)
    16141605
    1615           surf%pt_surface(m) = pt(k+koff,j,i)
     1606          surf%pt_surface(m) = pt(k+k_off,j,i)
    16161607
    16171608       ENDDO
     
    16191610    END SUBROUTINE calc_pt_surface
    16201611
     1612!
     1613!-- Calculate virtual potential temperature at first grid level
     1614!-- ( only for upward-facing surfs )
     1615    SUBROUTINE calc_vpt_surface
     1616
     1617       IMPLICIT NONE
     1618
     1619       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
     1620       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
     1621       
     1622       k_off = surf%koff
     1623       !$OMP PARALLEL DO PRIVATE( i, j, k )
     1624       DO  m = 1, surf%ns
     1625
     1626          i   = surf%i(m)           
     1627          j   = surf%j(m)
     1628          k   = surf%k(m)
     1629
     1630          surf%vpt_surface(m) = vpt(k+k_off,j,i)
     1631
     1632       ENDDO
     1633
     1634    END SUBROUTINE calc_vpt_surface
     1635   
    16211636!
    16221637!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
Note: See TracChangeset for help on using the changeset viewer.