Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 65


Ignore:
Timestamp:
Mar 13, 2007 12:11:43 PM (18 years ago)
Author:
raasch
Message:

bugfix for NEC because -C hopt gives runtime errors

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/diffusion_e.f90

    r57 r65  
    6464 
    6565
    66        DO  i = nxl, nxr
    67           DO  j = nys, nyn
    68 !
    69 !--          First, calculate phi-function for eventually adjusting the &
    70 !--          mixing length to the prandtl mixing length
    71              IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
    72                 IF ( rif(j,i) >= 0.0 )  THEN
    73                    phi_m = 1.0 + 5.0 * rif(j,i)
    74                 ELSE
    75                    phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
     66!
     67!--    This if clause must be outside the k-loop because otherwise
     68!--    runtime errors occur with -C hopt on NEC
     69       IF ( use_pt_reference )  THEN
     70
     71          DO  i = nxl, nxr
     72             DO  j = nys, nyn
     73!
     74!--             First, calculate phi-function for eventually adjusting the &
     75!--             mixing length to the prandtl mixing length
     76                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
     77                   IF ( rif(j,i) >= 0.0 )  THEN
     78                      phi_m = 1.0 + 5.0 * rif(j,i)
     79                   ELSE
     80                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
     81                   ENDIF
    7682                ENDIF
    77              ENDIF
    78 
    79              DO  k = nzb_s_inner(j,i)+1, nzt
    80 !
    81 !--             Calculate the mixing length (for dissipation)
    82                 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    83                 IF ( dpt_dz > 0.0 ) THEN
    84                    IF ( use_pt_reference )  THEN
     83
     84                DO  k = nzb_s_inner(j,i)+1, nzt
     85!
     86!--                Calculate the mixing length (for dissipation)
     87                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
     88                   IF ( dpt_dz > 0.0 ) THEN
    8589                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    8690                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
    8791                   ELSE
    88                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    89                                         SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
    90                    ENDIF
    91                 ELSE
    92                    l_stable = l_grid(k)
    93                 ENDIF
    94 !
    95 !--             Adjustment of the mixing length
    96                 IF ( wall_adjustment )  THEN
    97                    l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
    98                                   l_stable )
    99                    ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
    100                 ELSE
    101                    l(k,j)  = MIN( l_grid(k), l_stable )
    102                    ll(k,j) = l_grid(k)
    103                 ENDIF
    104                 IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
    105                    l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
    106                    ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
    107                 ENDIF
    108 
    109              ENDDO
    110           ENDDO
    111 !
    112 !--       Calculate the tendency terms
    113           DO  j = nys, nyn
    114              DO  k = nzb_s_inner(j,i)+1, nzt
    115 
    116                  dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
    117                                     e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
    118 
    119                  tend(k,j,i) = tend(k,j,i)                                     &
     92                      l_stable = l_grid(k)
     93                   ENDIF
     94!
     95!--                Adjustment of the mixing length
     96                   IF ( wall_adjustment )  THEN
     97                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
     98                                     l_stable )
     99                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
     100                   ELSE
     101                      l(k,j)  = MIN( l_grid(k), l_stable )
     102                      ll(k,j) = l_grid(k)
     103                   ENDIF
     104                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
     105                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
     106                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
     107                   ENDIF
     108
     109                ENDDO
     110             ENDDO
     111
     112!
     113!--          Calculate the tendency terms
     114             DO  j = nys, nyn
     115                DO  k = nzb_s_inner(j,i)+1, nzt
     116
     117                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
     118                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
     119
     120                    tend(k,j,i) = tend(k,j,i)                                  &
    120121                                        + (                                    &
    121122                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
    … …  
    132133                             - dissipation(k,j)
    133134
    134              ENDDO
     135                ENDDO
     136             ENDDO
     137
     138!
     139!--          Store dissipation if needed for calculating the sgs particle
     140!--          velocities
     141             IF ( use_sgs_for_particles )  THEN
     142                DO  j = nys, nyn
     143                   DO  k = nzb_s_inner(j,i)+1, nzt
     144                      diss(k,j,i) = dissipation(k,j)
     145                   ENDDO
     146                ENDDO
     147             ENDIF
     148
    135149          ENDDO
    136150
    137 !
    138 !--       Store dissipation if needed for calculating the sgs particle
    139 !--       velocities
    140           IF ( use_sgs_for_particles )  THEN
    141              DO  j = nys, nyn
     151       ELSE
     152
     153          DO  i = nxl, nxr
     154             DO  j = nys, nyn
     155!
     156!--             First, calculate phi-function for eventually adjusting the &
     157!--             mixing length to the prandtl mixing length
     158                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
     159                   IF ( rif(j,i) >= 0.0 )  THEN
     160                      phi_m = 1.0 + 5.0 * rif(j,i)
     161                   ELSE
     162                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
     163                   ENDIF
     164                ENDIF
     165
    142166                DO  k = nzb_s_inner(j,i)+1, nzt
    143                    diss(k,j,i) = dissipation(k,j)
    144                 ENDDO
    145              ENDDO
    146           ENDIF
    147 
    148        ENDDO
     167!
     168!--                Calculate the mixing length (for dissipation)
     169                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
     170                   IF ( dpt_dz > 0.0 ) THEN
     171                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
     172                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     173                   ELSE
     174                      l_stable = l_grid(k)
     175                   ENDIF
     176!
     177!--                Adjustment of the mixing length
     178                   IF ( wall_adjustment )  THEN
     179                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
     180                                     l_stable )
     181                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
     182                   ELSE
     183                      l(k,j)  = MIN( l_grid(k), l_stable )
     184                      ll(k,j) = l_grid(k)
     185                   ENDIF
     186                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
     187                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
     188                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
     189                   ENDIF
     190
     191                ENDDO
     192             ENDDO
     193
     194!
     195!--          Calculate the tendency terms
     196             DO  j = nys, nyn
     197                DO  k = nzb_s_inner(j,i)+1, nzt
     198
     199                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
     200                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
     201
     202                    tend(k,j,i) = tend(k,j,i)                                  &
     203                                        + (                                    &
     204                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
     205                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
     206                                          ) * ddx2                             &
     207                                        + (                                    &
     208                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
     209                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
     210                                          ) * ddy2                             &
     211                                        + (                                    &
     212               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
     213             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
     214                                          ) * ddzw(k)                          &
     215                             - dissipation(k,j)
     216
     217                ENDDO
     218             ENDDO
     219
     220!
     221!--          Store dissipation if needed for calculating the sgs particle
     222!--          velocities
     223             IF ( use_sgs_for_particles )  THEN
     224                DO  j = nys, nyn
     225                   DO  k = nzb_s_inner(j,i)+1, nzt
     226                      diss(k,j,i) = dissipation(k,j)
     227                   ENDDO
     228                ENDDO
     229             ENDIF
     230
     231          ENDDO
     232
     233       ENDIF
    149234
    150235!
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |