Changeset 709


Ignore:
Timestamp:
Mar 30, 2011 9:31:40 AM (13 years ago)
Author:
raasch
Message:

formatting adjustments

Location:
palm/trunk/SOURCE
Files:
12 edited

Legend:

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

    r706 r709  
    11 MODULE advec_ws
    2 
    32
    43!-----------------------------------------------------------------------------!
    54! Current revisions:
    65! -----------------
     6! formatting adjustments
    77!
    88! Former revisions:
     
    4444!-----------------------------------------------------------------------------!
    4545
    46 
    4746    PRIVATE
    4847    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
     
    5251       MODULE PROCEDURE ws_init
    5352    END INTERFACE ws_init
     53
    5454    INTERFACE ws_statistics
    5555       MODULE PROCEDURE ws_statistics
    5656    END INTERFACE ws_statistics
     57
    5758    INTERFACE advec_s_ws
    5859       MODULE PROCEDURE advec_s_ws
    5960       MODULE PROCEDURE advec_s_ws_ij
    6061    END INTERFACE advec_s_ws
     62
    6163    INTERFACE advec_u_ws
    6264       MODULE PROCEDURE advec_u_ws
    6365       MODULE PROCEDURE advec_u_ws_ij
    6466    END INTERFACE advec_u_ws
     67
    6568    INTERFACE advec_v_ws
    6669       MODULE PROCEDURE advec_v_ws
    6770       MODULE PROCEDURE advec_v_ws_ij
    6871    END INTERFACE advec_v_ws
     72
    6973    INTERFACE advec_w_ws
    7074       MODULE PROCEDURE advec_w_ws
    7175       MODULE PROCEDURE advec_w_ws_ij
    7276    END INTERFACE advec_w_ws
     77
    7378    INTERFACE local_diss
    7479       MODULE PROCEDURE local_diss
     
    7883 CONTAINS
    7984
    80 !
    81 !-- Initialisation
    82 
     85
     86!------------------------------------------------------------------------------!
     87! Initialization of WS-scheme
     88!------------------------------------------------------------------------------!
    8389    SUBROUTINE ws_init
     90
    8491       USE arrays_3d
    8592       USE constants
     
    9198!--    Allocate arrays needed for dissipation control.
    9299       IF ( dissipation_control )  THEN
    93    !          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    94     !                  var_y(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    95      !                 var_z(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    96       !                gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    97        !               gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    98         !              gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     100!          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     101!                   var_y(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     102!                   var_z(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     103!                   gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     104!                   gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     105!                   gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
    99106       ENDIF
    100              
     107
     108!
    101109!--    Set the appropriate factors for scalar and momentum advection.
    102110       adv_sca_5 = 0.016666666666666
    103        adv_sca_3 = 0.083333333333333 
     111       adv_sca_3 = 0.083333333333333
    104112       adv_mom_5 = 0.0083333333333333
    105113       adv_mom_3 = 0.041666666666666
     
    108116!--    Arrays needed for statical evaluation of fluxes.
    109117       IF ( ws_scheme_mom )  THEN
    110           ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),         &
    111               sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),               &
    112               sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),                &
    113               sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),                &
    114               sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
    115 
     118
     119          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),  &
     120                    sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),  &
     121                    sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),   &
     122                    sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),   &
     123                    sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
     124
     125          sums_wsus_ws_l = 0.0
     126          sums_wsvs_ws_l = 0.0
     127          sums_us2_ws_l  = 0.0
     128          sums_vs2_ws_l  = 0.0
     129          sums_ws2_ws_l  = 0.0
     130
     131       ENDIF
     132
     133       IF ( ws_scheme_sca )  THEN
     134
     135          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions) )
     136          sums_wspts_ws_l = 0.0
     137
     138          IF ( humidity .OR. passive_scalar )  THEN
     139             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions) )
     140             sums_wsqs_ws_l = 0.0
     141          ENDIF
     142
     143          IF ( ocean )  THEN
     144             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions) )
     145             sums_wssas_ws_l = 0.0
     146          ENDIF
     147
     148       ENDIF
     149
     150!
     151!--    Arrays needed for reasons of speed optimization for cache and noopt
     152!--    version. For the vector version the buffer arrays are not necessary,
     153!--    because the the fluxes can swapped directly inside the loops of the
     154!--    advection routines.
     155       IF ( loop_optimization /= 'vector' )  THEN
     156
     157          IF ( ws_scheme_mom )  THEN
     158
     159             ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),  &
     160                       flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),  &
     161                       diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) )
     162             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn),               &
     163                       flux_l_v(nzb+1:nzt,nys:nyn),               &
     164                       flux_l_w(nzb+1:nzt,nys:nyn),               &
     165                       diss_l_u(nzb+1:nzt,nys:nyn),               &
     166                       diss_l_v(nzb+1:nzt,nys:nyn),               &
     167                       diss_l_w(nzb+1:nzt,nys:nyn) )
     168
     169          ENDIF
     170
     171          IF ( ws_scheme_sca )  THEN
     172
     173             ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), &
     174                       diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) )
     175             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn),              &
     176                       flux_l_e(nzb+1:nzt,nys:nyn),               &
     177                       diss_l_pt(nzb+1:nzt,nys:nyn),              &
     178                       diss_l_e(nzb+1:nzt,nys:nyn) )
     179
     180             IF ( humidity .OR. passive_scalar )  THEN
     181                ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) )
     182                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn),            &
     183                          diss_l_q(nzb+1:nzt,nys:nyn) )
     184             ENDIF
     185
     186             IF ( ocean )  THEN
     187                ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) )
     188                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn),           &
     189                          diss_l_sa(nzb+1:nzt,nys:nyn) )
     190             ENDIF
     191
     192          ENDIF
     193
     194       ENDIF
     195
     196!
     197!--    Determine the flags where the order of the scheme for horizontal
     198!--    advection has to be degraded.
     199       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
     200       DO  i = nxl, nxr
     201          DO  j = nys, nyn
     202
     203             boundary_flags(j,i) = 0
     204             IF ( outflow_l )  THEN
     205                IF ( i == nxlu  )  boundary_flags(j,i) = 5
     206                IF ( i == nxl   )  boundary_flags(j,i) = 6
     207             ELSEIF ( outflow_r )  THEN
     208                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
     209                IF ( i == nxr   )  boundary_flags(j,i) = 2
     210             ELSEIF ( outflow_n )  THEN
     211                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
     212                IF ( j == nyn   )  boundary_flags(j,i) = 4
     213             ELSEIF ( outflow_s )  THEN
     214                IF ( j == nysv  )  boundary_flags(j,i) = 7
     215                IF ( j == nys   )  boundary_flags(j,i) = 8
     216             ENDIF
     217
     218          ENDDO
     219       ENDDO
     220       
     221    END SUBROUTINE ws_init
     222
     223
     224!------------------------------------------------------------------------------!
     225! Initialize variables used for storing statistic qauntities (fluxes, variances)
     226!------------------------------------------------------------------------------!
     227    SUBROUTINE ws_statistics
     228   
     229       USE control_parameters
     230       USE statistics
     231
     232       IMPLICIT NONE
     233
     234!       
     235!--    The arrays needed for statistical evaluation are set to to 0 at the
     236!--    begin of prognostic_equations.
     237       IF ( ws_scheme_mom )  THEN
    116238          sums_wsus_ws_l = 0.0
    117239          sums_wsvs_ws_l = 0.0
     
    122244
    123245       IF ( ws_scheme_sca )  THEN
    124           ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions))
    125246          sums_wspts_ws_l = 0.0
    126           IF ( humidity  .OR.  passive_scalar )  THEN
    127              ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions))
    128              sums_wsqs_ws_l = 0.0
    129           ENDIF
    130           IF ( ocean )  THEN
    131              ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions))
    132              sums_wssas_ws_l = 0.0
    133           ENDIF
     247          IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
     248          IF ( ocean )  sums_wssas_ws_l = 0.0
     249
    134250       ENDIF
    135 !
    136 !--    Arrays needed for reasons of speed optimization for cache and noopt
    137 !--    version. For the vector version the buffer arrays are not necessary,
    138 !--    because the the fluxes can swapped directly inside the loops of the
    139 !--    advection routines.
    140        IF ( loop_optimization /= 'vector' )  THEN
    141           IF ( ws_scheme_mom )  THEN
    142              ALLOCATE(flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),            &
    143                 flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),                  &
    144                 diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt))
    145              ALLOCATE(flux_l_u(nzb+1:nzt,nys:nyn),                         &
    146                 flux_l_v(nzb+1:nzt,nys:nyn), flux_l_w(nzb+1:nzt,nys:nyn),  &
    147                 diss_l_u(nzb+1:nzt,nys:nyn), diss_l_v(nzb+1:nzt,nys:nyn),  &
    148                 diss_l_w(nzb+1:nzt,nys:nyn))
    149           ENDIF
    150           IF ( ws_scheme_sca )  THEN
    151              ALLOCATE(flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt),           &
    152                 diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt))
    153              ALLOCATE(flux_l_pt(nzb+1:nzt,nys:nyn),                        &
    154                 flux_l_e(nzb+1:nzt,nys:nyn), diss_l_pt(nzb+1:nzt,nys:nyn), &
    155                 diss_l_e(nzb+1:nzt,nys:nyn))
    156              IF ( humidity  .OR.  passive_scalar )  THEN
    157                 ALLOCATE(flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt))
    158                 ALLOCATE(flux_l_q(nzb+1:nzt,nys:nyn),                      &
    159                    diss_l_q(nzb+1:nzt,nys:nyn))
    160              ENDIF
    161              IF ( ocean )  THEN
    162                 ALLOCATE(flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt))
    163                 ALLOCATE(flux_l_sa(nzb+1:nzt,nys:nyn),                     &
    164                    diss_l_sa(nzb+1:nzt,nys:nyn))
    165              ENDIF
    166           ENDIF
    167        ENDIF
    168 !
    169 !--    Determine the flags where the order of the scheme for horizontal
    170 !--    advection should be degraded.
    171        ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
    172        DO  i = nxl, nxr
    173           DO  j = nys, nyn
    174              boundary_flags(j,i) = 0
    175              IF ( outflow_l )  THEN
    176                 IF ( i == nxlu )  boundary_flags(j,i) = 5
    177                 IF ( i == nxl )  boundary_flags(j,i) = 6
    178              ELSEIF ( outflow_r )  THEN
    179                 IF ( i == nxr-1 )  boundary_flags(j,i) = 1
    180                 IF ( i == nxr )  boundary_flags(j,i) = 2
    181              ELSEIF ( outflow_n )  THEN
    182                 IF ( j == nyn-1 )  boundary_flags(j,i) = 3
    183                 IF ( j == nyn )  boundary_flags(j,i) = 4
    184              ELSEIF ( outflow_s )  THEN
    185                 IF ( j == nysv )  boundary_flags(j,i) = 7
    186                 IF ( j == nys )  boundary_flags(j,i) = 8
    187              ENDIF
    188           ENDDO
    189        ENDDO
    190        
    191     END SUBROUTINE ws_init
    192    
    193 
    194    
    195     SUBROUTINE ws_statistics
    196    
    197        USE control_parameters
    198        USE statistics
    199 
    200 !       
    201 !--      The arrays needed for statistical evaluation are set to to 0 at the
    202 !--      begin of prognostic_equations.
    203           IF ( ws_scheme_mom )  THEN
    204              sums_wsus_ws_l = 0.0
    205              sums_wsvs_ws_l = 0.0
    206              sums_us2_ws_l = 0.0
    207              sums_vs2_ws_l = 0.0
    208              sums_ws2_ws_l = 0.0
    209           ENDIF
    210           IF ( ws_scheme_sca )  THEN
    211              sums_wspts_ws_l = 0.0
    212              IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0
    213              IF ( ocean )  sums_wssas_ws_l = 0.0
    214           ENDIF
    215251
    216252    END SUBROUTINE ws_statistics
    217 
    218253
    219254
     
    222257!------------------------------------------------------------------------------!
    223258    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
    224                swap_diss_y_local, swap_flux_x_local, swap_diss_x_local  )
     259                              swap_diss_y_local, swap_flux_x_local, &
     260                              swap_diss_x_local  )
    225261
    226262       USE arrays_3d
     
    239275       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
    240276                                             flux_n, diss_n
    241        REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &       
     277       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &
    242278                                             swap_diss_y_local
    243279       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
    244280                                             swap_diss_x_local
    245281       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
    246        
     282
     283
    247284       degraded_l = .FALSE.
    248285       degraded_s = .FALSE.
     
    250287       IF ( boundary_flags(j,i) /= 0 )  THEN
    251288!       
    252 !--      Degrade the order for Dirichlet bc. at the outflow boundary.
    253          SELECT CASE ( boundary_flags(j,i) )
    254          
    255             CASE ( 1 )
    256                DO  k = nzb_s_inner(j,i) + 1, nzt
    257                   u_comp    = u(k,j,i+1) - u_gtrans
    258                   flux_r(k) = u_comp * (                                      &
    259                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    260                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    261                   diss_r(k) = - abs(u_comp) * (                               &
    262                               3. * ( sk(k,j,i+1)-sk(k,j,i)     )              &
    263                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     289!--       Degrade the order for Dirichlet bc. at the outflow boundary
     290          SELECT CASE ( boundary_flags(j,i) )
     291
     292             CASE ( 1 )
     293
     294                DO  k = nzb_s_inner(j,i)+1, nzt
     295                   u_comp    = u(k,j,i+1) - u_gtrans
     296                   flux_r(k) = u_comp * (                                      &
     297                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     298                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     299                   diss_r(k) = -ABS( u_comp ) * (                              &
     300                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     301                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     302                ENDDO
     303
     304             CASE ( 2 )
     305
     306                DO  k = nzb_s_inner(j,i)+1, nzt
     307                   u_comp    = u(k,j,i+1) - u_gtrans
     308                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
     309                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
     310                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
     311                                         0.5, ddx )
     312                ENDDO
     313
     314             CASE ( 3 )
     315
     316                DO  k = nzb_s_inner(j,i) + 1, nzt
     317                   v_comp = v(k,j+1,i) - v_gtrans
     318                   flux_n(k) = v_comp * (                                      &
     319                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     320                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     321                   diss_n(k) = -ABS( v_comp ) * (                              &
     322                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)     )           &
     323                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    264324               ENDDO
    265                
    266             CASE ( 2 )
    267                DO  k = nzb_s_inner(j,i) + 1, nzt
    268                   u_comp    = u(k,j,i+1) - u_gtrans
    269                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
    270                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
    271                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
    272                                         0.5, ddx )
    273               ENDDO
    274              
    275             CASE ( 3 )
    276                DO  k = nzb_s_inner(j,i) + 1, nzt
    277                   v_comp = v(k,j+1,i) - v_gtrans
    278                   flux_n(k) = v_comp * (                                      &
    279                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    280                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    281                   diss_n(k) = - abs(v_comp) * (                               &
    282                               3. * ( sk(k,j+1,i)-sk(k,j,i)    )               &
    283                               -    (sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    284               ENDDO
    285             CASE ( 4 )
    286                DO  k = nzb_s_inner(j,i) + 1, nzt
    287                   v_comp    = v(k,j+1,i) - v_gtrans
    288                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
    289                   diss_n(k) =  diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
     325
     326             CASE ( 4 )
     327
     328                DO  k = nzb_s_inner(j,i)+1, nzt
     329                   v_comp    = v(k,j+1,i) - v_gtrans
     330                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
     331                   diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
    290332                                         sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
    291333                                         0.5, ddy )     
    292                ENDDO
     334                ENDDO
     335
     336             CASE ( 5 )
     337
     338                DO  k = nzb_w_inner(j,i)+1, nzt
     339                   u_comp    = u(k,j,i+1) - u_gtrans
     340                   flux_r(k) = u_comp  * (                                     &
     341                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     342                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     343                   diss_r(k) = -ABS( u_comp ) * (                              &
     344                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     345                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     346                ENDDO 
     347
     348             CASE ( 6 )
     349
     350                DO  k = nzb_s_inner(j,i)+1, nzt
     351                   u_comp    = u(k,j,i+1) - u_gtrans
     352                   flux_r(k) = u_comp * (                                      &
     353                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     354                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     355                   diss_r(k) = -ABS( u_comp ) * (                              &
     356                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     357                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     358!                             
     359!--                Compute leftside fluxes for the left boundary of PE domain
     360                   u_comp                 = u(k,j,i) - u_gtrans
     361                   swap_flux_x_local(k,j) = u_comp * (                         &
     362                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
     363                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
     364                                                      sk(k,j,i), sk(k,j,i-1),  &
     365                                                      sk(k,j,i-1), u_comp,     &
     366                                                      0.5, ddx )
     367                ENDDO
     368                degraded_l = .TRUE.
     369
     370             CASE ( 7 )
     371
     372                DO  k = nzb_s_inner(j,i)+1, nzt
     373                   v_comp    = v(k,j+1,i)-v_gtrans
     374                   flux_n(k) = v_comp * (                                      &
     375                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     376                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     377                   diss_n(k) = -ABS( v_comp ) * (                              &
     378                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
     379                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     380                ENDDO
     381
     382             CASE ( 8 )
     383
     384                DO  k = nzb_s_inner(j,i)+1, nzt
     385                   v_comp    = v(k,j+1,i) - v_gtrans
     386                   flux_n(k) = v_comp * (                                      &
     387                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     388                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     389                   diss_n(k) = -ABS( v_comp ) * (                              &
     390                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
     391                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     392!                             
     393!--                Compute southside fluxes for the south boundary of PE domain
     394                   v_comp               = v(k,j,i) - v_gtrans
     395                   swap_flux_y_local(k) = v_comp *                             &
     396                                          ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
     397                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
     398                                                    sk(k,j,i), sk(k,j-1,i),    &
     399                                                    sk(k,j-1,i), v_comp,       &
     400                                                    0.5, ddy )
     401                ENDDO
     402                degraded_s = .TRUE.
    293403               
    294             CASE ( 5 )
    295                DO  k = nzb_w_inner(j,i) + 1, nzt
    296                   u_comp    = u(k,j,i+1) - u_gtrans
    297                   flux_r(k) = u_comp  * (                                     &
    298                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    299                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    300                   diss_r(k) = - abs(u_comp) * (                               &
    301                               3. * ( sk(k,j,i+1)-sk(k,j,i) )                  &
    302                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    303                ENDDO 
    304                
    305             CASE ( 6 )
    306                DO  k = nzb_s_inner(j,i) + 1, nzt
    307                   u_comp    = u(k,j,i+1) - u_gtrans
    308                   flux_r(k) = u_comp * (                                      &
    309                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    310                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    311                   diss_r(k) = - abs(u_comp) * (                               &
    312                               3. * ( sk(k,j,i+1)-sk(k,j,i)     )              &
    313                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    314 !                             
    315 !--               Compute leftside fluxes for the left boundary of PE domain.
    316                   u_comp                 = u(k,j,i) - u_gtrans
    317                   swap_flux_x_local(k,j) = u_comp * (                         &
    318                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
    319                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
    320                                                      sk(k,j,i), sk(k,j,i-1),  &
    321                                                      sk(k,j,i-1), u_comp,     &
    322                                                      0.5, ddx )
    323                ENDDO
    324                degraded_l = .TRUE.
    325                
    326             CASE ( 7 )
    327                DO  k = nzb_s_inner(j,i) + 1, nzt
    328                   v_comp    = v(k,j+1,i)-v_gtrans
    329                   flux_n(k) = v_comp * (                                      &
    330                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    331                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    332                   diss_n(k) = - abs(v_comp) * (                               &
    333                               3. * ( sk(k,j+1,i) - sk(k,j,i)   )              &
    334                               -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    335                ENDDO
    336                
    337             CASE ( 8 )
    338                DO  k = nzb_s_inner(j,i) + 1, nzt
    339                   v_comp    = v(k,j+1,i) - v_gtrans
    340                   flux_n(k) = v_comp * (                                      &
    341                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    342                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    343                   diss_n(k) = - abs(v_comp) * (                               &
    344                               3. * ( sk(k,j+1,i) - sk(k,j,i) )                &
    345                               -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    346 !                             
    347 !--               Compute southside fluxes for the south boundary of PE domain           
    348                   v_comp               = v(k,j,i) - v_gtrans
    349                   swap_flux_y_local(k) = v_comp * (                           &
    350                                          sk(k,j,i)+ sk(k,j-1,i) ) * 0.5
    351                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
    352                                                    sk(k,j,i), sk(k,j-1,i),    &
    353                                                    sk(k,j-1,i), v_comp,       &
    354                                                    0.5, ddy )
    355                ENDDO
    356                degraded_s = .TRUE.
    357                
    358             CASE DEFAULT
     404             CASE DEFAULT
    359405           
    360          END SELECT
     406          END SELECT
     407
    361408!         
    362 !--      Compute the crosswise 5th order fluxes at the outflow
    363          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
    364          .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    365             DO  k = nzb_s_inner(j,i) + 1, nzt
    366                v_comp    = v(k,j+1,i) - v_gtrans
    367                flux_n(k) = v_comp * (                                         &
    368                            37.  * ( sk(k,j+1,i) + sk(k,j,i)   )               &
    369                            - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
    370                            +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    371                diss_n(k) = - abs(v_comp) * (                                  &
    372                            10. *  ( sk(k,j+1,i) - sk(k,j,i)   )               &
    373                            - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
    374                            +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    375             ENDDO
     409!--       Compute the crosswise 5th order fluxes at the outflow
     410          IF ( boundary_flags(j,i) == 1  .OR.  boundary_flags(j,i) == 2  .OR. &
     411               boundary_flags(j,i) == 5  .OR.  boundary_flags(j,i) == 6 )  THEN
     412
     413             DO  k = nzb_s_inner(j,i)+1, nzt
     414                v_comp    = v(k,j+1,i) - v_gtrans
     415                flux_n(k) = v_comp * (                                         &
     416                            37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )               &
     417                          -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
     418                          +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     419                diss_n(k) = -ABS( v_comp ) * (                                 &
     420                            10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )               &
     421                          -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
     422                          +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     423             ENDDO
    376424           
    377          ELSE
     425          ELSE
    378426         
    379             DO  k = nzb_s_inner(j,i) + 1, nzt
    380                u_comp    = u(k,j,i+1) - u_gtrans 
    381                flux_r(k) = u_comp * (                                         &
    382                            37. * ( sk(k,j,i+1) + sk(k,j,i)   )               &
    383                            - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
    384                            +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    385                diss_r(k) = - abs(u_comp) * (                                  &
    386                            10. * ( sk(k,j,i+1) - sk(k,j,i)   )               &
    387                            - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
    388                            +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    389             ENDDO
     427             DO  k = nzb_s_inner(j,i)+1, nzt
     428                u_comp    = u(k,j,i+1) - u_gtrans 
     429                flux_r(k) = u_comp * (                                         &
     430                            37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )               &
     431                          -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
     432                                ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     433                diss_r(k) = -ABS( u_comp ) * (                                 &
     434                            10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )               &
     435                          -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
     436                                ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     437             ENDDO
    390438           
    391          ENDIF
    392 !
    393 !--    Compute the fifth order fluxes for the interior of PE domain.
     439          ENDIF
     440
    394441       ELSE
    395442               
    396           DO  k = nzb_u_inner(j,i) + 1, nzt
     443!
     444!--       Compute the fifth order fluxes for the interior of PE domain.
     445          DO  k = nzb_u_inner(j,i)+1, nzt
    397446             u_comp    = u(k,j,i+1) - u_gtrans
    398447             flux_r(k) = u_comp * (                                           &
    399                          37.  * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
    400                          - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
    401                          +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    402              diss_r(k) = - abs(u_comp) * (                                    &
    403                          10.  * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
    404                          - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
    405                          +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     448                         37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     449                       -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
     450                             ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     451             diss_r(k) = -ABS( u_comp ) * (                                   &
     452                         10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
     453                       -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
     454                             ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    406455
    407456             v_comp    = v(k,j+1,i) - v_gtrans
    408457             flux_n(k) = v_comp * (                                           &
    409                          37.  * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
    410                          - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
    411                          +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    412              diss_n(k) = - abs(v_comp) * (                                    &
    413                          10.  * ( sk(k,j+1,i) - sk(k,j,i)   )                 &                     
    414                          - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    415                          +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5               
     458                         37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
     459                       -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
     460                             ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     461             diss_n(k) = -ABS( v_comp ) * (                                   &
     462                         10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
     463                       -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
     464                       +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    416465          ENDDO
    417466         
     
    424473              v_comp               = v(k,j,i) - v_gtrans
    425474              swap_flux_y_local(k) = v_comp * (                               &
    426                                      37. *    ( sk(k,j,i) + sk(k,j-1,i)   )   &
    427                                      - 8. *   ( sk(k,j+1,i) + sk(k,j-2,i) )   &
    428                                      +        ( sk(k,j+2,i) + sk(k,j-3,i) )   &
    429                                      ) * adv_sca_5
    430               swap_diss_y_local(k) = - abs(v_comp) * (                        &
    431                                      10. *  ( sk(k,j,i) - sk(k,j-1,i)  )     &
    432                                      - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
    433                                      +        sk(k,j+2,i) - sk(k,j-3,i)       &
    434                                      ) * adv_sca_5
     475                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
     476                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
     477                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
     478                                              ) * adv_sca_5
     479              swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
     480                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
     481                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     482                                           sk(k,j+2,i) - sk(k,j-3,i)       &
     483                                                      ) * adv_sca_5
    435484           ENDDO           
    436485         
    437486        ENDIF
    438487       
    439         IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
     488        IF ( i == nxl  .AND.  .NOT. degraded_l )  THEN
    440489       
    441            DO  k = nzb_s_inner(j,i) + 1, nzt
     490           DO  k = nzb_s_inner(j,i)+1, nzt
    442491              u_comp                 = u(k,j,i) - u_gtrans
    443492              swap_flux_x_local(k,j) = u_comp * (                             &
    444                                        37.  * ( sk(k,j,i) + sk(k,j,i-1)  )   &
    445                                        - 8. * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
    446                                        +      ( sk(k,j,i+2) + sk(k,j,i-3) )   &
    447                                        ) * adv_sca_5
    448               swap_diss_x_local(k,j) = - abs(u_comp) * (                      &
    449                                        10.  * ( sk(k,j,i) - sk(k,j,i-1)  )   &
    450                                        - 5. * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
    451                                        +      ( sk(k,j,i+2) - sk(k,j,i-3) )   &
    452                                        ) * adv_sca_5
     493                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
     494                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
     495                                           ( sk(k,j,i+2) + sk(k,j,i-3) )   &
     496                                                ) * adv_sca_5
     497              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     498                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
     499                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
     500                                           ( sk(k,j,i+2) - sk(k,j,i-3) )   &
     501                                                        ) * adv_sca_5
    453502           ENDDO
    454503           
    455504        ENDIF
     505
    456506!       
    457 !--    Now compute the tendency terms for the horizontal parts.
    458        DO  k = nzb_s_inner(j,i) + 1, nzt
     507!--    Now compute the tendency terms for the horizontal parts
     508       DO  k = nzb_s_inner(j,i)+1, nzt
    459509         
    460           tend(k,j,i) = tend(k,j,i) - (                                       &
    461                ( flux_r(k) + diss_r(k)                                        &
    462              - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx    &
    463              + (flux_n(k)+diss_n(k)                                           &
    464              - (swap_flux_y_local(k)    + swap_diss_y_local(k) ) )   * ddy )
     510          tend(k,j,i) = tend(k,j,i) - (                                      &
     511                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
     512                            swap_diss_x_local(k,j) ) * ddx                   &
     513                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   - &
     514                            swap_diss_y_local(k)   ) * ddy                   &
     515                                      )
    465516
    466517          swap_flux_y_local(k)   = flux_n(k)
     
    470521         
    471522       ENDDO
    472 !
    473 !--    Vertical advection, degradation of order near surface and top.
    474 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    475 !--    statistical evaluation the top flux at the surface should be 0
    476 
    477        flux_t(nzb_s_inner(j,i)) = 0.
    478        diss_t(nzb_s_inner(j,i)) = 0.
     523
     524!
     525!--    Vertical advection, degradation of order near bottom and top.
     526!--    The fluxes flux_d and diss_d at the surface are 0. Due to later
     527!--    calculation of statistics the top flux at the surface should be 0.
     528       flux_t(nzb_s_inner(j,i)) = 0.0
     529       diss_t(nzb_s_inner(j,i)) = 0.0
     530
    479531!       
    480 !--    2nd order scheme       
     532!--    2nd-order scheme (bottom)
    481533       k         = nzb_s_inner(j,i) + 1
    482534       flux_d    = flux_t(k-1)
    483535       diss_d    = diss_t(k-1)
    484536       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
     537
    485538!       
    486 !--    sk(k,j,i) is referenced three times to avoid a access below surface.
     539!--    sk(k,j,i) is referenced three times to avoid an access below surface
    487540       diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
    488541                             sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
    489542                   
    490        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                   &
    491                                  - ( flux_d + diss_d )     ) * ddzw(k)
    492 !
    493 !--    WS3 as an intermediate step      
     543       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     544                                   * ddzw(k)
     545!
     546!--    WS3 as an intermediate step (bottom)
    494547       k         = nzb_s_inner(j,i) + 2
    495548       flux_d    = flux_t(k-1)
    496549       diss_d    = diss_t(k-1)
    497550       flux_t(k) = w(k,j,i) * (                                               &
    498                    7. * ( sk(k+1,j,i) + sk(k,j,i)   )                         &
    499                    -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    500        diss_t(k) = - abs(w(k,j,i)) * (                                        &
    501                    3. * ( sk(k+1,j,i) - sk(k,j,i)   )                         &
    502                    -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    503 
    504        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    505                                  - ( flux_d + diss_d ) ) * ddzw(k)
     551                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )           &
     552                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )           &
     553                              ) * adv_sca_3
     554       diss_t(k) = -ABS( w(k,j,i) ) * (                                       &
     555                                3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )           &
     556                                    - ( sk(k+2,j,i) - sk(k-1,j,i) )           &
     557                                      ) * adv_sca_3
     558
     559       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     560                                   * ddzw(k)
    506561!
    507562!--    WS5
    508        DO  k = nzb_s_inner(j,i) + 3, nzt - 2
     563       DO  k = nzb_s_inner(j,i)+3, nzt-2
    509564       
    510565          flux_d    = flux_t(k-1)
    511           diss_d    = diss_t(k-1)
    512          
    513           flux_t(k) = w(k,j,i) * (                                            & 
    514                       37.  * ( sk(k+1,j,i) + sk(k,j,i)   )                    &
    515                       - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )                    &
    516                       +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) *adv_sca_5
    517           diss_t(k) = - abs(w(k,j,i)) * (                                     &
    518                       10.  * ( sk(k+1,j,i) - sk(k,j,i)   )                    &     
    519                       - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )                    &
    520                       +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
     566          diss_d    = diss_t(k-1)         
     567          flux_t(k) = w(k,j,i) * (                                            &
     568                                   37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )       &
     569                                 -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )       &
     570                                 +        ( sk(k+3,j,i) + sk(k-2,j,i) )       &
     571                                 ) * adv_sca_5
     572          diss_t(k) = -ABS( w(k,j,i) ) * (                                     &
     573                                           10.0 * ( sk(k+1,j,i) - sk(k,j,i)   )&
     574                                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
     575                                         +        ( sk(k+3,j,i) - sk(k-2,j,i) )&
     576                                         ) * adv_sca_5
    521577
    522578          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    523579                                    - ( flux_d + diss_d ) ) * ddzw(k)
     580
    524581       ENDDO
    525 !
    526 !--    WS3 as an intermediate step
     582
     583!
     584!--    WS3 as an intermediate step (top)
    527585       k         = nzt - 1
    528586       flux_d    = flux_t(k-1)
    529587       diss_d    = diss_t(k-1)
    530        flux_t(k) = w(k,j,i) * (                                               &
    531                    7. * ( sk(k+1,j,i) + sk(k,j,i) )                           &
    532                    -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    533        diss_t(k) = - abs(w(k,j,i)) * (                                        &
    534                    3. * ( sk(k+1,j,i) - sk(k,j,i) )                           &
    535                    -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    536 
    537        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    538                                  - ( flux_d + diss_d ) ) * ddzw(k)
     588       flux_t(k) = w(k,j,i) * (                                                &
     589                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )            &
     590                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )            &
     591                              ) * adv_sca_3
     592       diss_t(k) = -ABS( w(k,j,i) ) * (                                        &
     593                                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )    &
     594                                            - ( sk(k+2,j,i) - sk(k-1,j,i) )    &
     595                                      ) * adv_sca_3
     596
     597       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     598                                   * ddzw(k)
    539599!                                 
    540 !--    2nd
     600!--    2nd-order scheme (top)
    541601       k         = nzt
    542602       flux_d    = flux_t(k-1)
    543603       diss_d    = diss_t(k-1)
    544        flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.5
     604       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
     605
    545606!       
    546 !--    sk(k+1) is referenced two times to avoid a segmentation fault at top.
    547        diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i),&
     607!--    sk(k+1) is referenced two times to avoid a segmentation fault at top
     608       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
    548609                             sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
    549610
    550        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                   &
    551                                  - ( flux_d + diss_d ) ) * ddzw(k)
     611       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     612                                   * ddzw(k)
    552613!       
    553 !--    evaluation of statistics
     614!--    Evaluation of statistics
    554615       SELECT CASE ( sk_char )
    555616
    556617          CASE ( 'pt' )
    557             DO  k = nzb_s_inner(j,i), nzt
    558                sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                     &
    559                  + (flux_t(k)+diss_t(k))                                       &
    560                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     618
     619             DO  k = nzb_s_inner(j,i), nzt
     620                sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) +                  &
     621                                       ( flux_t(k) + diss_t(k) )               &
     622                                 * weight_substep(intermediate_timestep_count) &
     623                                       * rmask(j,i,:)
    561624             ENDDO
    562625             
    563626          CASE ( 'sa' )
    564             DO  k = nzb_s_inner(j,i), nzt
    565                sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                     &
    566                  +(flux_t(k)+diss_t(k))                                        &
    567                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     627
     628             DO  k = nzb_s_inner(j,i), nzt
     629                sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) +                  &
     630                                       ( flux_t(k) + diss_t(k) )               &
     631                                 * weight_substep(intermediate_timestep_count) &
     632                                       * rmask(j,i,:)
    568633             ENDDO
    569634             
    570635          CASE ( 'q' )
     636
    571637             DO  k = nzb_s_inner(j,i), nzt
    572                 sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                      &
    573                   + (flux_t(k)+diss_t(k))                                      &
    574                   * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     638                sums_wsqs_ws_l(k,:)  = sums_wsqs_ws_l(k,:) +                   &
     639                                      ( flux_t(k) + diss_t(k) )                &
     640                                 * weight_substep(intermediate_timestep_count) &
     641                                      * rmask(j,i,:)
    575642             ENDDO
    576643
  • palm/trunk/SOURCE/exchange_horiz.f90

    r708 r709  
    44! Current revisions:
    55! -----------------
    6 !
     6! formatting adjustments
    77!
    88! Former revisions:
     
    6969
    7070!
    71 !-- Exchange of lateral boundary values for parallel computers
     71!-- Exchange of lateral boundary values
    7272    IF ( pdims(1) == 1  .OR.  mg_switch_to_pe0 )  THEN
    7373!
  • palm/trunk/SOURCE/flow_statistics.f90

    r700 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
     
    114115    USE arrays_3d
    115116    USE cloud_parameters
     117    USE control_parameters
    116118    USE cpulog
    117119    USE grid_variables
     
    120122    USE pegrid
    121123    USE statistics
    122     USE control_parameters
    123124
    124125    IMPLICIT NONE
     
    738739          ENDDO
    739740       ENDDO
    740 !-     for reasons of speed optimization the loop is splitted, to avoid if-else
    741 !-     statements inside the loop
    742 !-     Fluxes which have been computed in part directly inside the advection routines
    743 !-     treated seperatly.
    744 !-     First treat the momentum fluxes
     741!
     742!--    For speed optimization fluxes which have been computed in part directly
     743!--    inside the WS advection routines are treated seperatly
     744!--    Momentum fluxes first:
    745745       IF ( .NOT. ws_scheme_mom )  THEN
    746746         !$OMP DO
     
    771771         DO  i = nxl, nxr
    772772            DO  j = nys, nyn
    773                DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
    774 !-                vertical heat flux
     773               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
     774!
     775!--               Vertical heat flux
    775776                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
    776777                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
     
    12211222
    12221223 END SUBROUTINE flow_statistics
    1223 
    1224 
    1225 
  • palm/trunk/SOURCE/inflow_turbulence.f90

    r668 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
     
    5455!
    5556!-- First, local averaging within the recycling domain
    56 
    5757    i = recycling_plane
    5858
     
    7979!-- Now, averaging over all PEs
    8080    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    81     CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  &
    82                     MPI_REAL, MPI_SUM, comm2d, ierr )
     81    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr, MPI_REAL, &
     82                        MPI_SUM, comm2d, ierr )
    8383
    8484#else
     
    165165          DO  k = nzb, nzt + 1
    166166
    167               u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
     167              u(k,j,-nbgp+1:0) = mean_inflow_profiles(k,1) + &
    168168                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
    169169              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
    170170                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
    171               w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
     171              w(k,j,-nbgp:-1)  =                             &
     172                           inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
    172173              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
    173174                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
     
    181182    ENDIF
    182183
    183 !
    184 !-- Conserve the volume flow at the inflow in order to avoid generation of
    185 !-- waves in the stable layer
    186 !    IF ( conserve_volume_flow  .AND.  inflow_l )  THEN
    187 
    188 !       volume_flow(1)   = 0.0
    189 !       volume_flow_l(1) = 0.0
    190 
    191 !       i = 0
    192 
    193 !       DO  j = nys, nyn
    194 !
    195 !--       Sum up the volume flow through the south/north boundary
    196 !          DO  k = nzb_2d(j,i) + 1, nzt
    197 !             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
    198 !          ENDDO
    199 !       ENDDO
    200 
    201 #if defined( __parallel )   
    202 !       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    203 !       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
    204 !                           MPI_SUM, comm1dy, ierr )   
    205 #else
    206 !       volume_flow = volume_flow_l 
    207 #endif
    208 !       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
    209 !                               / volume_flow_area(1)
    210 
    211 !       DO  j = nys-1, nyn+1
    212 !          DO  k = nzb_v_inner(j,i) + 1, nzt
    213 !             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
    214 !          ENDDO
    215 !       ENDDO
    216 
    217 !    ENDIF
    218 
    219184    CALL cpu_log( log_point(40), 'inflow_turbulence', 'stop' )
    220185
  • palm/trunk/SOURCE/init_3d_model.f90

    r708 r709  
    77! Current revisions:
    88! -----------------
    9 !
     9! formatting adjustments
    1010!
    1111! Former revisions:
     
    2929! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
    3030! allocation of arrays. Calls of exchange_horiz are modified.
    31 ! Call ws_init to initialize arrays needed for statistical evaluation and
     31! Call ws_init to initialize arrays needed for calculating statisticas and for
    3232! optimization when ws-scheme is used.
    3333! Initial volume flow is now calculated by using the variable hom_sum.
     
    3535! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
    3636! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
    37 ! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
    38 ! representative for the height z0 
     37! mirror to Dirichlet boundary conditions (u=v=0), so that k=nzb is
     38! representative for the height z0.
    3939! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
    4040!
     
    453453   
    454454!
    455 !-- Allocate arrays containing the RK coefficient for right evaluation of
    456 !-- perturbation pressure and turbulent fluxes. At this point it is needed
    457 !-- for right pressure correction during initialization. Further below
    458 !-- the real values will be set.
    459     ALLOCATE (weight_substep(1:intermediate_timestep_count_max),           &
     455!-- Allocate arrays containing the RK coefficient for calculation of
     456!-- perturbation pressure and turbulent fluxes. At this point values are
     457!-- set for pressure calculation during initialization (where no timestep
     458!-- is done). Further below the values needed within the timestep scheme
     459!-- will be set.
     460    ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
    460461              weight_pres(1:intermediate_timestep_count_max) )
    461     weight_substep = 1.
    462     weight_pres = 1.
    463     intermediate_timestep_count = 1 ! needed for simulated_time=0
     462    weight_substep = 1.0
     463    weight_pres    = 1.0
     464    intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
    464465       
    465466!
     
    10111012!
    10121013!--    Read binary data from restart file
    1013 
    10141014       CALL read_3d_binary
    10151015
     
    11191119!
    11201120!-- Calculate the initial volume flow at the right and north boundary
    1121     IF ( conserve_volume_flow ) THEN
     1121    IF ( conserve_volume_flow )  THEN
    11221122
    11231123       volume_flow_initial_l = 0.0
     
    11281128          IF ( nxr == nx )  THEN
    11291129             DO  j = nys, nyn
    1130                 DO  k = nzb_2d(j,nx) + 1, nzt
     1130                DO  k = nzb_2d(j,nx)+1, nzt
    11311131                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    11321132                                              hom_sum(k,1,0) * dzw(k)
     
    11381138          IF ( nyn == ny )  THEN
    11391139             DO  i = nxl, nxr
    1140                 DO  k = nzb_2d(ny,i) + 1, nzt 
     1140                DO  k = nzb_2d(ny,i)+1, nzt 
    11411141                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    1142                                                hom_sum(k,2,0) * dzw(k)
     1142                                              hom_sum(k,2,0) * dzw(k)
    11431143                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
    11441144                ENDDO
     
    11501150          IF ( nxr == nx )  THEN
    11511151             DO  j = nys, nyn
    1152                 DO  k = nzb_2d(j,nx) + 1, nzt
     1152                DO  k = nzb_2d(j,nx)+1, nzt
    11531153                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    1154                                                u(k,j,nx) * dzw(k)
     1154                                              u(k,j,nx) * dzw(k)
    11551155                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
    11561156                ENDDO
     
    11601160          IF ( nyn == ny )  THEN
    11611161             DO  i = nxl, nxr
    1162                 DO  k = nzb_2d(ny,i) + 1, nzt 
     1162                DO  k = nzb_2d(ny,i)+1, nzt 
    11631163                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    11641164                                              v(k,ny,i) * dzw(k)
     
    11711171
    11721172#if defined( __parallel )
    1173        CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
     1173       CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), &
    11741174                           2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1175        CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
     1175       CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),       &
    11761176                           2, MPI_REAL, MPI_SUM, comm2d, ierr )
    11771177
     
    11821182
    11831183!
    1184 !--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
    1185 !--       and calculated from u|v_bulk instead.
     1184!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
     1185!--    from u|v_bulk instead
    11861186       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
    11871187          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
     
    12081208!
    12091209!-- Initialization of the leaf area density
    1210      IF ( plant_canopy ) THEN
     1210    IF ( plant_canopy ) THEN
    12111211 
    12121212       SELECT CASE ( TRIM( canopy_mode ) )
     
    12181218                   lad_s(:,j,i) = lad(:)
    12191219                   cdc(:,j,i)   = drag_coefficient
    1220                    IF ( passive_scalar ) THEN
     1220                   IF ( passive_scalar )  THEN
    12211221                      sls(:,j,i) = leaf_surface_concentration
    12221222                      sec(:,j,i) = scalar_exchange_coefficient
     
    12401240       CALL exchange_horiz( cdc, nbgp )
    12411241
    1242        IF ( passive_scalar ) THEN
     1242       IF ( passive_scalar )  THEN
    12431243          CALL exchange_horiz( sls, nbgp )
    12441244          CALL exchange_horiz( sec, nbgp )
     
    12551255          DO  j = nys, nyn
    12561256             DO  k = nzb, nzt+1
    1257                 IF ( lad_s(k,j,i) > 0.0 ) THEN
     1257                IF ( lad_s(k,j,i) > 0.0 )  THEN
    12581258                   lad_u(k,j,i)   = lad_s(k,j,i)
    12591259                   lad_u(k,j,i+1) = lad_s(k,j,i)
     
    12771277!
    12781278!--    Initialisation of the canopy heat source distribution
    1279        IF ( cthf /= 0.0 ) THEN
     1279       IF ( cthf /= 0.0 )  THEN
    12801280!
    12811281!--       Piecewise evaluation of the leaf area index by
     
    13131313          shf(:,:) = canopy_heat_flux(0,:,:)
    13141314
    1315           IF ( ASSOCIATED( shf_m ) ) shf_m = shf
     1315          IF ( ASSOCIATED( shf_m ) )  shf_m = shf
    13161316
    13171317       ENDIF
     
    13461346
    13471347!
    1348 !-- Setting weighting factors for right evaluation of perturbation pressure
    1349 !-- and turbulent quantities during the RK substeps.               
    1350     IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN    ! RK3       
     1348!-- Setting weighting factors for calculation of perturbation pressure
     1349!-- and turbulent quantities from the RK substeps               
     1350    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
     1351
    13511352       weight_substep(1) = 0.166666666666666
    13521353       weight_substep(2) = 0.3
    13531354       weight_substep(3) = 0.533333333333333
    1354          
    1355        weight_pres(1) = 0.333333333333333
    1356        weight_pres(2) = 0.416666666666666
    1357        weight_pres(3) = 0.25         
    1358     ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! RK2       
     1355
     1356       weight_pres(1)    = 0.333333333333333
     1357       weight_pres(2)    = 0.416666666666666
     1358       weight_pres(3)    = 0.25
     1359
     1360    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
     1361
    13591362       weight_substep(1) = 0.5
    13601363       weight_substep(2) = 0.5
    13611364         
    1362        weight_pres(1) = 0.5
    1363        weight_pres(2) = 0.5         
    1364     ELSE                                                  ! Euler and Leapfrog       
     1365       weight_pres(1)    = 0.5
     1366       weight_pres(2)    = 0.5         
     1367
     1368    ELSE                                     ! for Euler- and leapfrog-method
     1369
    13651370       weight_substep(1) = 1.0     
    1366        weight_pres(1) = 1.0                   
     1371       weight_pres(1)    = 1.0                   
     1372
    13671373    ENDIF
    13681374
     
    14761482
    14771483!
    1478 !-- Initialize local summation arrays for UP flow_statistics. This is necessary
    1479 !-- because they may not yet have been initialized when they are called from
    1480 !-- flow_statistics (or - depending on the chosen model run - are never
    1481 !-- initialized)
     1484!-- Initialize local summation arrays for routine flow_statistics.
     1485!-- This is necessary because they may not yet have been initialized when they
     1486!-- are called from flow_statistics (or - depending on the chosen model run -
     1487!-- are never initialized)
    14821488    sums_divnew_l      = 0.0
    14831489    sums_divold_l      = 0.0
     
    14921498!
    14931499!-- User-defined initializing actions. Check afterwards, if maximum number
    1494 !-- of allowed timeseries is not exceeded
     1500!-- of allowed timeseries is exceeded
    14951501    CALL user_init
    14961502
  • palm/trunk/SOURCE/init_coupling.f90

    r692 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
  • palm/trunk/SOURCE/init_grid.f90

    r708 r709  
    44! Current revisions:
    55! -----------------
    6 !
     6! formatting adjustments
    77!
    88! Former revisions:
     
    8989   
    9090!
    91 !   Computation of the array bounds.
     91!-- Calculation of horizontal array bounds including ghost layers
    9292    nxlg = nxl - nbgp
    9393    nxrg = nxr + nbgp
    9494    nysg = nys - nbgp
    9595    nyng = nyn + nbgp
     96
    9697!
    9798!-- Allocate grid arrays
     
    196197   
    197198!   
    198 !-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and
    199 !-- modify the lowest entry. This is necessary for atmosphere runs, because
    200 !-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays
    201 !-- pressure solver uses wrong values and this causes kinks in the profiles
    202 !-- of turbulent quantities. 
    203     IF ( psolver /= 'multigrid' ) THEN
     199!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
     200!-- everywhere. For the actual grid, the grid spacing at the lowest level
     201!-- is only dz/2, but should be dz. Therefore, an additional array
     202!-- containing with appropriate grid information is created for these
     203!-- solvers.
     204    IF ( psolver /= 'multigrid' )  THEN
    204205       ALLOCATE( ddzu_pres(1:nzt+1) )
    205206       ddzu_pres = ddzu
    206        IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)
     207       IF( .NOT. ocean )  ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
    207208    ENDIF   
    208209
     
    221222       dzu_mg(:,maximum_grid_level) = dzu
    222223!       
    223 !--    To ensure a equally spaced grid. For ocean runs this is not necessary,
     224!--    Next line to ensure an equally spaced grid. For ocean runs this is not
     225!--    necessary,
    224226!--    because zu(0) does not changed so far. Also this would cause errors
    225227!--    if a vertical stretching for ocean runs is used.
    226228       IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2)
     229
    227230       dzw_mg(:,maximum_grid_level) = dzw
    228231       nzt_l = nzt
     
    11321135
    11331136!
    1134 !-- Need to set lateral boundary conditions for l_wall
    1135 
     1137!-- Set lateral boundary conditions for l_wall
    11361138    CALL exchange_horiz( l_wall, nbgp )
    11371139
  • palm/trunk/SOURCE/init_pegrid.f90

    r708 r709  
    44! Current revisions:
    55! -----------------
    6 !
     6! formatting adjustments
    77!
    88! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
     
    577577
    578578          CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
    579 !
    580 !--       TEST OUTPUT (TO BE REMOVED)
    581           WRITE(9,*)  TRIM( coupling_mode ),  &
    582                ', ierr after MPI_OPEN_PORT: ', ierr
    583           CALL LOCAL_FLUSH( 9 )
    584579
    585580          CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
    586581                                 ierr )
    587 !
    588 !--       TEST OUTPUT (TO BE REMOVED)
    589           WRITE(9,*)  TRIM( coupling_mode ),  &
    590                ', ierr after MPI_PUBLISH_NAME: ', ierr
    591           CALL LOCAL_FLUSH( 9 )
    592582
    593583!
     
    614604
    615605          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
    616 !
    617 !--       TEST OUTPUT (TO BE REMOVED)
    618           WRITE(9,*)  TRIM( coupling_mode ),  &
    619                ', ierr after MPI_LOOKUP_NAME: ', ierr
    620           CALL LOCAL_FLUSH( 9 )
    621 
    622606
    623607       ENDIF
     
    631615    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
    632616
    633        PRINT*, '... before COMM_ACCEPT'
    634617       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
    635618                             comm_inter, ierr )
    636        PRINT*, '--- ierr = ', ierr
    637        PRINT*, '--- comm_inter atmosphere = ', comm_inter
    638 
    639619       coupling_mode_remote = 'ocean_to_atmosphere'
    640620
    641621    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
    642622
    643        IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
    644        PRINT*, '... before COMM_CONNECT'
    645623       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
    646624                              comm_inter, ierr )
    647        PRINT*, '--- ierr = ', ierr
    648        PRINT*, '--- comm_inter ocean      = ', comm_inter
    649 
    650625       coupling_mode_remote = 'atmosphere_to_ocean'
    651626
     
    654629
    655630!
    656 !-- Determine the number of ghost points
    657     IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN
     631!-- Determine the number of ghost point layers
     632    IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) THEN
    658633       nbgp = 3
    659634    ELSE
    660635       nbgp = 1
    661     END IF
    662 
    663 !
    664 !-- In case of coupled runs, create a new MPI derived datatype for the
    665 !-- exchange of surface (xy) data .
    666 !-- Gridpoint number for the exchange of ghost points (xy-plane)
    667 
     636    ENDIF
     637
     638!
     639!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
     640!-- which is needed for coupled atmosphere-ocean runs.
     641!-- First, calculate number of grid points of an xy-plane.
    668642    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
    669 
    670 !
    671 !-- Define a new MPI derived datatype for the exchange of ghost points in
    672 !-- y-direction for 2D-arrays (line)
    673643    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
    674644    CALL MPI_TYPE_COMMIT( type_xy, ierr )
    675645
    676 
    677     IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN
     646    IF ( TRIM( coupling_mode ) /= 'uncoupled' )  THEN
    678647   
    679648!
     
    685654          ny_a = ny
    686655
    687           IF ( myid == 0 ) THEN
    688              CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &
    689                             comm_inter, ierr )
    690              CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &
    691                             comm_inter, ierr )
    692              CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &
    693                             comm_inter, ierr )
    694              CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &
    695                             comm_inter, status, ierr )
    696              CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &
    697                             comm_inter, status, ierr )
    698              CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
     656          IF ( myid == 0 )  THEN
     657
     658             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter,  &
     659                            ierr )
     660             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter,  &
     661                            ierr )
     662             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
     663                            ierr )
     664             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter,  &
     665                            status, ierr )
     666             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter,  &
     667                            status, ierr )
     668             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6,      &
    699669                            comm_inter, status, ierr )
    700670          ENDIF
    701671
    702           CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)
    703           CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr)
    704           CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
     672          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
     673          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
     674          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
    705675       
    706676       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
     
    710680
    711681          IF ( myid == 0 ) THEN
    712              CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &
    713                             comm_inter, status, ierr )
    714              CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &
    715                             comm_inter, status, ierr )
    716              CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &
    717                             comm_inter, status, ierr )
    718              CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &
    719                             comm_inter, ierr )
    720              CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &
    721                             comm_inter, ierr )
    722              CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &
    723                             comm_inter, ierr )
     682
     683             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
     684                            ierr )
     685             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
     686                            ierr )
     687             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
     688                            status, ierr )
     689             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
     690             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
     691             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
    724692          ENDIF
    725693
     
    730698       ENDIF
    731699 
    732        ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)
    733        ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)
    734 
    735 !
    736 !--    determine if the horizontal grid and the number of PEs
    737 !--    in ocean and atmosphere is same or not
    738 !--    (different number of PEs still not implemented)
    739        IF ( nx_o == nx_a .AND. ny_o == ny_a .AND.  &
     700       ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
     701       ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
     702
     703!
     704!--    Determine if the horizontal grid and the number of PEs in ocean and
     705!--    atmosphere is same or not
     706       IF ( nx_o == nx_a  .AND.  ny_o == ny_a  .AND.  &
    740707            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
    741708       THEN
     
    748715!--    Determine the target PEs for the exchange between ocean and
    749716!--    atmosphere (comm2d)
    750        IF ( coupling_topology == 0) THEN
    751           IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
     717       IF ( coupling_topology == 0 )  THEN
     718!
     719!--       In case of identical topologies, every atmosphere PE has exactly one
     720!--       ocean PE counterpart and vice versa
     721          IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
    752722             target_id = myid + numprocs
    753723          ELSE
     
    756726
    757727       ELSE
    758 
    759728!
    760729!--       In case of nonequivalent topology in ocean and atmosphere only for
    761730!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
    762 !--       data echxchange between ocean and atmosphere will be done only by
    763 !--       those PEs.   
    764           IF ( myid == 0 ) THEN
    765              IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
     731!--       data echxchange between ocean and atmosphere will be done only
     732!--       between these PEs.   
     733          IF ( myid == 0 )  THEN
     734
     735             IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' )  THEN
    766736                target_id = numprocs
    767737             ELSE
    768738                target_id = 0
    769739             ENDIF
    770  print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs
     740
    771741          ENDIF
     742
    772743       ENDIF
    773744
     
    861832!
    862833!--    Find out, if the total domain allows more levels. These additional
    863 !--    levels are processed on PE0 only.
     834!--    levels are identically processed on all PEs.
    864835       IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
     836
    865837          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
     838
    866839             mg_switch_to_pe0_level_l = maximum_grid_level
    867840
     
    889862                mg_switch_to_pe0_level_l = 0
    890863             ENDIF
     864
    891865          ELSE
     866
    892867             mg_switch_to_pe0_level_l = 0
    893868             maximum_grid_level_l = maximum_grid_level
     869
    894870          ENDIF
    895871
     
    920896
    921897             ENDIF
     898
    922899          ENDIF
    923900
     
    939916!--          Save the grid size of the subdomain at the switch level, because
    940917!--          it is needed in poismg.
    941 !--          Array bounds of the local subdomain grids are gathered on PE0
    942918             ind(1) = nxl_l; ind(2) = nxr_l
    943919             ind(3) = nys_l; ind(4) = nyn_l
     
    953929             DEALLOCATE( ind_all )
    954930!
    955 !--          Calculate the grid size of the total domain gathered on PE0
     931!--          Calculate the grid size of the total domain
    956932             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
    957933             nxl_l = 0
     
    1006982
    1007983!
    1008 !-- Define a new MPI derived datatype for the exchange of ghost points in
    1009 !-- y-direction for 2D-arrays (line)
    1010     CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, ierr )
     984!-- Define new MPI derived datatypes for the exchange of ghost points in
     985!-- x- and y-direction for 2D-arrays (line)
     986    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
     987                          ierr )
    1011988    CALL MPI_TYPE_COMMIT( type_x, ierr )
    1012     CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, type_x_int, ierr )
     989    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
     990                          type_x_int, ierr )
    1013991    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
    1014992
     
    10291007
    10301008    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
     1009
    10311010!
    10321011!-- Discern between the model grid, which needs nbgp ghost points and
    10331012!-- grid levels for the multigrid scheme. In the latter case only one
    10341013!-- ghost point is necessary.
    1035 !-- First definition of mpi-vectors for exchange of ghost layers on normal
     1014!-- First definition of MPI-datatypes for exchange of ghost layers on normal
    10361015!-- grid. The following loop is needed for data exchange in poismg.f90.
    10371016!
    10381017!-- Determine number of grid points of yz-layer for exchange
    10391018    ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
    1040 !
    1041 !-- Define a new mpi datatype for the exchange of left - right boundaries.
    1042 !-- Indeed the data are connected in the physical memory and no mpi-vector
    1043 !-- is necessary, but the data exchange between left and right PE's using
    1044 !-- mpi-vectors is 10% faster than without.
     1019
     1020!
     1021!-- Define an MPI-datatype for the exchange of left/right boundaries.
     1022!-- Although data are contiguous in physical memory (which does not
     1023!-- necessarily require an MPI-derived datatype), the data exchange between
     1024!-- left and right PE's using the MPI-derived type is 10% faster than without.
    10451025    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
    1046                              MPI_REAL, type_xz(0), ierr )
     1026                          MPI_REAL, type_xz(0), ierr )
    10471027    CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
    10481028
    1049     CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), ierr)
     1029    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
     1030                          ierr )
    10501031    CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
    1051 !
    1052 !-- Definition of mpi-vectors for multigrid
     1032
     1033!
     1034!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
    10531035    IF ( psolver == 'multigrid' )  THEN
    10541036!   
    1055 !--   The definition of mpi-vectors as aforementioned, but only 1 ghost point is used.
    1056        DO i = maximum_grid_level, 1 , -1
     1037!--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
     1038       DO  i = maximum_grid_level, 1 , -1
     1039
    10571040          ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
    10581041
    10591042          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
    1060                              MPI_REAL, type_xz(i), ierr )
     1043                                MPI_REAL, type_xz(i), ierr )
    10611044          CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
    10621045
    1063           CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), ierr)
     1046          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
     1047                                ierr )
    10641048          CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
    10651049
     
    10691053          nyn_l = nyn_l / 2
    10701054          nzt_l = nzt_l / 2
     1055
    10711056       ENDDO
    1072     END IF
     1057
     1058    ENDIF
    10731059#endif
    10741060
  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r668 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
     
    1011!
    1112! 667 2010-12-23 12:06:00Z suehring/gryschka
    12 ! Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
    13 ! therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
     13! Changed surface boundary conditions for u and v from mirror to Dirichlet.
     14! Therefore u(uzb,:,:) and v(nzb,:,:) are now representative for height z0.
    1415! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    1516!
  • palm/trunk/SOURCE/pres.f90

    r708 r709  
    44! Current revisions:
    55! -----------------
    6 !
     6! formatting adjustments
    77!
    88! Former revisions:
     
    104104
    105105    ddt_3d = 1.0 / dt_3d
    106     d_weight_pres = 1. / weight_pres(intermediate_timestep_count)
     106    d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
    107107
    108108!
     
    146146!
    147147!-- Left/right
    148 
    149     IF ( conserve_volume_flow  .AND.  ( outflow_l  .OR. outflow_r ) )  THEN
     148    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
    150149
    151150       volume_flow(1)   = 0.0
     
    161160!
    162161!--       Sum up the volume flow through the south/north boundary
    163           DO  k = nzb_2d(j,i) + 1, nzt
     162          DO  k = nzb_2d(j,i)+1, nzt
    164163             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
    165164          ENDDO
     
    173172       volume_flow = volume_flow_l 
    174173#endif
    175        volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
     174       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
    176175                               / volume_flow_area(1)
    177176
    178177       DO  j = nysg, nyng
    179           DO  k = nzb_2d(j,i) + 1, nzt
     178          DO  k = nzb_2d(j,i)+1, nzt
    180179             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
    181180          ENDDO
     
    186185!
    187186!-- South/north
    188     IF ( conserve_volume_flow  .AND.  ( outflow_n  .OR. outflow_s ) )  THEN
     187    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
    189188
    190189       volume_flow(2)   = 0.0
     
    200199!
    201200!--       Sum up the volume flow through the south/north boundary
    202           DO  k = nzb_2d(j,i) + 1, nzt
     201          DO  k = nzb_2d(j,i)+1, nzt
    203202             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
    204203          ENDDO
     
    216215
    217216       DO  i = nxlg, nxrg
    218           DO  k = nzb_v_inner(j,i) + 1, nzt
     217          DO  k = nzb_v_inner(j,i)+1, nzt
    219218             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
    220219          ENDDO
     
    226225!-- Remove mean vertical velocity
    227226    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    228        IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner is not yet known
     227       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
    229228          w_l = 0.0;  w_l_l = 0.0
    230229          DO  i = nxl, nxr
     
    237236#if defined( __parallel )   
    238237          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    239           CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
    240                               ierr )
     238          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
     239                              comm2d, ierr )
    241240#else
    242241          w_l = w_l_l 
     
    574573!-- Correction of the provisional velocities with the current perturbation
    575574!-- pressure just computed
    576     IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .OR. bc_ns_cyc ) )  THEN
     575    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
    577576       volume_flow_l(1) = 0.0
    578577       volume_flow_l(2) = 0.0
  • palm/trunk/SOURCE/prognostic_equations.f90

    r674 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
     
    152153    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
    153154    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
    154     IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
    155        intermediate_timestep_count == 1 )  CALL ws_statistics
     155    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND. &
     156         intermediate_timestep_count == 1 )  CALL ws_statistics
    156157
    157158!
     
    927928    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
    928929    IF ( .NOT. constant_diffusion )  CALL production_e_init
    929     IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
    930        intermediate_timestep_count == 1 )  CALL ws_statistics
     930    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND. &
     931         intermediate_timestep_count == 1 )  CALL ws_statistics
    931932
    932933
     
    14591460    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
    14601461    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
    1461     IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
    1462        intermediate_timestep_count == 1 )  CALL ws_statistics
     1462    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND. &
     1463         intermediate_timestep_count == 1 )  CALL ws_statistics
    14631464
    14641465!
  • palm/trunk/SOURCE/surface_coupler.f90

    r668 r709  
    44! Current revisions:
    55! -----------------
     6! formatting adjustments
    67!
    78! Former revisions:
     
    1011!
    1112! 667 2010-12-23 12:06:00Z suehring/gryschka
    12 ! additional case for nonequivalent processor and grid topopolgy in ocean and
    13 ! atmosphere added (coupling_topology = 1)
     13! Additional case for nonequivalent processor and grid topopolgy in ocean and
     14! atmosphere added (coupling_topology = 1).
    1415! Added exchange of u and v from Ocean to Atmosphere
    1516!
     
    6061
    6162    IF ( coupling_topology == 0 ) THEN
    62        CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id,  &
    63                           0,                                                    &
    64                           terminate_coupled_remote, 1, MPI_INTEGER, target_id,  &
     63       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id, &
     64                          0,                                                   &
     65                          terminate_coupled_remote, 1, MPI_INTEGER, target_id, &
    6566                          0, comm_inter, status, ierr )
    6667    ELSE
     
    7273                             comm_inter, status, ierr )
    7374       ENDIF
    74        CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
     75       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
     76                       ierr )
    7577
    7678       ALLOCATE( total_2d_a(-nbgp:ny_a+nbgp,-nbgp:nx_a+nbgp),       &
     
    9799!-- Exchange the current simulated time between the models,
    98100!-- currently just for total_2ding
    99     IF ( coupling_topology == 0 ) THEN   
    100        CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, &
    101                       target_id, 11, comm_inter, ierr )
    102        CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
    103                       target_id, 11, comm_inter, status, ierr )
     101    IF ( coupling_topology == 0 ) THEN
     102   
     103       CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, 11, &
     104                      comm_inter, ierr )
     105       CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, target_id, &
     106                      11, comm_inter, status, ierr )
    104107    ELSE
     108
    105109       IF ( myid == 0 ) THEN
    106           CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, &
    107                          target_id, 11, comm_inter, ierr )
    108           CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
     110
     111          CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
     112                         11, comm_inter, ierr )
     113          CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
    109114                         target_id, 11, comm_inter, status, ierr )
     115
    110116       ENDIF
    111        CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, &
    112                        0, comm2d, ierr )
    113     ENDIF
    114     WRITE ( 9, * ) 'simulated time: ', simulated_time
    115     WRITE ( 9, * ) 'time since start of coupling: ', &
    116                   time_since_reference_point, ' remote: ', &
    117                   time_since_reference_point_rem
    118    CALL local_flush( 9 )
    119  
     117
     118       CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
     119                       ierr )
     120
     121    ENDIF
    120122
    121123!
     
    124126   
    125127!
    126 !--    Horizontal grid size and number of processors is equal
    127 !--    in ocean and atmosphere
    128        IF ( coupling_topology == 0 ) THEN
    129 
    130 !
    131 !--       Send heat flux at bottom surface to the ocean model
    132           CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, &
    133                          target_id, 12, comm_inter, ierr )
    134 
    135 !
    136 !--       Send humidity flux at bottom surface to the ocean model
     128!--    Horizontal grid size and number of processors is equal in ocean and
     129!--    atmosphere
     130       IF ( coupling_topology == 0 )  THEN
     131
     132!
     133!--       Send heat flux at bottom surface to the ocean
     134          CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
     135                         comm_inter, ierr )
     136!
     137!--       Send humidity flux at bottom surface to the ocean
    137138          IF ( humidity )  THEN
    138              CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, &
    139                             target_id, 13, comm_inter, ierr )
    140           ENDIF
    141 
    142 !
    143 !--       Receive temperature at the bottom surface from the ocean model
    144           WRITE ( 9, * )  '*** receive pt from ocean'
    145           CALL local_flush( 9 )
    146           CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, &
    147                          target_id, 14, comm_inter, status, ierr )
    148 
    149 !
    150 !--       Send the momentum flux (u) at bottom surface to the ocean model
    151           CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, &
    152                          target_id, 15, comm_inter, ierr )
    153 
    154 !
    155 !--       Send the momentum flux (v) at bottom surface to the ocean model
    156           CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, &
    157                          target_id, 16, comm_inter, ierr )
    158 
    159 !
    160 !--       Receive u at the bottom surface from the ocean model
    161           CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, &
    162                          target_id, 17, comm_inter, status, ierr )
    163 
    164 !
    165 !--       Receive v at the bottom surface from the ocean model
    166           CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, &
    167                          target_id, 18,  comm_inter, status, ierr )
    168 
     139             CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 13, &
     140                            comm_inter, ierr )
     141          ENDIF
     142!
     143!--       Receive temperature at the bottom surface from the ocean
     144          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14, &
     145                         comm_inter, status, ierr )
     146!
     147!--       Send the momentum flux (u) at bottom surface to the ocean
     148          CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
     149                         comm_inter, ierr )
     150!
     151!--       Send the momentum flux (v) at bottom surface to the ocean
     152          CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
     153                         comm_inter, ierr )
     154!
     155!--       Receive u at the bottom surface from the ocean
     156          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17, &
     157                         comm_inter, status, ierr )
     158!
     159!--       Receive v at the bottom surface from the ocean
     160          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18, &
     161                         comm_inter, status, ierr )
    169162!
    170163!--    Horizontal grid size or number of processors differs between
     
    173166     
    174167!
    175 !--       Send heat flux at bottom surface to the ocean model
     168!--       Send heat flux at bottom surface to the ocean
    176169          total_2d_a = 0.0
    177           total_2d = 0.0
     170          total_2d   = 0.0
    178171          total_2d(nys:nyn,nxl:nxr) = shf(nys:nyn,nxl:nxr)
    179           CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
    180                            MPI_SUM, 0, comm2d, ierr )
    181           CALL interpolate_to_ocean(12)
    182    
    183 !
    184 !--       Send humidity flux at bottom surface to the ocean model
    185           IF ( humidity ) THEN
     172
     173          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
     174                           comm2d, ierr )
     175          CALL interpolate_to_ocean( 12 )   
     176!
     177!--       Send humidity flux at bottom surface to the ocean
     178          IF ( humidity )  THEN
    186179             total_2d_a = 0.0
    187              total_2d = 0.0
     180             total_2d   = 0.0
    188181             total_2d(nys:nyn,nxl:nxr) = qsws(nys:nyn,nxl:nxr)
    189              CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
    190                               MPI_SUM, 0, comm2d, ierr )
    191              CALL interpolate_to_ocean(13)
    192           ENDIF
    193 
    194 !
    195 !--       Receive temperature at the bottom surface from the ocean model
    196           IF ( myid == 0 ) THEN
     182
     183             CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, &
     184                              0, comm2d, ierr )
     185             CALL interpolate_to_ocean( 13 )
     186          ENDIF
     187!
     188!--       Receive temperature at the bottom surface from the ocean
     189          IF ( myid == 0 )  THEN
    197190             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    198191                            target_id, 14, comm_inter, status, ierr )   
    199192          ENDIF
    200193          CALL MPI_BARRIER( comm2d, ierr )
    201           CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    202                           0, comm2d, ierr )
     194          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
     195                          ierr )
    203196          pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
    204 
    205 !
    206 !--       Send momentum flux (u) at bottom surface to the ocean model
     197!
     198!--       Send momentum flux (u) at bottom surface to the ocean
    207199          total_2d_a = 0.0
    208           total_2d = 0.0
     200          total_2d   = 0.0
    209201          total_2d(nys:nyn,nxl:nxr) = usws(nys:nyn,nxl:nxr)
    210           CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
    211                            MPI_SUM, 0, comm2d, ierr )
    212           CALL interpolate_to_ocean(15)
    213 
    214 !
    215 !--       Send momentum flux (v) at bottom surface to the ocean model
     202          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
     203                           comm2d, ierr )
     204          CALL interpolate_to_ocean( 15 )
     205!
     206!--       Send momentum flux (v) at bottom surface to the ocean
    216207          total_2d_a = 0.0
    217           total_2d = 0.0
     208          total_2d   = 0.0
    218209          total_2d(nys:nyn,nxl:nxr) = vsws(nys:nyn,nxl:nxr)
    219           CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, &
    220                            MPI_SUM, 0, comm2d, ierr )
    221           CALL interpolate_to_ocean(16)
    222 
    223 !
    224 !--       Receive u at the bottom surface from the ocean model
    225           IF ( myid == 0 ) THEN
     210          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
     211                           comm2d, ierr )
     212          CALL interpolate_to_ocean( 16 )
     213!
     214!--       Receive u at the bottom surface from the ocean
     215          IF ( myid == 0 )  THEN
    226216             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    227                             target_id, 17, comm_inter, status, ierr )           
     217                            target_id, 17, comm_inter, status, ierr )
    228218          ENDIF
    229219          CALL MPI_BARRIER( comm2d, ierr )
    230           CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    231                           0, comm2d, ierr )
     220          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
     221                          ierr )
    232222          u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
    233    
    234 !
    235 !--       Receive v at the bottom surface from the ocean model
    236           IF ( myid == 0 ) THEN
     223!
     224!--       Receive v at the bottom surface from the ocean
     225          IF ( myid == 0 )  THEN
    237226             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    238                             target_id, 18, comm_inter, status, ierr )           
     227                            target_id, 18, comm_inter, status, ierr )
    239228          ENDIF
    240229          CALL MPI_BARRIER( comm2d, ierr )
    241           CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    242                           0, comm2d, ierr )
     230          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
     231                          ierr )
    243232          v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
    244233
     
    252241       IF ( coupling_topology == 0 ) THEN
    253242!
    254 !--       Receive heat flux at the sea surface (top) from the atmosphere model
    255           CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, &
    256                          target_id, 12, comm_inter, status, ierr )
    257 
    258 
    259 !
    260 !--       Receive humidity flux from the atmosphere model (bottom)
     243!--       Receive heat flux at the sea surface (top) from the atmosphere
     244          CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
     245                         comm_inter, status, ierr )
     246!
     247!--       Receive humidity flux from the atmosphere (bottom)
    261248!--       and add it to the heat flux at the sea surface (top)...
    262249          IF ( humidity_remote )  THEN
    263250             CALL MPI_RECV( qswst_remote(nysg,nxlg), ngp_xy, MPI_REAL, &
    264251                            target_id, 13, comm_inter, status, ierr )
    265 
    266           ENDIF
    267 
     252          ENDIF
    268253!
    269254!--       Send sea surface temperature to the atmosphere model
    270           CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, &
    271                          target_id, 14, comm_inter, ierr )
    272 
     255          CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, target_id, 14, &
     256                         comm_inter, ierr )
    273257!
    274258!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
    275 !--       model
    276           WRITE ( 9, * )  '*** receive uswst from atmosphere'
    277           CALL local_flush( 9 )
    278           CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, &
    279                          target_id, 15, comm_inter, status, ierr )
    280 
     259          CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
     260                         comm_inter, status, ierr )
    281261!
    282262!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
    283 !--       model
    284           CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, &
    285                          target_id, 16, comm_inter, status, ierr )
    286 
    287 !--       Send u to the atmosphere model
    288           CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, &
    289                          target_id, 17, comm_inter, ierr )
    290 
    291 !
    292 !--       Send v to the atmosphere model
    293           CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, &
    294                          target_id, 18, comm_inter, ierr )
    295 
     263          CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
     264                         comm_inter, status, ierr )
     265!
     266!--       Send u to the atmosphere
     267          CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, target_id, 17, &
     268                         comm_inter, ierr )
     269!
     270!--       Send v to the atmosphere
     271          CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, target_id, 18, &
     272                         comm_inter, ierr )
    296273!
    297274!--    Horizontal gridsize or number of processors differs between
    298275!--    ocean and atmosphere
    299276       ELSE
    300 
    301 !
    302 !--       Receive heat flux at the sea surface (top) from the atmosphere model
    303           IF ( myid == 0 ) THEN
     277!
     278!--       Receive heat flux at the sea surface (top) from the atmosphere
     279          IF ( myid == 0 )  THEN
    304280             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    305                             target_id, 12, comm_inter, status, ierr )           
     281                            target_id, 12, comm_inter, status, ierr )
     282          ENDIF
     283          CALL MPI_BARRIER( comm2d, ierr )
     284          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
     285                          ierr )
     286          tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
     287!
     288!--       Receive humidity flux at the sea surface (top) from the atmosphere
     289          IF ( humidity_remote )  THEN
     290             IF ( myid == 0 )  THEN
     291                CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
     292                               target_id, 13, comm_inter, status, ierr )
     293             ENDIF
     294             CALL MPI_BARRIER( comm2d, ierr )
     295             CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, &
     296                             comm2d, ierr)
     297             qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
     298          ENDIF
     299!
     300!--       Send surface temperature to atmosphere
     301          total_2d_o = 0.0
     302          total_2d   = 0.0
     303          total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
     304
     305          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
     306                           comm2d, ierr)
     307          CALL interpolate_to_atmos( 14 )
     308!
     309!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
     310          IF ( myid == 0 )  THEN
     311             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
     312                            target_id, 15, comm_inter, status, ierr )
    306313          ENDIF
    307314          CALL MPI_BARRIER( comm2d, ierr )
    308315          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    309                           0, comm2d, ierr)
    310           tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
    311 
    312 !
    313 !--       Receive humidity flux at the sea surface (top) from the
    314 !--       atmosphere model
    315           IF ( humidity_remote ) THEN
    316              IF ( myid == 0 ) THEN
    317                 CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    318                                target_id, 13, comm_inter, status, ierr )           
    319              ENDIF
    320              CALL MPI_BARRIER( comm2d, ierr )
    321              CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    322                              0, comm2d, ierr)
    323              qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
    324           ENDIF
    325 
    326 !
    327 !--       Send surface temperature to atmosphere
    328           total_2d_o = 0.0
    329           total_2d = 0.0
    330           total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
    331 
    332           CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, &
    333                           MPI_REAL, MPI_SUM, 0, comm2d, ierr)
    334 
    335           CALL interpolate_to_atmos(14)
    336 
    337 !
    338 !--       Receive momentum flux (u) at the sea surface (top) from the
    339 !--       atmosphere model
    340           IF ( myid == 0 ) THEN
     316                          0, comm2d, ierr )
     317          uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
     318!
     319!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
     320          IF ( myid == 0 )  THEN
    341321             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    342                             target_id, 15, comm_inter, status, ierr )           
     322                            target_id, 16, comm_inter, status, ierr )
    343323          ENDIF
    344324          CALL MPI_BARRIER( comm2d, ierr )
    345           CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    346                           0, comm2d, ierr)
    347           uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
    348 
    349 !
    350 !--       Receive momentum flux (v) at the sea surface (top) from the
    351 !--       atmosphere model
    352           IF ( myid == 0 ) THEN
    353              CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    354                             target_id, 16, comm_inter, status, ierr )           
    355           ENDIF
    356           CALL MPI_BARRIER( comm2d, ierr )
    357           CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
    358                           0, comm2d, ierr)
     325          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
     326                          ierr )
    359327          vswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
    360 
    361328!
    362329!--       Send u to atmosphere
    363330          total_2d_o = 0.0 
    364           total_2d = 0.0
     331          total_2d   = 0.0
    365332          total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr)
    366           CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, &
    367                           MPI_SUM, 0, comm2d, ierr) 
    368           CALL interpolate_to_atmos(17)
    369 
     333          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
     334                           comm2d, ierr )
     335          CALL interpolate_to_atmos( 17 )
    370336!
    371337!--       Send v to atmosphere
    372338          total_2d_o = 0.0
    373           total_2d = 0.0
     339          total_2d   = 0.0
    374340          total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr)
    375           CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, &
    376                           MPI_SUM, 0, comm2d, ierr) 
    377           CALL interpolate_to_atmos(18)
     341          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
     342                           comm2d, ierr )
     343          CALL interpolate_to_atmos( 18 )
    378344       
    379345       ENDIF
     
    382348!--    Conversions of fluxes received from atmosphere
    383349       IF ( humidity_remote )  THEN
    384           !here tswst is still the sum of atmospheric bottom heat fluxes
    385           tswst = tswst + qswst_remote * 2.2626108e6 / 1005.0
    386           !*latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
    387           !/(rho_atm(=1.0)*c_p)
     350!
     351!--       Here tswst is still the sum of atmospheric bottom heat fluxes,
     352!--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
     353!--       /(rho_atm(=1.0)*c_p)
     354          tswst = tswst + qswst_remote * 2.2626108E6 / 1005.0
    388355!
    389356!--        ...and convert it to a salinity flux at the sea surface (top)
     
    405372       vswst = vswst / rho(nzt,:,:)
    406373
    407 
    408     ENDIF
    409 
    410     IF ( coupling_topology == 1 ) THEN
     374    ENDIF
     375
     376    IF ( coupling_topology == 1 )  THEN
    411377       DEALLOCATE( total_2d_o, total_2d_a )
    412378    ENDIF
     
    420386
    421387
    422   SUBROUTINE interpolate_to_atmos(tag)
     388  SUBROUTINE interpolate_to_atmos( tag )
    423389
    424390    USE arrays_3d
     
    430396    IMPLICIT NONE
    431397
    432  
    433398    INTEGER             ::  dnx, dnx2, dny, dny2, i, ii, j, jj
    434399    INTEGER, intent(in) ::  tag
     
    436401    CALL MPI_BARRIER( comm2d, ierr )
    437402
    438     IF ( myid == 0 ) THEN
    439 
    440 !
    441 !--    cyclic boundary conditions for the total 2D-grid
     403    IF ( myid == 0 )  THEN
     404!
     405!--    Cyclic boundary conditions for the total 2D-grid
    442406       total_2d_o(-nbgp:-1,:) = total_2d_o(ny+1-nbgp:ny,:)
    443407       total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx+1-nbgp:nx)
     
    452416
    453417!
    454 !--    Distance for interpolation around coarse grid points within the fine grid
    455 !--    (note: 2*dnx2 must not be equal with dnx) 
     418!--    Distance for interpolation around coarse grid points within the fine
     419!--    grid (note: 2*dnx2 must not be equal with dnx)
    456420       dnx2 = 2 * ( dnx / 2 )
    457421       dny2 = 2 * ( dny / 2 )
     
    472436       ENDDO
    473437!
    474 !--    cyclic boundary conditions for atmosphere grid
     438!--    Cyclic boundary conditions for atmosphere grid
    475439       total_2d_a(-nbgp:-1,:) = total_2d_a(ny_a+1-nbgp:ny_a,:)
    476440       total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx_a+1-nbgp:nx_a)
     
    480444!
    481445!--    Transfer of the atmosphere-grid-layer to the atmosphere
    482        CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
    483                       target_id, tag, comm_inter, ierr )
     446       CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, target_id, &
     447                      tag, comm_inter, ierr )
    484448
    485449    ENDIF
     
    490454
    491455
    492   SUBROUTINE interpolate_to_ocean(tag)
     456  SUBROUTINE interpolate_to_ocean( tag )
    493457
    494458    USE arrays_3d
     
    500464    IMPLICIT NONE
    501465
    502     REAL                ::  fl, fr, myl, myr
    503466    INTEGER             ::  dnx, dny, i, ii, j, jj
    504467    INTEGER, intent(in) ::  tag
     468    REAL                ::  fl, fr, myl, myr
     469
    505470
    506471    CALL MPI_BARRIER( comm2d, ierr )
    507472
    508     IF ( myid == 0 ) THEN   
    509 
    510 !
    511 !      Number of gridpoints of the fine grid within one mesh of the coarse grid
     473    IF ( myid == 0 )  THEN   
     474
     475!
     476!--    Number of gridpoints of the fine grid within one mesh of the coarse grid
    512477       dnx = ( nx_o + 1 ) / ( nx_a + 1 )
    513478       dny = ( ny_o + 1 ) / ( ny_a + 1 )
    514479
    515480!
    516 !--    cyclic boundary conditions for atmosphere grid
     481!--    Cyclic boundary conditions for atmosphere grid
    517482       total_2d_a(-nbgp:-1,:) = total_2d_a(ny+1-nbgp:ny,:)
    518483       total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx+1-nbgp:nx)
     
    521486       total_2d_a(:,nx+1:nx+nbgp) = total_2d_a(:,0:nbgp-1)
    522487!
    523 !--    Bilinear Interpolation from atmosphere-grid-layer to ocean-grid-layer
     488!--    Bilinear Interpolation from atmosphere grid-layer to ocean grid-layer
    524489       DO  j = 0, ny
    525490          DO  i = 0, nx
     
    527492             myr = ( total_2d_a(j+1,i+1) - total_2d_a(j,i+1) ) / dny
    528493             DO  jj = 0, dny-1
    529                 fl = myl*jj  + total_2d_a(j,i) 
    530                 fr = myr*jj  + total_2d_a(j,i+1) 
     494                fl = myl*jj + total_2d_a(j,i) 
     495                fr = myr*jj + total_2d_a(j,i+1) 
    531496                DO  ii = 0, dnx-1
    532497                   total_2d_o(j*dny+jj,i*dnx+ii) = ( fr - fl ) / dnx * ii + fl
     
    536501       ENDDO
    537502!
    538 !--    cyclic boundary conditions for ocean grid
     503!--    Cyclic boundary conditions for ocean grid
    539504       total_2d_o(-nbgp:-1,:) = total_2d_o(ny_o+1-nbgp:ny_o,:)
    540505       total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx_o+1-nbgp:nx_o)
     
    542507       total_2d_o(ny_o+1:ny_o+nbgp,:) = total_2d_o(0:nbgp-1,:)
    543508       total_2d_o(:,nx_o+1:nx_o+nbgp) = total_2d_o(:,0:nbgp-1)
    544        
    545509
    546510       CALL MPI_SEND( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
Note: See TracChangeset for help on using the changeset viewer.