Ignore:
Timestamp:
Sep 13, 2012 2:08:46 PM (9 years ago)
Author:
raasch
Message:

leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code,
reading of dt_fixed from restart file removed

File:
1 edited

Legend:

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

    r667 r1001  
    44! Current revisions:
    55! -----------------
    6 !
     6! all actions concerning upstream-spline-method removed
    77!
    88! Former revisions:
     
    8080    ENDIF
    8181
    82     IF ( momentum_advec == 'ups-scheme' .OR. scalar_advec  == 'ups-scheme' )  &
    83     THEN
    84 
    85 !
    86 !--    Provide the constant parameters for the Upstream-Spline advection scheme.
    87 !--    In x- und y-direction the Sherman-Morrison formula is applied
    88 !--    (cf. Press et al, 1986 (Numerical Recipes)).
    89 !
    90 !--    Allocate nonlocal arrays
    91        ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), &
    92                  spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1),      &
    93                  spl_tri_zw(5,nzb:nzt+1) )
    94 
    95 !
    96 !--    Provide diagonal elements of the tridiagonal matrices for all
    97 !--    directions
    98        spl_tri_x(1,:)  = 2.0
    99        spl_tri_y(1,:)  = 2.0
    100        spl_tri_zu(1,:) = 2.0
    101        spl_tri_zw(1,:) = 2.0
    102 
    103 !
    104 !--    Elements of the cyclic tridiagonal matrix
    105 !--    (same for all horizontal directions)
    106        spl_alpha = 0.5
    107        spl_beta  = 0.5
    108 
    109 !
    110 !--    Sub- and superdiagonal elements, x-direction
    111        spl_tri_x(2,0:nx) = 0.5
    112        spl_tri_x(3,0:nx) = 0.5
    113 
    114 !
    115 !--    mMdify the diagonal elements (Sherman-Morrison)
    116        spl_gamma_x    = -spl_tri_x(1,0)
    117        spl_tri_x(1,0)  = spl_tri_x(1,0) - spl_gamma_x
    118        spl_tri_x(1,nx) = spl_tri_x(1,nx) - spl_alpha * spl_beta / spl_gamma_x
    119 
    120 !
    121 !--    Split the tridiagonal matrix for Thomas algorithm
    122        spl_tri_x(4,0) = spl_tri_x(1,0)
    123        DO  i = 1, nx
    124           spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i-1)
    125           spl_tri_x(4,i) = spl_tri_x(1,i) - spl_tri_x(5,i) * spl_tri_x(3,i-1)
    126        ENDDO
    127 
    128 !
    129 !--    Allocate arrays required locally
    130        ALLOCATE( temp(0:nx), spl_u(0:nx) )
    131 
    132 !
    133 !--    Provide "corrective vector", x-direction
    134        spl_u(0)      = spl_gamma_x
    135        spl_u(1:nx-1) = 0.0
    136        spl_u(nx)     = spl_alpha
    137 
    138 !
    139 !--    Solve the system of equations for the corrective vector
    140 !--    (Sherman-Morrison)
    141        temp(0) = spl_u(0)
    142        DO  i = 1, nx
    143           temp(i) = spl_u(i) - spl_tri_x(5,i) * temp(i-1)
    144        ENDDO
    145        spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx)
    146        DO  i = nx-1, 0, -1
    147           spl_z_x(i) = ( temp(i) - spl_tri_x(3,i) * spl_z_x(i+1) ) / &
    148                        spl_tri_x(4,i)
    149        ENDDO
    150 
    151 !
    152 !--    Deallocate local arrays, for they are allocated in a different way for
    153 !--    operations in y-direction
    154        DEALLOCATE( temp, spl_u )   
    155    
    156 !
    157 !--    Provide sub- and superdiagonal elements, y-direction
    158        spl_tri_y(2,0:ny) = 0.5
    159        spl_tri_y(3,0:ny) = 0.5
    160 
    161 !
    162 !--    Modify the diagonal elements (Sherman-Morrison)
    163        spl_gamma_y    = -spl_tri_y(1,0)
    164        spl_tri_y(1,0)  = spl_tri_y(1,0) - spl_gamma_y
    165        spl_tri_y(1,ny) = spl_tri_y(1,ny) - spl_alpha * spl_beta / spl_gamma_y
    166 
    167 !
    168 !--    Split the tridiagonal matrix for Thomas algorithm
    169        spl_tri_y(4,0) = spl_tri_y(1,0)
    170        DO  j = 1, ny
    171           spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j-1)
    172           spl_tri_y(4,j) = spl_tri_y(1,j) - spl_tri_y(5,j) * spl_tri_y(3,j-1)
    173        ENDDO
    174 
    175 !
    176 !--    Allocate arrays required locally
    177        ALLOCATE( temp(0:ny), spl_u(0:ny) )
    178 
    179 !
    180 !--    Provide "corrective vector", y-direction
    181        spl_u(0)      = spl_gamma_y
    182        spl_u(1:ny-1) = 0.0
    183        spl_u(ny)     = spl_alpha
    184    
    185 !
    186 !--    Solve the system of equations for the corrective vector
    187 !--    (Sherman-Morrison)
    188        temp = 0.0
    189        spl_z_y = 0.0
    190        temp(0) = spl_u(0)
    191        DO  j = 1, ny
    192           temp(j) = spl_u(j) - spl_tri_y(5,j) * temp(j-1)
    193        ENDDO
    194        spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny)
    195        DO  j = ny-1, 0, -1
    196           spl_z_y(j) = ( temp(j) - spl_tri_y(3,j) * spl_z_y(j+1) ) / &
    197                        spl_tri_y(4,j)
    198        ENDDO
    199 
    200 !
    201 !--    deallocate local arrays, for they are no longer required
    202        DEALLOCATE( temp, spl_u )
    203 
    204 !
    205 !--    provide sub- and superdiagonal elements, z-direction
    206        spl_tri_zu(2,nzb)   = 0.0
    207        spl_tri_zu(2,nzt+1) = 1.0
    208        spl_tri_zw(2,nzb)   = 0.0
    209        spl_tri_zw(2,nzt+1) = 1.0
    210 
    211        spl_tri_zu(3,nzb)   = 1.0
    212        spl_tri_zu(3,nzt+1) = 0.0
    213        spl_tri_zw(3,nzb)   = 1.0
    214        spl_tri_zw(3,nzt+1) = 0.0
    215 
    216        DO  k = nzb+1, nzt
    217           spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) )
    218           spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) )
    219           spl_tri_zu(3,k) = 1.0 - spl_tri_zu(2,k)
    220           spl_tri_zw(3,k) = 1.0 - spl_tri_zw(2,k)
    221        ENDDO
    222    
    223        spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb)
    224        spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb)
    225        DO  k = nzb+1, nzt+1
    226           spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k-1)
    227           spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k-1)
    228           spl_tri_zu(4,k) = spl_tri_zu(1,k) - spl_tri_zu(5,k) * &
    229                             spl_tri_zu(3,k-1)
    230           spl_tri_zw(4,k) = spl_tri_zw(1,k) - spl_tri_zw(5,k) * &
    231                             spl_tri_zw(3,k-1)
    232        ENDDO
    233 
    234     ENDIF
    235 
    23682 END SUBROUTINE init_advec
Note: See TracChangeset for help on using the changeset viewer.