Ignore:
Timestamp:
Aug 9, 2012 8:28:32 AM (9 years ago)
Author:
fricke
Message:

merge fricke branch back into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/diffusion_w.f90

    r668 r978  
    44! Current revisions:
    55! -----------------
     6! outflow damping layer removed
     7! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
     8! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
    69!
    710! Former revisions:
     
    5457! Call for all grid points
    5558!------------------------------------------------------------------------------!
    56     SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
    57                             w )
     59    SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w )
    5860
    5961       USE control_parameters
     
    6466
    6567       INTEGER ::  i, j, k
    66        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    67                    kmyp_z
    68        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
    69                    km_damp_y(nysg:nyng)
     68       REAL    ::  kmxm, kmxp, kmym, kmyp
     69       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    7070       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7171       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     
    8888!
    8989!--             Interpolate eddy diffusivities on staggered gridpoints
    90                 kmxp_x = 0.25 * &
    91                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    92                 kmxm_x = 0.25 * &
    93                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    94                 kmxp_z = kmxp_x
    95                 kmxm_z = kmxm_x
    96                 kmyp_y = 0.25 * &
    97                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    98                 kmym_y = 0.25 * &
    99                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    100                 kmyp_z = kmyp_y
    101                 kmym_z = kmym_y
    102 !
    103 !--             Increase diffusion at the outflow boundary in case of
    104 !--             non-cyclic lateral boundaries. Damping is only needed for
    105 !--             velocity components parallel to the outflow boundary in
    106 !--             the direction normal to the outflow boundary.
    107                 IF ( .NOT. bc_lr_cyc )  THEN
    108                    kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    109                    kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    110                 ENDIF
    111                 IF ( .NOT. bc_ns_cyc )  THEN
    112                    kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    113                    kmym_y = MAX( kmym_y, km_damp_y(j) )
    114                 ENDIF
     90                kmxp = 0.25 * &
     91                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     92                kmxm = 0.25 * &
     93                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     94                kmyp = 0.25 * &
     95                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     96                kmym = 0.25 * &
     97                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    11598
    11699                tend(k,j,i) = tend(k,j,i)                                      &
    117                       & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    118                       &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
    119                       &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    120                       &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     100                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     101                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
     102                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
     103                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    121104                      &   ) * ddx                                              &
    122                       & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    123                       &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
    124                       &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    125                       &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     105                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     106                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
     107                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
     108                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    126109                      &   ) * ddy                                              &
    127110                      & + 2.0 * (                                              &
     
    138121!
    139122!--                Interpolate eddy diffusivities on staggered gridpoints
    140                    kmxp_x = 0.25 * &
    141                             ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    142                    kmxm_x = 0.25 * &
    143                             ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    144                    kmxp_z = kmxp_x
    145                    kmxm_z = kmxm_x
    146                    kmyp_y = 0.25 * &
    147                             ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    148                    kmym_y = 0.25 * &
    149                             ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    150                    kmyp_z = kmyp_y
    151                    kmym_z = kmym_y
    152 !
    153 !--                Increase diffusion at the outflow boundary in case of
    154 !--                non-cyclic lateral boundaries. Damping is only needed for
    155 !--                velocity components parallel to the outflow boundary in
    156 !--                the direction normal to the outflow boundary.
    157                    IF ( .NOT. bc_lr_cyc )  THEN
    158                       kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    159                       kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    160                    ENDIF
    161                    IF ( .NOT. bc_ns_cyc )  THEN
    162                       kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    163                       kmym_y = MAX( kmym_y, km_damp_y(j) )
    164                    ENDIF
     123                   kmxp = 0.25 * &
     124                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     125                   kmxm = 0.25 * &
     126                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     127                   kmyp = 0.25 * &
     128                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     129                   kmym = 0.25 * &
     130                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    165131
    166132                   tend(k,j,i) = tend(k,j,i)                                   &
    167133                                 + (   fwxp(j,i) * (                           &
    168                             kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    169                           + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
     134                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     135                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
    170136                                                   )                           &
    171137                                     - fwxm(j,i) * (                           &
    172                             kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
    173                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
     138                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     139                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)  &
    174140                                                   )                           &
    175141                                     + wall_w_x(j,i) * wsus(k,j,i)             &
    176142                                   ) * ddx                                     &
    177143                                 + (   fwyp(j,i) * (                           &
    178                             kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    179                           + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
     144                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     145                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
    180146                                                   )                           &
    181147                                     - fwym(j,i) * (                           &
    182                             kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
    183                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
     148                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     149                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)  &
    184150                                                   )                           &
    185151                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
     
    201167! Call for grid point i,j
    202168!------------------------------------------------------------------------------!
    203     SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    204                                tend, u, v, w )
     169    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w )
    205170
    206171       USE control_parameters
     
    211176
    212177       INTEGER ::  i, j, k
    213        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    214                    kmyp_z
    215        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
    216                    km_damp_y(nysg:nyng)
     178       REAL    ::  kmxm, kmxp, kmym, kmyp
     179       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    217180       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    218181       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
     
    223186!
    224187!--       Interpolate eddy diffusivities on staggered gridpoints
    225           kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    226           kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    227           kmxp_z = kmxp_x
    228           kmxm_z = kmxm_x
    229           kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    230           kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    231           kmyp_z = kmyp_y
    232           kmym_z = kmym_y
    233 !
    234 !--       Increase diffusion at the outflow boundary in case of non-cyclic
    235 !--       lateral boundaries. Damping is only needed for velocity components
    236 !--       parallel to the outflow boundary in the direction normal to the
    237 !--       outflow boundary.
    238           IF ( .NOT. bc_lr_cyc )  THEN
    239              kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    240              kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    241           ENDIF
    242           IF ( .NOT. bc_ns_cyc )  THEN
    243              kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    244              kmym_y = MAX( kmym_y, km_damp_y(j) )
    245           ENDIF
     188          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     189          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     190          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     191          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    246192
    247193          tend(k,j,i) = tend(k,j,i)                                            &
    248                       & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    249                       &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
    250                       &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    251                       &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     194                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     195                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
     196                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
     197                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    252198                      &   ) * ddx                                              &
    253                       & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    254                       &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
    255                       &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    256                       &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     199                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     200                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
     201                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
     202                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    257203                      &   ) * ddy                                              &
    258204                      & + 2.0 * (                                              &
     
    285231!
    286232!--          Interpolate eddy diffusivities on staggered gridpoints
    287              kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    288              kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    289              kmxp_z = kmxp_x
    290              kmxm_z = kmxm_x
    291              kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    292              kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    293              kmyp_z = kmyp_y
    294              kmym_z = kmym_y
    295 !
    296 !--          Increase diffusion at the outflow boundary in case of
    297 !--          non-cyclic lateral boundaries. Damping is only needed for
    298 !--          velocity components parallel to the outflow boundary in
    299 !--          the direction normal to the outflow boundary.
    300              IF ( .NOT. bc_lr_cyc )  THEN
    301                 kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    302                 kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    303              ENDIF
    304              IF ( .NOT. bc_ns_cyc )  THEN
    305                 kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    306                 kmym_y = MAX( kmym_y, km_damp_y(j) )
    307              ENDIF
     233             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     234             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     235             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     236             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    308237
    309238             tend(k,j,i) = tend(k,j,i)                                         &
    310239                                 + (   fwxp(j,i) * (                           &
    311                             kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    312                           + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
     240                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     241                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
    313242                                                   )                           &
    314243                                     - fwxm(j,i) * (                           &
    315                             kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
    316                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
     244                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     245                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)  &
    317246                                                   )                           &
    318247                                     + wall_w_x(j,i) * wsus(k)                 &
    319248                                   ) * ddx                                     &
    320249                                 + (   fwyp(j,i) * (                           &
    321                             kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    322                           + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
     250                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     251                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
    323252                                                   )                           &
    324253                                     - fwym(j,i) * (                           &
    325                             kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
    326                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
     254                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     255                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)  &
    327256                                                   )                           &
    328257                                     + wall_w_y(j,i) * wsvs(k)                 &
Note: See TracChangeset for help on using the changeset viewer.