Changeset 3802


Ignore:
Timestamp:
Mar 17, 2019 1:33:42 PM (3 years ago)
Author:
raasch
Message:

unused variables removed, unused subroutines commented out, type conversion added to avoid compiler warning about constant integer division truncation, script document_changes made bash compatible

Location:
palm/trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/document_changes

    r3665 r3802  
    1 #!/bin/ksh
     1#!/bin/bash
    22#------------------------------------------------------------------------------#
    33# This file is part of the PALM model system.
     
    2525# -----------------
    2626# $Id$
     27# script made bash compatible
     28#
     29# 3665 2019-01-10 08:28:24Z raasch
    2730# Corrected "Former revisions" section
    2831#
     
    117120                           -and -not -name '*.eps'  -and -not -name '*.png'    \
    118121                           -and -not -name '*.svn*' -and -not -name '*~'       \
    119                            -and -not -name '*.tar'); do
     122                           -and -not -name '*.tar'  -and -not -type d ); do
    120123
    121124#
     
    198201            timestamp_string="$comment_char $timestamp"
    199202          fi
    200 
    201          
    202          
    203203         
    204204       done <"$fn"
  • palm/trunk/SOURCE/random_function_mod.f90

    r3655 r3802  
    2525! -----------------
    2626! $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
    2731! Corrected "Former revisions" section
    2832!
     
    131135
    132136       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 )
    134138
    135139       IF ( idum .le. 0  .or.  random_iy .eq. 0 )  THEN
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3769 r3802  
    2828! -----------------
    2929! $Id$
     30! unused subroutine commented out
     31!
     32! 3769 2019-02-28 10:16:49Z moh.hefny
    3033! removed unused variables
    3134!
     
    86298632 
    86308633 
    8631      CONTAINS
    8632  !------------------------------------------------------------------------------!
    8633  ! Description:
    8634  ! ------------
    8635  !> Calculation of specific humidity of the skin layer (surface). It is assumend
    8636  !> that the skin is always saturated.
    8637  !------------------------------------------------------------------------------!
    8638         SUBROUTINE calc_q_surface_usm
    8639 
    8640            IMPLICIT NONE
    8641 
    8642            REAL(wp) :: resistance    !< aerodynamic and soil resistance term
    8643 
    8644            DO  m = 1, surf_usm_h%ns
    8645 
    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 saturation
    8652               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 saturation
    8659               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 zero
    8665 !              IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp )  surf_usm_h%r_a_green(m) = 1.0E-10_wp
    8666 
    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 surface
    8671               IF ( bulk_cloud_model )  THEN
    8672                  q(k,j,i) = resistance * q_s +                   &
    8673                                             ( 1.0_wp - resistance ) *              &
    8674                                             ( q(k,j,i) - ql(k,j,i) )
    8675               ELSE
    8676                  q(k,j,i) = resistance * q_s +                   &
    8677                                             ( 1.0_wp - resistance ) *              &
    8678                                               q(k,j,i)
    8679               ENDIF
    8680 
    8681 !
    8682 !--          Update virtual potential temperature
    8683               vpt(k,j,i) = pt(k,j,i) *         &
    8684                          ( 1.0_wp + 0.61_wp * q(k,j,i) )
    8685 
    8686            ENDDO
    8687 
    8688 !
    8689 !--       Now, treat vertical surface elements
    8690            DO  l = 0, 3
    8691               DO  m = 1, surf_usm_v(l)%ns
    8692 !
    8693 !--             Get indices of respective grid point
    8694                  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 saturation
    8700                  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 saturation
    8707                  q_s = 0.622_wp * e_s / ( surface_pressure -e_s )
    8708 
    8709 !
    8710 !--             Calculate specific humidity at surface
    8711                  IF ( bulk_cloud_model )  THEN
    8712                     q(k,j,i) = ( q(k,j,i) - ql(k,j,i) )
    8713                  ELSE
    8714                     q(k,j,i) = q(k,j,i)
    8715                  ENDIF
    8716 !
    8717 !--             Update virtual potential temperature
    8718                  vpt(k,j,i) = pt(k,j,i) *         &
    8719                             ( 1.0_wp + 0.61_wp * q(k,j,i) )
    8720 
    8721               ENDDO
    8722 
    8723            ENDDO
    8724 
    8725         END SUBROUTINE calc_q_surface_usm
     8634!     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
    87268729       
    87278730     END SUBROUTINE usm_surface_energy_balance
  • palm/trunk/SOURCE/vertical_nesting_mod.f90

    r3655 r3802  
    2626! -----------------
    2727! $Id$
     28! unused subroutines commented out
     29!
     30! 3655 2019-01-07 16:51:22Z knoop
    2831! unused variables removed
    2932!
     
    15131516     
    15141517     
    1515        SUBROUTINE interpolate_to_fine_flux
    1516      
    1517      
    1518            USE arrays_3d
    1519            USE control_parameters
    1520            USE grid_variables
    1521            USE indices
    1522            USE pegrid
    1523            
    1524      
    1525            IMPLICIT NONE
    1526 
    1527            INTEGER(iwp)                       :: i
    1528            INTEGER(iwp)                       :: j
    1529            INTEGER(iwp)                       :: iif
    1530            INTEGER(iwp)                       :: jjf
    1531            INTEGER(iwp)                       :: bottomx
    1532            INTEGER(iwp)                       :: bottomy
    1533            INTEGER(iwp)                       :: topx
    1534            INTEGER(iwp)                       :: topy
    1535            REAL(wp)                           :: eps
    1536            REAL(wp)                           :: alpha
    1537            REAL(wp)                           :: eminus
    1538            REAL(wp)                           :: edot
    1539            REAL(wp)                           :: eplus
    1540            REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uswspr
    1541            REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: vswspr
    1542            REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: tspr
    1543            REAL(wp), DIMENSION(:,:), ALLOCATABLE   :: uspr
    1544      
    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)+1
    1557                DO i = bdims_rem(1,1), bdims_rem(1,2)
    1558              
    1559                    bottomx = (nxf+1)/(nxc+1) * i
    1560                    topx    = (nxf+1)/(nxc+1) * (i+1) - 1
    1561      
    1562                    DO iif = bottomx, topx
    1563      
    1564                        eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
    1565                        alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
    1566                        eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
    1567                        edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    1568                        eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    1569      
    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 DO
    1587      
    1588                END DO
    1589            END DO
    1590      
    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) * j
    1597                topy    = (nyf+1)/(nyc+1) * (j+1) - 1
    1598      
    1599                DO iif = nxl, nxr
    1600                    DO jjf = bottomy, topy
    1601      
    1602                        eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
    1603                        alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
    1604                        eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
    1605                        edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    1606                        eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    1607 
     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!
    16081611!!
    16091612!!--   TODO
     
    16251628!                        + edot  * uspr(j,iif)   &
    16261629!                        + eplus  * uspr(j+1,iif)
    1627      
    1628                    END DO
    1629                END DO
    1630      
    1631            END DO
    1632      
    1633      
    1634            DEALLOCATE( uswspr, vswspr )
    1635            DEALLOCATE( tspr, uspr )
    1636      
    1637      
    1638        END SUBROUTINE interpolate_to_fine_flux
     1630!     
     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
    16391642   
    16401643   
     
    23982401   
    23992402   
    2400     CONTAINS
    2401 
    2402        SUBROUTINE vnest_set_topbc_kh
    2403        
    2404        
    2405            USE arrays_3d
    2406            USE control_parameters
    2407            USE grid_variables
    2408            USE indices
    2409            USE pegrid
    2410            
    2411        
    2412            IMPLICIT NONE
    2413 
    2414            INTEGER(iwp)                       :: i
    2415            INTEGER(iwp)                       :: j
    2416            INTEGER(iwp)                       :: k
    2417            INTEGER(iwp)                       :: iif
    2418            INTEGER(iwp)                       :: jjf
    2419            INTEGER(iwp)                       :: bottomx
    2420            INTEGER(iwp)                       :: bottomy
    2421            INTEGER(iwp)                       :: topx
    2422            INTEGER(iwp)                       :: topy
    2423            REAL(wp)                           :: eps
    2424            REAL(wp)                           :: alpha
    2425            REAL(wp)                           :: eminus
    2426            REAL(wp)                           :: edot
    2427            REAL(wp)                           :: eplus
    2428            REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
    2429            REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
    2430        
    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-direction
    2442        
    2443            DO k = bdims_rem(3,1), bdims_rem(3,2)
    2444        
    2445                DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
    2446        
    2447                    DO i = bdims_rem(1,1), bdims_rem(1,2)
    2448        
    2449                        bottomx = (nxf+1)/(nxc+1) * i
    2450                        topx    = (nxf+1)/(nxc+1) *(i+1) - 1
    2451        
    2452                        DO iif = bottomx, topx
    2453        
    2454                            eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
    2455        
    2456                            alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
    2457        
    2458                            eminus = eps * (eps - 1.0 ) / 2.0 + alpha
    2459        
    2460                            edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2461        
    2462                            eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
    2463        
    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 DO
    2468        
    2469                    END DO
    2470        
    2471                END DO
    2472        
    2473            END DO
    2474        
    2475            !
    2476            !-- Interpolation in y-direction
    2477        
    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) * j
    2483                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
    2484        
    2485                    DO iif = nxl, nxr
    2486        
    2487                        DO jjf = bottomy, topy
    2488        
    2489                            eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
    2490        
    2491                            alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
    2492                    
    2493                            eminus = eps * (eps - 1.0) / 2.0 + alpha
    2494        
    2495                            edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2496        
    2497                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    2498        
    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 DO
    2503        
    2504                    END DO
    2505        
    2506                END DO
    2507        
    2508            END DO
    2509        
    2510            !
    2511            !-- Interpolation in z-direction
    2512        
    2513            DO jjf = nys, nyn
    2514                DO iif = nxl, nxr
    2515        
    2516                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
    2517        
    2518                    alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
    2519        
    2520                    eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
    2521        
    2522                    edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2523        
    2524                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    2525        
    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 DO
    2531            END DO         
    2532            
    2533            DEALLOCATE ( ptprf, ptprs )
    2534        
    2535        
    2536        
    2537        END SUBROUTINE vnest_set_topbc_kh
    2538        
    2539        SUBROUTINE vnest_set_topbc_km
    2540        
    2541        
    2542            USE arrays_3d
    2543            USE control_parameters
    2544            USE grid_variables
    2545            USE indices
    2546            USE pegrid
    2547            
    2548        
    2549            IMPLICIT NONE
    2550 
    2551            INTEGER(iwp)                       :: i
    2552            INTEGER(iwp)                       :: j
    2553            INTEGER(iwp)                       :: k
    2554            INTEGER(iwp)                       :: iif
    2555            INTEGER(iwp)                       :: jjf
    2556            INTEGER(iwp)                       :: bottomx
    2557            INTEGER(iwp)                       :: bottomy
    2558            INTEGER(iwp)                       :: topx
    2559            INTEGER(iwp)                       :: topy
    2560            REAL(wp)                           :: eps
    2561            REAL(wp)                           :: alpha
    2562            REAL(wp)                           :: eminus
    2563            REAL(wp)                           :: edot
    2564            REAL(wp)                           :: eplus
    2565            REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
    2566            REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
    2567        
    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-direction
    2579        
    2580            DO k = bdims_rem(3,1), bdims_rem(3,2)
    2581        
    2582                DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
    2583        
    2584                    DO i = bdims_rem(1,1), bdims_rem(1,2)
    2585        
    2586                        bottomx = (nxf+1)/(nxc+1) * i
    2587                        topx    = (nxf+1)/(nxc+1) *(i+1) - 1
    2588        
    2589                        DO iif = bottomx, topx
    2590        
    2591                            eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
    2592        
    2593                            alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
    2594        
    2595                            eminus = eps * (eps - 1.0 ) / 2.0 + alpha
    2596        
    2597                            edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2598        
    2599                            eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
    2600        
    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 DO
    2605        
    2606                    END DO
    2607        
    2608                END DO
    2609        
    2610            END DO
    2611        
    2612            !
    2613            !-- Interpolation in y-direction
    2614        
    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) * j
    2620                    topy    = (nyf+1)/(nyc+1) * (j+1) - 1
    2621        
    2622                    DO iif = nxl, nxr
    2623        
    2624                        DO jjf = bottomy, topy
    2625        
    2626                            eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
    2627        
    2628                            alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
    2629                    
    2630                            eminus = eps * (eps - 1.0) / 2.0 + alpha
    2631        
    2632                            edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2633        
    2634                            eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    2635        
    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 DO
    2640        
    2641                    END DO
    2642        
    2643                END DO
    2644        
    2645            END DO
    2646        
    2647            !
    2648            !-- Interpolation in z-direction
    2649        
    2650            DO jjf = nys, nyn
    2651                DO iif = nxl, nxr
    2652        
    2653                    eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
    2654        
    2655                    alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
    2656        
    2657                    eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
    2658        
    2659                    edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
    2660        
    2661                    eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
    2662        
    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 DO
    2668            END DO         
    2669            
    2670            DEALLOCATE ( ptprf, ptprs )
    2671        
    2672        
    2673        
    2674        END SUBROUTINE vnest_set_topbc_km
     2403!    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
    26752678   
    26762679 
  • palm/trunk/UTIL/inifor/src/inifor_grid.f90

    r3801 r3802  
    2626! -----------------
    2727! $Id$
     28! unused variable removed
     29!
     30! 3801 2019-03-15 17:14:25Z eckhard
    2831! Read COSMO rotated pole from HHL file
    2932! Check for PALM and averaging domains extending outside COSMO domain
     
    353356!------------------------------------------------------------------------------!
    354357    SUBROUTINE setup_parameters()
    355        INTEGER ::  k !< loop index
    356358
    357359!
Note: See TracChangeset for help on using the changeset viewer.