Changeset 3802 for palm/trunk/SOURCE
- Timestamp:
- Mar 17, 2019 1:33:42 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
random_function_mod.f90 (modified) (2 diffs)
-
urban_surface_mod.f90 (modified) (2 diffs)
-
vertical_nesting_mod.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/random_function_mod.f90
r3655 r3802 25 25 ! ----------------- 26 26 ! $Id$ 27 ! type conversion added to avoid compiler warning about constant integer 28 ! division truncation 29 ! 30 ! 3655 2019-01-07 16:51:22Z knoop 27 31 ! Corrected "Former revisions" section 28 32 ! … … 131 135 132 136 PARAMETER ( ia=16807, im=2147483647, am=1.0_wp/im, iq=127773, ir=2836, & 133 ntab=32, ndiv=1+ (im-1)/ntab, eps=1.2e-7_wp, rnmx=1.0_wp-eps )137 ntab=32, ndiv=1+INT(REAL(im-1)/ntab), eps=1.2e-7_wp, rnmx=1.0_wp-eps ) 134 138 135 139 IF ( idum .le. 0 .or. random_iy .eq. 0 ) THEN -
palm/trunk/SOURCE/urban_surface_mod.f90
r3769 r3802 28 28 ! ----------------- 29 29 ! $Id$ 30 ! unused subroutine commented out 31 ! 32 ! 3769 2019-02-28 10:16:49Z moh.hefny 30 33 ! removed unused variables 31 34 ! … … 8629 8632 8630 8633 8631 CONTAINS8632 !------------------------------------------------------------------------------!8633 ! Description:8634 ! ------------8635 !> Calculation of specific humidity of the skin layer (surface). It is assumend8636 !> that the skin is always saturated.8637 !------------------------------------------------------------------------------!8638 SUBROUTINE calc_q_surface_usm8639 8640 IMPLICIT NONE8641 8642 REAL(wp) :: resistance !< aerodynamic and soil resistance term8643 8644 DO m = 1, surf_usm_h%ns8645 8646 i = surf_usm_h%i(m)8647 j = surf_usm_h%j(m)8648 k = surf_usm_h%k(m)8649 8650 ! 8651 ! -- Calculate water vapour pressure at saturation8652 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * &8653 ( t_surf_green_h_p(m) - 273.16_wp ) / &8654 ( t_surf_green_h_p(m) - 35.86_wp ) &8655 )8656 8657 ! 8658 ! -- Calculate specific humidity at saturation8659 q_s = 0.622_wp * e_s / ( surface_pressure - e_s )8660 8661 ! surf_usm_h%r_a_green(m) = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exner(k) ) / &8662 ! ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )8663 ! 8664 ! !-- make sure that the resistance does not drop to zero8665 ! IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp ) surf_usm_h%r_a_green(m) = 1.0E-10_wp8666 8667 resistance = surf_usm_h%r_a_green(m) / ( surf_usm_h%r_a_green(m) + surf_usm_h%r_s(m) + 1E-5_wp )8668 8669 ! 8670 ! -- Calculate specific humidity at surface8671 IF ( bulk_cloud_model ) THEN8672 q(k,j,i) = resistance * q_s + &8673 ( 1.0_wp - resistance ) * &8674 ( q(k,j,i) - ql(k,j,i) )8675 ELSE8676 q(k,j,i) = resistance * q_s + &8677 ( 1.0_wp - resistance ) * &8678 q(k,j,i)8679 ENDIF8680 8681 ! 8682 ! -- Update virtual potential temperature8683 vpt(k,j,i) = pt(k,j,i) * &8684 ( 1.0_wp + 0.61_wp * q(k,j,i) )8685 8686 ENDDO8687 8688 ! 8689 ! -- Now, treat vertical surface elements8690 DO l = 0, 38691 DO m = 1, surf_usm_v(l)%ns8692 ! 8693 ! -- Get indices of respective grid point8694 i = surf_usm_v(l)%i(m)8695 j = surf_usm_v(l)%j(m)8696 k = surf_usm_v(l)%k(m)8697 8698 ! 8699 ! -- Calculate water vapour pressure at saturation8700 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * &8701 ( t_surf_green_v_p(l)%t(m) - 273.16_wp ) / &8702 ( t_surf_green_v_p(l)%t(m) - 35.86_wp ) &8703 )8704 8705 ! 8706 ! -- Calculate specific humidity at saturation8707 q_s = 0.622_wp * e_s / ( surface_pressure -e_s )8708 8709 ! 8710 ! -- Calculate specific humidity at surface8711 IF ( bulk_cloud_model ) THEN8712 q(k,j,i) = ( q(k,j,i) - ql(k,j,i) )8713 ELSE8714 q(k,j,i) = q(k,j,i)8715 ENDIF8716 ! 8717 ! -- Update virtual potential temperature8718 vpt(k,j,i) = pt(k,j,i) * &8719 ( 1.0_wp + 0.61_wp * q(k,j,i) )8720 8721 ENDDO8722 8723 ENDDO8724 8725 END SUBROUTINE calc_q_surface_usm8634 ! CONTAINS 8635 ! !------------------------------------------------------------------------------! 8636 ! ! Description: 8637 ! ! ------------ 8638 ! !> Calculation of specific humidity of the skin layer (surface). It is assumend 8639 ! !> that the skin is always saturated. 8640 ! !------------------------------------------------------------------------------! 8641 ! SUBROUTINE calc_q_surface_usm 8642 ! 8643 ! IMPLICIT NONE 8644 ! 8645 ! REAL(wp) :: resistance !< aerodynamic and soil resistance term 8646 ! 8647 ! DO m = 1, surf_usm_h%ns 8648 ! 8649 ! i = surf_usm_h%i(m) 8650 ! j = surf_usm_h%j(m) 8651 ! k = surf_usm_h%k(m) 8652 ! 8653 !! 8654 !!-- Calculate water vapour pressure at saturation 8655 ! e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * & 8656 ! ( t_surf_green_h_p(m) - 273.16_wp ) / & 8657 ! ( t_surf_green_h_p(m) - 35.86_wp ) & 8658 ! ) 8659 ! 8660 !! 8661 !!-- Calculate specific humidity at saturation 8662 ! q_s = 0.622_wp * e_s / ( surface_pressure - e_s ) 8663 ! 8664 !! surf_usm_h%r_a_green(m) = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exner(k) ) / & 8665 !! ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp ) 8666 !! 8667 !! !-- make sure that the resistance does not drop to zero 8668 !! IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp ) surf_usm_h%r_a_green(m) = 1.0E-10_wp 8669 ! 8670 ! resistance = surf_usm_h%r_a_green(m) / ( surf_usm_h%r_a_green(m) + surf_usm_h%r_s(m) + 1E-5_wp ) 8671 ! 8672 !! 8673 !!-- Calculate specific humidity at surface 8674 ! IF ( bulk_cloud_model ) THEN 8675 ! q(k,j,i) = resistance * q_s + & 8676 ! ( 1.0_wp - resistance ) * & 8677 ! ( q(k,j,i) - ql(k,j,i) ) 8678 ! ELSE 8679 ! q(k,j,i) = resistance * q_s + & 8680 ! ( 1.0_wp - resistance ) * & 8681 ! q(k,j,i) 8682 ! ENDIF 8683 ! 8684 !! 8685 !!-- Update virtual potential temperature 8686 ! vpt(k,j,i) = pt(k,j,i) * & 8687 ! ( 1.0_wp + 0.61_wp * q(k,j,i) ) 8688 ! 8689 ! ENDDO 8690 ! 8691 !! 8692 !!-- Now, treat vertical surface elements 8693 ! DO l = 0, 3 8694 ! DO m = 1, surf_usm_v(l)%ns 8695 !! 8696 !!-- Get indices of respective grid point 8697 ! i = surf_usm_v(l)%i(m) 8698 ! j = surf_usm_v(l)%j(m) 8699 ! k = surf_usm_v(l)%k(m) 8700 ! 8701 !! 8702 !!-- Calculate water vapour pressure at saturation 8703 ! e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * & 8704 ! ( t_surf_green_v_p(l)%t(m) - 273.16_wp ) / & 8705 ! ( t_surf_green_v_p(l)%t(m) - 35.86_wp ) & 8706 ! ) 8707 ! 8708 !! 8709 !!-- Calculate specific humidity at saturation 8710 ! q_s = 0.622_wp * e_s / ( surface_pressure -e_s ) 8711 ! 8712 !! 8713 !!-- Calculate specific humidity at surface 8714 ! IF ( bulk_cloud_model ) THEN 8715 ! q(k,j,i) = ( q(k,j,i) - ql(k,j,i) ) 8716 ! ELSE 8717 ! q(k,j,i) = q(k,j,i) 8718 ! ENDIF 8719 !! 8720 !!-- Update virtual potential temperature 8721 ! vpt(k,j,i) = pt(k,j,i) * & 8722 ! ( 1.0_wp + 0.61_wp * q(k,j,i) ) 8723 ! 8724 ! ENDDO 8725 ! 8726 ! ENDDO 8727 ! 8728 ! END SUBROUTINE calc_q_surface_usm 8726 8729 8727 8730 END SUBROUTINE usm_surface_energy_balance -
palm/trunk/SOURCE/vertical_nesting_mod.f90
r3655 r3802 26 26 ! ----------------- 27 27 ! $Id$ 28 ! unused subroutines commented out 29 ! 30 ! 3655 2019-01-07 16:51:22Z knoop 28 31 ! unused variables removed 29 32 ! … … 1513 1516 1514 1517 1515 SUBROUTINE interpolate_to_fine_flux1516 1517 1518 USE arrays_3d1519 USE control_parameters1520 USE grid_variables1521 USE indices1522 USE pegrid1523 1524 1525 IMPLICIT NONE1526 1527 INTEGER(iwp) :: i1528 INTEGER(iwp) :: j1529 INTEGER(iwp) :: iif1530 INTEGER(iwp) :: jjf1531 INTEGER(iwp) :: bottomx1532 INTEGER(iwp) :: bottomy1533 INTEGER(iwp) :: topx1534 INTEGER(iwp) :: topy1535 REAL(wp) :: eps1536 REAL(wp) :: alpha1537 REAL(wp) :: eminus1538 REAL(wp) :: edot1539 REAL(wp) :: eplus1540 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uswspr1541 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vswspr1542 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tspr1543 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uspr1544 1545 ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )1546 ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )1547 ALLOCATE( tspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )1548 ALLOCATE( uspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )1549 1550 !1551 !-- Initialisation of scalar variables (2D)1552 1553 !1554 !-- Interpolation in x-direction (quadratic, Clark and Farley)1555 1556 DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+11557 DO i = bdims_rem(1,1), bdims_rem(1,2)1558 1559 bottomx = (nxf+1)/(nxc+1) * i1560 topx = (nxf+1)/(nxc+1) * (i+1) - 11561 1562 DO iif = bottomx, topx1563 1564 eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc1565 alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.01566 eminus = eps * ( eps - 1.0 ) / 2.0 + alpha1567 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha1568 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha1569 1570 uswspr(j,iif) = eminus * work2dusws(j,i-1) &1571 + edot * work2dusws(j,i) &1572 + eplus * work2dusws(j,i+1)1573 1574 vswspr(j,iif) = eminus * work2dvsws(j,i-1) &1575 + edot * work2dvsws(j,i) &1576 + eplus * work2dvsws(j,i+1)1577 1578 tspr(j,iif) = eminus * work2dts(j,i-1) &1579 + edot * work2dts(j,i) &1580 + eplus * work2dts(j,i+1)1581 1582 uspr(j,iif) = eminus * work2dus(j,i-1) &1583 + edot * work2dus(j,i) &1584 + eplus * work2dus(j,i+1)1585 1586 END DO1587 1588 END DO1589 END DO1590 1591 !1592 !-- Interpolation in y-direction (quadratic, Clark and Farley)1593 1594 DO j = bdims_rem(2,1), bdims_rem(2,2)1595 1596 bottomy = (nyf+1)/(nyc+1) * j1597 topy = (nyf+1)/(nyc+1) * (j+1) - 11598 1599 DO iif = nxl, nxr1600 DO jjf = bottomy, topy1601 1602 eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc1603 alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.01604 eminus = eps * ( eps - 1.0 ) / 2.0 + alpha1605 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha1606 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha1607 1518 ! SUBROUTINE interpolate_to_fine_flux 1519 ! 1520 ! 1521 ! USE arrays_3d 1522 ! USE control_parameters 1523 ! USE grid_variables 1524 ! USE indices 1525 ! USE pegrid 1526 ! 1527 ! 1528 ! IMPLICIT NONE 1529 ! 1530 ! INTEGER(iwp) :: i 1531 ! INTEGER(iwp) :: j 1532 ! INTEGER(iwp) :: iif 1533 ! INTEGER(iwp) :: jjf 1534 ! INTEGER(iwp) :: bottomx 1535 ! INTEGER(iwp) :: bottomy 1536 ! INTEGER(iwp) :: topx 1537 ! INTEGER(iwp) :: topy 1538 ! REAL(wp) :: eps 1539 ! REAL(wp) :: alpha 1540 ! REAL(wp) :: eminus 1541 ! REAL(wp) :: edot 1542 ! REAL(wp) :: eplus 1543 ! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uswspr 1544 ! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vswspr 1545 ! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tspr 1546 ! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uspr 1547 ! 1548 ! ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 1549 ! ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 1550 ! ALLOCATE( tspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 1551 ! ALLOCATE( uspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 1552 ! 1553 ! ! 1554 ! !-- Initialisation of scalar variables (2D) 1555 ! 1556 ! ! 1557 ! !-- Interpolation in x-direction (quadratic, Clark and Farley) 1558 ! 1559 ! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1 1560 ! DO i = bdims_rem(1,1), bdims_rem(1,2) 1561 ! 1562 ! bottomx = (nxf+1)/(nxc+1) * i 1563 ! topx = (nxf+1)/(nxc+1) * (i+1) - 1 1564 ! 1565 ! DO iif = bottomx, topx 1566 ! 1567 ! eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc 1568 ! alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0 1569 ! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha 1570 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 1571 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 1572 ! 1573 ! uswspr(j,iif) = eminus * work2dusws(j,i-1) & 1574 ! + edot * work2dusws(j,i) & 1575 ! + eplus * work2dusws(j,i+1) 1576 ! 1577 ! vswspr(j,iif) = eminus * work2dvsws(j,i-1) & 1578 ! + edot * work2dvsws(j,i) & 1579 ! + eplus * work2dvsws(j,i+1) 1580 ! 1581 ! tspr(j,iif) = eminus * work2dts(j,i-1) & 1582 ! + edot * work2dts(j,i) & 1583 ! + eplus * work2dts(j,i+1) 1584 ! 1585 ! uspr(j,iif) = eminus * work2dus(j,i-1) & 1586 ! + edot * work2dus(j,i) & 1587 ! + eplus * work2dus(j,i+1) 1588 ! 1589 ! END DO 1590 ! 1591 ! END DO 1592 ! END DO 1593 ! 1594 ! ! 1595 ! !-- Interpolation in y-direction (quadratic, Clark and Farley) 1596 ! 1597 ! DO j = bdims_rem(2,1), bdims_rem(2,2) 1598 ! 1599 ! bottomy = (nyf+1)/(nyc+1) * j 1600 ! topy = (nyf+1)/(nyc+1) * (j+1) - 1 1601 ! 1602 ! DO iif = nxl, nxr 1603 ! DO jjf = bottomy, topy 1604 ! 1605 ! eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc 1606 ! alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0 1607 ! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha 1608 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 1609 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 1610 ! 1608 1611 !! 1609 1612 !!-- TODO … … 1625 1628 ! + edot * uspr(j,iif) & 1626 1629 ! + eplus * uspr(j+1,iif) 1627 1628 END DO1629 END DO1630 1631 END DO1632 1633 1634 DEALLOCATE( uswspr, vswspr )1635 DEALLOCATE( tspr, uspr )1636 1637 1638 END SUBROUTINE interpolate_to_fine_flux1630 ! 1631 ! END DO 1632 ! END DO 1633 ! 1634 ! END DO 1635 ! 1636 ! 1637 ! DEALLOCATE( uswspr, vswspr ) 1638 ! DEALLOCATE( tspr, uspr ) 1639 ! 1640 ! 1641 ! END SUBROUTINE interpolate_to_fine_flux 1639 1642 1640 1643 … … 2398 2401 2399 2402 2400 CONTAINS2401 2402 SUBROUTINE vnest_set_topbc_kh2403 2404 2405 USE arrays_3d2406 USE control_parameters2407 USE grid_variables2408 USE indices2409 USE pegrid2410 2411 2412 IMPLICIT NONE2413 2414 INTEGER(iwp) :: i2415 INTEGER(iwp) :: j2416 INTEGER(iwp) :: k2417 INTEGER(iwp) :: iif2418 INTEGER(iwp) :: jjf2419 INTEGER(iwp) :: bottomx2420 INTEGER(iwp) :: bottomy2421 INTEGER(iwp) :: topx2422 INTEGER(iwp) :: topy2423 REAL(wp) :: eps2424 REAL(wp) :: alpha2425 REAL(wp) :: eminus2426 REAL(wp) :: edot2427 REAL(wp) :: eplus2428 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf2429 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs2430 2431 2432 2433 ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )2434 ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )2435 2436 !2437 !-- Determination of a boundary condition for the potential temperature pt:2438 !-- The scheme derived by Clark and Farley can be used in all three dimensions.2439 2440 !2441 !-- Interpolation in x-direction2442 2443 DO k = bdims_rem(3,1), bdims_rem(3,2)2444 2445 DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+12446 2447 DO i = bdims_rem(1,1), bdims_rem(1,2)2448 2449 bottomx = (nxf+1)/(nxc+1) * i2450 topx = (nxf+1)/(nxc+1) *(i+1) - 12451 2452 DO iif = bottomx, topx2453 2454 eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc2455 2456 alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.02457 2458 eminus = eps * (eps - 1.0 ) / 2.0 + alpha2459 2460 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2461 2462 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2463 2464 ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &2465 + edot * work3d(k,j,i) &2466 + eplus * work3d(k,j,i+1)2467 END DO2468 2469 END DO2470 2471 END DO2472 2473 END DO2474 2475 !2476 !-- Interpolation in y-direction2477 2478 DO k = bdims_rem(3,1), bdims_rem(3,2)2479 2480 DO j = bdims_rem(2,1), bdims_rem(2,2)2481 2482 bottomy = (nyf+1)/(nyc+1) * j2483 topy = (nyf+1)/(nyc+1) * (j+1) - 12484 2485 DO iif = nxl, nxr2486 2487 DO jjf = bottomy, topy2488 2489 eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc2490 2491 alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.02492 2493 eminus = eps * (eps - 1.0) / 2.0 + alpha2494 2495 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2496 2497 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2498 2499 ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &2500 + edot * ptprf(k,j,iif) &2501 + eplus * ptprf(k,j+1,iif)2502 END DO2503 2504 END DO2505 2506 END DO2507 2508 END DO2509 2510 !2511 !-- Interpolation in z-direction2512 2513 DO jjf = nys, nyn2514 DO iif = nxl, nxr2515 2516 eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc2517 2518 alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.02519 2520 eminus = eps * ( eps - 1.0 ) / 2.0 + alpha2521 2522 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2523 2524 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2525 2526 kh (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) &2527 + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) &2528 + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif)2529 2530 END DO2531 END DO2532 2533 DEALLOCATE ( ptprf, ptprs )2534 2535 2536 2537 END SUBROUTINE vnest_set_topbc_kh2538 2539 SUBROUTINE vnest_set_topbc_km2540 2541 2542 USE arrays_3d2543 USE control_parameters2544 USE grid_variables2545 USE indices2546 USE pegrid2547 2548 2549 IMPLICIT NONE2550 2551 INTEGER(iwp) :: i2552 INTEGER(iwp) :: j2553 INTEGER(iwp) :: k2554 INTEGER(iwp) :: iif2555 INTEGER(iwp) :: jjf2556 INTEGER(iwp) :: bottomx2557 INTEGER(iwp) :: bottomy2558 INTEGER(iwp) :: topx2559 INTEGER(iwp) :: topy2560 REAL(wp) :: eps2561 REAL(wp) :: alpha2562 REAL(wp) :: eminus2563 REAL(wp) :: edot2564 REAL(wp) :: eplus2565 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf2566 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs2567 2568 2569 2570 ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )2571 ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )2572 2573 !2574 !-- Determination of a boundary condition for the potential temperature pt:2575 !-- The scheme derived by Clark and Farley can be used in all three dimensions.2576 2577 !2578 !-- Interpolation in x-direction2579 2580 DO k = bdims_rem(3,1), bdims_rem(3,2)2581 2582 DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+12583 2584 DO i = bdims_rem(1,1), bdims_rem(1,2)2585 2586 bottomx = (nxf+1)/(nxc+1) * i2587 topx = (nxf+1)/(nxc+1) *(i+1) - 12588 2589 DO iif = bottomx, topx2590 2591 eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc2592 2593 alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.02594 2595 eminus = eps * (eps - 1.0 ) / 2.0 + alpha2596 2597 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2598 2599 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2600 2601 ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &2602 + edot * work3d(k,j,i) &2603 + eplus * work3d(k,j,i+1)2604 END DO2605 2606 END DO2607 2608 END DO2609 2610 END DO2611 2612 !2613 !-- Interpolation in y-direction2614 2615 DO k = bdims_rem(3,1), bdims_rem(3,2)2616 2617 DO j = bdims_rem(2,1), bdims_rem(2,2)2618 2619 bottomy = (nyf+1)/(nyc+1) * j2620 topy = (nyf+1)/(nyc+1) * (j+1) - 12621 2622 DO iif = nxl, nxr2623 2624 DO jjf = bottomy, topy2625 2626 eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc2627 2628 alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.02629 2630 eminus = eps * (eps - 1.0) / 2.0 + alpha2631 2632 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2633 2634 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2635 2636 ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &2637 + edot * ptprf(k,j,iif) &2638 + eplus * ptprf(k,j+1,iif)2639 END DO2640 2641 END DO2642 2643 END DO2644 2645 END DO2646 2647 !2648 !-- Interpolation in z-direction2649 2650 DO jjf = nys, nyn2651 DO iif = nxl, nxr2652 2653 eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc2654 2655 alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.02656 2657 eminus = eps * ( eps - 1.0 ) / 2.0 + alpha2658 2659 edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha2660 2661 eplus = eps * ( eps + 1.0 ) / 2.0 + alpha2662 2663 km (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) &2664 + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) &2665 + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif)2666 2667 END DO2668 END DO2669 2670 DEALLOCATE ( ptprf, ptprs )2671 2672 2673 2674 END SUBROUTINE vnest_set_topbc_km2403 ! CONTAINS 2404 ! 2405 ! SUBROUTINE vnest_set_topbc_kh 2406 ! 2407 ! 2408 ! USE arrays_3d 2409 ! USE control_parameters 2410 ! USE grid_variables 2411 ! USE indices 2412 ! USE pegrid 2413 ! 2414 ! 2415 ! IMPLICIT NONE 2416 ! 2417 ! INTEGER(iwp) :: i 2418 ! INTEGER(iwp) :: j 2419 ! INTEGER(iwp) :: k 2420 ! INTEGER(iwp) :: iif 2421 ! INTEGER(iwp) :: jjf 2422 ! INTEGER(iwp) :: bottomx 2423 ! INTEGER(iwp) :: bottomy 2424 ! INTEGER(iwp) :: topx 2425 ! INTEGER(iwp) :: topy 2426 ! REAL(wp) :: eps 2427 ! REAL(wp) :: alpha 2428 ! REAL(wp) :: eminus 2429 ! REAL(wp) :: edot 2430 ! REAL(wp) :: eplus 2431 ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf 2432 ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs 2433 ! 2434 ! 2435 ! 2436 ! ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 2437 ! ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) ) 2438 ! 2439 ! ! 2440 ! !-- Determination of a boundary condition for the potential temperature pt: 2441 ! !-- The scheme derived by Clark and Farley can be used in all three dimensions. 2442 ! 2443 ! ! 2444 ! !-- Interpolation in x-direction 2445 ! 2446 ! DO k = bdims_rem(3,1), bdims_rem(3,2) 2447 ! 2448 ! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1 2449 ! 2450 ! DO i = bdims_rem(1,1), bdims_rem(1,2) 2451 ! 2452 ! bottomx = (nxf+1)/(nxc+1) * i 2453 ! topx = (nxf+1)/(nxc+1) *(i+1) - 1 2454 ! 2455 ! DO iif = bottomx, topx 2456 ! 2457 ! eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc 2458 ! 2459 ! alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0 2460 ! 2461 ! eminus = eps * (eps - 1.0 ) / 2.0 + alpha 2462 ! 2463 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2464 ! 2465 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2466 ! 2467 ! ptprf(k,j,iif) = eminus * work3d(k,j,i-1) & 2468 ! + edot * work3d(k,j,i) & 2469 ! + eplus * work3d(k,j,i+1) 2470 ! END DO 2471 ! 2472 ! END DO 2473 ! 2474 ! END DO 2475 ! 2476 ! END DO 2477 ! 2478 ! ! 2479 ! !-- Interpolation in y-direction 2480 ! 2481 ! DO k = bdims_rem(3,1), bdims_rem(3,2) 2482 ! 2483 ! DO j = bdims_rem(2,1), bdims_rem(2,2) 2484 ! 2485 ! bottomy = (nyf+1)/(nyc+1) * j 2486 ! topy = (nyf+1)/(nyc+1) * (j+1) - 1 2487 ! 2488 ! DO iif = nxl, nxr 2489 ! 2490 ! DO jjf = bottomy, topy 2491 ! 2492 ! eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc 2493 ! 2494 ! alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0 2495 ! 2496 ! eminus = eps * (eps - 1.0) / 2.0 + alpha 2497 ! 2498 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2499 ! 2500 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2501 ! 2502 ! ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) & 2503 ! + edot * ptprf(k,j,iif) & 2504 ! + eplus * ptprf(k,j+1,iif) 2505 ! END DO 2506 ! 2507 ! END DO 2508 ! 2509 ! END DO 2510 ! 2511 ! END DO 2512 ! 2513 ! ! 2514 ! !-- Interpolation in z-direction 2515 ! 2516 ! DO jjf = nys, nyn 2517 ! DO iif = nxl, nxr 2518 ! 2519 ! eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc 2520 ! 2521 ! alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0 2522 ! 2523 ! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha 2524 ! 2525 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2526 ! 2527 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2528 ! 2529 ! kh (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) & 2530 ! + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) & 2531 ! + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif) 2532 ! 2533 ! END DO 2534 ! END DO 2535 ! 2536 ! DEALLOCATE ( ptprf, ptprs ) 2537 ! 2538 ! 2539 ! 2540 ! END SUBROUTINE vnest_set_topbc_kh 2541 2542 ! SUBROUTINE vnest_set_topbc_km 2543 ! 2544 ! 2545 ! USE arrays_3d 2546 ! USE control_parameters 2547 ! USE grid_variables 2548 ! USE indices 2549 ! USE pegrid 2550 ! 2551 ! 2552 ! IMPLICIT NONE 2553 ! 2554 ! INTEGER(iwp) :: i 2555 ! INTEGER(iwp) :: j 2556 ! INTEGER(iwp) :: k 2557 ! INTEGER(iwp) :: iif 2558 ! INTEGER(iwp) :: jjf 2559 ! INTEGER(iwp) :: bottomx 2560 ! INTEGER(iwp) :: bottomy 2561 ! INTEGER(iwp) :: topx 2562 ! INTEGER(iwp) :: topy 2563 ! REAL(wp) :: eps 2564 ! REAL(wp) :: alpha 2565 ! REAL(wp) :: eminus 2566 ! REAL(wp) :: edot 2567 ! REAL(wp) :: eplus 2568 ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf 2569 ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs 2570 ! 2571 ! 2572 ! 2573 ! ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) ) 2574 ! ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) ) 2575 ! 2576 ! ! 2577 ! !-- Determination of a boundary condition for the potential temperature pt: 2578 ! !-- The scheme derived by Clark and Farley can be used in all three dimensions. 2579 ! 2580 ! ! 2581 ! !-- Interpolation in x-direction 2582 ! 2583 ! DO k = bdims_rem(3,1), bdims_rem(3,2) 2584 ! 2585 ! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1 2586 ! 2587 ! DO i = bdims_rem(1,1), bdims_rem(1,2) 2588 ! 2589 ! bottomx = (nxf+1)/(nxc+1) * i 2590 ! topx = (nxf+1)/(nxc+1) *(i+1) - 1 2591 ! 2592 ! DO iif = bottomx, topx 2593 ! 2594 ! eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc 2595 ! 2596 ! alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0 2597 ! 2598 ! eminus = eps * (eps - 1.0 ) / 2.0 + alpha 2599 ! 2600 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2601 ! 2602 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2603 ! 2604 ! ptprf(k,j,iif) = eminus * work3d(k,j,i-1) & 2605 ! + edot * work3d(k,j,i) & 2606 ! + eplus * work3d(k,j,i+1) 2607 ! END DO 2608 ! 2609 ! END DO 2610 ! 2611 ! END DO 2612 ! 2613 ! END DO 2614 ! 2615 ! ! 2616 ! !-- Interpolation in y-direction 2617 ! 2618 ! DO k = bdims_rem(3,1), bdims_rem(3,2) 2619 ! 2620 ! DO j = bdims_rem(2,1), bdims_rem(2,2) 2621 ! 2622 ! bottomy = (nyf+1)/(nyc+1) * j 2623 ! topy = (nyf+1)/(nyc+1) * (j+1) - 1 2624 ! 2625 ! DO iif = nxl, nxr 2626 ! 2627 ! DO jjf = bottomy, topy 2628 ! 2629 ! eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc 2630 ! 2631 ! alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0 2632 ! 2633 ! eminus = eps * (eps - 1.0) / 2.0 + alpha 2634 ! 2635 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2636 ! 2637 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2638 ! 2639 ! ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) & 2640 ! + edot * ptprf(k,j,iif) & 2641 ! + eplus * ptprf(k,j+1,iif) 2642 ! END DO 2643 ! 2644 ! END DO 2645 ! 2646 ! END DO 2647 ! 2648 ! END DO 2649 ! 2650 ! ! 2651 ! !-- Interpolation in z-direction 2652 ! 2653 ! DO jjf = nys, nyn 2654 ! DO iif = nxl, nxr 2655 ! 2656 ! eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc 2657 ! 2658 ! alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0 2659 ! 2660 ! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha 2661 ! 2662 ! edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha 2663 ! 2664 ! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha 2665 ! 2666 ! km (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) & 2667 ! + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) & 2668 ! + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif) 2669 ! 2670 ! END DO 2671 ! END DO 2672 ! 2673 ! DEALLOCATE ( ptprf, ptprs ) 2674 ! 2675 ! 2676 ! 2677 ! END SUBROUTINE vnest_set_topbc_km 2675 2678 2676 2679
Note: See TracChangeset
for help on using the changeset viewer.
