Changeset 3146 for palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
- Timestamp:
- Jul 18, 2018 10:36:19 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r3130 r3146 21 21 ! Current revisions: 22 22 ! ------------------ 23 ! 23 ! Major bugfix in calculation of Obukhov length 24 24 ! 25 25 ! Former revisions: … … 329 329 surf => surf_def_h(l) 330 330 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 332 337 ENDIF 333 338 ENDDO … … 354 359 ENDIF 355 360 356 ! 357 !-- Call for urban-type vertical surfaces 361 !-- Call for urban-type vertical surfaces 358 362 IF ( surf_usm_v(l)%ns >= 1 ) THEN 359 363 surf => surf_usm_v(l) … … 1184 1188 IF ( ibc_pt_b /= 1 ) THEN 1185 1189 IF ( humidity ) THEN 1186 !$OMP PARALLEL DO PRIVATE( i, j, k,z_mo )1190 !$OMP PARALLEL DO PRIVATE( z_mo ) 1187 1191 DO m = 1, surf%ns 1188 1192 1189 i = surf%i(m)1190 j = surf%j(m)1191 k = surf%k(m)1192 1193 1193 z_mo = surf%z_mo(m) 1194 1194 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 ) 1198 1199 ENDDO 1199 1200 ELSE 1200 !$OMP PARALLEL DO PRIVATE( i, j, k,z_mo )1201 !$OMP PARALLEL DO PRIVATE( z_mo ) 1201 1202 DO m = 1, surf%ns 1202 1203 1203 i = surf%i(m)1204 j = surf%j(m)1205 k = surf%k(m)1206 1207 1204 z_mo = surf%z_mo(m) 1208 1205 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 ) 1212 1210 ENDDO 1213 1211 ENDIF 1214 1212 ELSE 1215 1213 IF ( humidity ) THEN 1216 !$OMP PARALLEL DO PRIVATE( i, j,k, z_mo )1214 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1217 1215 DO m = 1, surf%ns 1218 1216 1219 i = surf%i(m)1220 j = surf%j(m)1221 1217 k = surf%k(m) 1222 1218 … … 1224 1220 1225 1221 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) ) * & 1228 1224 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 & 1230 1226 + 1.0E-20_wp ) 1231 1227 ENDDO 1232 1228 ELSE 1233 !$OMP PARALLEL DO PRIVATE( i, j,k, z_mo )1229 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1234 1230 DO m = 1, surf%ns 1235 1231 1236 i = surf%i(m)1237 j = surf%j(m)1238 1232 k = surf%k(m) 1239 1233 … … 1242 1236 surf%rib(m) = - g * z_mo * surf%shf(m) * & 1243 1237 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 & 1245 1239 + 1.0E-20_wp ) 1246 1240 ENDDO … … 1418 1412 1419 1413 IF ( .NOT. humidity ) THEN 1420 !$OMP PARALLEL DO PRIVATE( i, j, k,z_mo )1414 !$OMP PARALLEL DO PRIVATE( z_mo ) 1421 1415 DO m = 1, surf%ns 1422 1416 1423 i = surf%i(m)1424 j = surf%j(m)1425 k = surf%k(m)1426 1427 1417 z_mo = surf%z_mo(m) 1428 1418 1429 surf%ol(m) = ( pt(k,j,i) * surf%us(m)**2 ) /&1419 surf%ol(m) = ( surf%pt1(m) * surf%us(m)**2 ) / & 1430 1420 ( kappa * g * & 1431 1421 surf%ts(m) + 1E-30_wp ) … … 1453 1443 1454 1444 1455 surf%ol(m) = ( vpt(k,j,i) * surf%us(m)**2 ) /&1445 surf%ol(m) = ( surf%vpt1(m) * surf%us(m)**2 ) / & 1456 1446 ( kappa * g * ( surf%ts(m) + & 1457 1447 0.61_wp * surf%pt1(m) * surf%us(m) & … … 1472 1462 ELSE 1473 1463 1474 !$OMP PARALLEL DO PRIVATE( i, j, k,z_mo )1464 !$OMP PARALLEL DO PRIVATE( z_mo ) 1475 1465 DO m = 1, surf%ns 1476 1466 1477 i = surf%i(m)1478 j = surf%j(m)1479 k = surf%k(m)1480 1481 1467 z_mo = surf%z_mo(m) 1482 1468 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) * & 1486 1472 surf%ts(m) ) + 1E-30_wp ) 1487 1473 … … 1560 1546 1561 1547 ! 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 1563 1550 !-- ( only for upward-facing surfs ) 1564 1551 SUBROUTINE calc_pt_q … … 1590 1577 ENDIF 1591 1578 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 1592 1583 ENDDO 1593 1584 … … 1596 1587 1597 1588 ! 1598 !-- Calculate potential temperature a nd specific humidity at first grid level1589 !-- Calculate potential temperature at first grid level 1599 1590 !-- ( only for upward-facing surfs ) 1600 1591 SUBROUTINE calc_pt_surface … … 1602 1593 IMPLICIT NONE 1603 1594 1604 INTEGER(iwp) :: k off !< 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) 1605 1596 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1606 1597 1607 k off = surf%koff1598 k_off = surf%koff 1608 1599 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1609 1600 DO m = 1, surf%ns … … 1613 1604 k = surf%k(m) 1614 1605 1615 surf%pt_surface(m) = pt(k+k off,j,i)1606 surf%pt_surface(m) = pt(k+k_off,j,i) 1616 1607 1617 1608 ENDDO … … 1619 1610 END SUBROUTINE calc_pt_surface 1620 1611 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 1621 1636 ! 1622 1637 !-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
Note: See TracChangeset
for help on using the changeset viewer.