Changeset 63


Ignore:
Timestamp:
Mar 13, 2007 3:52:49 AM (18 years ago)
Author:
raasch
Message:

preliminary changes concerning update of BC-scheme

Location:
palm/trunk
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/.mrun.config.default

    r40 r63  
    5959%cpumax            1000                                       ibmb parallel
    6060%numprocs          4                                          ibmb parallel
    61 #%tmp_data_catalog  /fastfs/work/<replace by your HLRN username>/palm_restart_data    ibmh parallel debug
     61#%tmp_data_catalog  /fastfs/work/<replace by your HLRN username>/palm_restart_data    ibmb parallel debug
    6262#
    6363#%remote_username   <replace by your DKRZ username>           nech parallel
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r61 r63  
    99
    1010new d3par-parameter dt_max to define the maximum value for the allowed timestep
     11
     12new inipar-parameter loop_optimization to control the loop optimization method
    1113
    1214new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).
  • palm/trunk/SOURCE/advec_particles.f90

    r60 r63  
    23152315          ENDIF
    23162316
    2317           CALL MPI_SENDRECV( trlp, trlp_count, mpi_particle_type, pleft, 1,    &
     2317          CALL MPI_SENDRECV( trlp(1), trlp_count, mpi_particle_type, pleft, 1, &
    23182318                             particles(number_of_particles+1), trrp_count_recv,&
    23192319                             mpi_particle_type, pright, 1,                     &
     
    23432343             ENDIF
    23442344
    2345              CALL MPI_SENDRECV( trlpt, trlpt_count*tlength, MPI_REAL, pleft, 1,&
     2345             CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,   &
     2346                                pleft, 1,                                      &
    23462347                             particle_tail_coordinates(1,1,number_of_tails+1), &
    23472348                                trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
     
    23912392          ENDIF
    23922393
    2393           CALL MPI_SENDRECV( trrp, trrp_count, mpi_particle_type, pright, 1,   &
     2394          CALL MPI_SENDRECV( trrp(1), trrp_count, mpi_particle_type, pright, 1,&
    23942395                             particles(number_of_particles+1), trlp_count_recv,&
    23952396                             mpi_particle_type, pleft, 1,                      &
     
    24192420             ENDIF
    24202421
    2421              CALL MPI_SENDRECV( trrpt, trrpt_count*tlength, MPI_REAL, pright,  &
    2422                           1, particle_tail_coordinates(1,1,number_of_tails+1), &
     2422             CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,   &
     2423                                pright, 1,                                     &
     2424                             particle_tail_coordinates(1,1,number_of_tails+1), &
    24232425                                trlpt_count_recv*tlength, MPI_REAL, pleft, 1,  &
    24242426                                comm2d, status, ierr )
     
    27542756          ENDIF
    27552757
    2756           CALL MPI_SENDRECV( trsp, trsp_count, mpi_particle_type, psouth, 1,   &
     2758          CALL MPI_SENDRECV( trsp(1), trsp_count, mpi_particle_type, psouth, 1,&
    27572759                             particles(number_of_particles+1), trnp_count_recv,&
    27582760                             mpi_particle_type, pnorth, 1,                     &
     
    27822784             ENDIF
    27832785
    2784              CALL MPI_SENDRECV( trspt, trspt_count*tlength, MPI_REAL, psouth,  &
    2785                           1, particle_tail_coordinates(1,1,number_of_tails+1), &
     2786             CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,   &
     2787                                psouth, 1,                                     &
     2788                             particle_tail_coordinates(1,1,number_of_tails+1), &
    27862789                                trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
    27872790                                comm2d, status, ierr )
     
    28302833          ENDIF
    28312834
    2832           CALL MPI_SENDRECV( trnp, trnp_count, mpi_particle_type, pnorth, 1,   &
     2835          CALL MPI_SENDRECV( trnp(1), trnp_count, mpi_particle_type, pnorth, 1,&
    28332836                             particles(number_of_particles+1), trsp_count_recv,&
    28342837                             mpi_particle_type, psouth, 1,                     &
     
    28582861             ENDIF
    28592862
    2860              CALL MPI_SENDRECV( trnpt, trnpt_count*tlength, MPI_REAL, pnorth,  &
    2861                           1, particle_tail_coordinates(1,1,number_of_tails+1), &
     2863             CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,   &
     2864                                pnorth, 1,                                     &
     2865                             particle_tail_coordinates(1,1,number_of_tails+1), &
    28622866                                trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
    28632867                                comm2d, status, ierr )
  • palm/trunk/SOURCE/advec_s_bc.f90

    r4 r63  
    11 SUBROUTINE advec_s_bc( sk, sk_char )
    22
    3 !-------------------------------------------------------------------------------!
     3!------------------------------------------------------------------------------!
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
     
    2929! stretched grids.
    3030! NOTE: This is a provisional, non-optimised version!
    31 !-------------------------------------------------------------------------------!
     31!------------------------------------------------------------------------------!
    3232
    3333    USE advection
    3434    USE arrays_3d
     35    USE control_parameters
    3536    USE cpulog
    3637    USE grid_variables
     
    3940    USE pegrid
    4041    USE statistics
    41     USE control_parameters
    4242
    4343    IMPLICIT NONE
     
    5757    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  sk_p
    5858
    59 #if defined( __t3eb ) || defined( __t3eh ) || defined( __t3ej2 ) || defined( __t3ej5 )
     59#if defined( __nec )
    6060    REAL (kind=4) ::  m1n, m1z  !Wichtig: Division
    6161    REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw
     
    6868!
    6969!-- Array sk_p requires 2 extra elements for each dimension
    70     ALLOCATE( sk_p(nzb-2:nzt+2,nys-3:nyn+3,nxl-3:nxr+3) )
     70    ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
    7171    sk_p = 0.0
    7272
     
    8989    DO  i = nxl-1, nxr+1
    9090       DO  j = nys, nyn
    91           DO  k = nzb+1, nzt-1
     91          DO  k = nzb, nzt+1
    9292             sk_p(k,j,i) = sk(k,j,i)
    9393          ENDDO
     
    9595    ENDDO
    9696#if defined( __parallel )
    97     ngp = 2 * ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
     97    ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
    9898    CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
    9999!
     
    142142    DO  i = nxl, nxr
    143143       DO  j = nys, nyn
    144           DO  k = nzb+1, nzt-1
     144          DO  k = nzb+1, nzt
    145145             zaehler = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
    146146             nenner  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     
    160160!
    161161!-- Allocate temporary arrays
    162     ALLOCATE( a0(nzb+1:nzt-1,nxl-1:nxr+1),   a1(nzb+1:nzt-1,nxl-1:nxr+1),   &
    163               a2(nzb+1:nzt-1,nxl-1:nxr+1),   a12(nzb+1:nzt-1,nxl-1:nxr+1),  &
    164               a22(nzb+1:nzt-1,nxl-1:nxr+1),  immb(nzb+1:nzt-1,nxl-1:nxr+1), &
    165               imme(nzb+1:nzt-1,nxl-1:nxr+1), impb(nzb+1:nzt-1,nxl-1:nxr+1), &
    166               impe(nzb+1:nzt-1,nxl-1:nxr+1), ipmb(nzb+1:nzt-1,nxl-1:nxr+1), &
    167               ipme(nzb+1:nzt-1,nxl-1:nxr+1), ippb(nzb+1:nzt-1,nxl-1:nxr+1), &
    168               ippe(nzb+1:nzt-1,nxl-1:nxr+1), m1(nzb+1:nzt-1,nxl-2:nxr+2),   &
    169               sw(nzb+1:nzt-1,nxl-1:nxr+1)                                   &
     162    ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),   &
     163              a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),  &
     164              a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1), &
     165              imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &
     166              impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &
     167              ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &
     168              ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),   &
     169              sw(nzb+1:nzt,nxl-1:nxr+1)                                 &
    170170            )
    171171    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     
    184184!--    Compute polynomial coefficients
    185185       DO  i = nxl-1, nxr+1
    186           DO  k = nzb+1, nzt-1
     186          DO  k = nzb+1, nzt
    187187             a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    188188             a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) &
     
    204204!--    *VOCL LOOP,UNROLL(2)
    205205       DO  i = nxl, nxr
    206           DO  k = nzb+1, nzt-1
     206          DO  k = nzb+1, nzt
    207207             cip  =  MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    208208             cim  = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
     
    240240!--    Compute monitor function m1
    241241       DO  i = nxl-2, nxr+2
    242           DO  k = nzb+1, nzt-1
     242          DO  k = nzb+1, nzt
    243243             m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
    244244             m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     
    258258       sw = 0.0
    259259       DO  i = nxl-1, nxr+1
    260           DO  k = nzb+1, nzt-1
     260          DO  k = nzb+1, nzt
    261261             m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / &
    262262                  MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 )
     
    284284       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    285285       DO  i = nxl, nxr
    286           DO  k = nzb+1, nzt-1
     286          DO  k = nzb+1, nzt
    287287
    288288!--          *VOCL STMT,IF(10)
     
    374374!--    Prognostic equation
    375375       DO  i = nxl, nxr
    376           DO  k = nzb+1, nzt-1
     376          DO  k = nzb+1, nzt
    377377             fplus  = ( 1.0 - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i) &
    378378                    - ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
     
    404404!-- Enlarge boundary of local array cyclically in y-direction
    405405#if defined( __parallel )
    406     ngp = ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
    407     CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+5), ngp, MPI_REAL, &
     406    ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
     407    CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, &
    408408                          type_xz_2, ierr )
    409409    CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
     
    423423#else
    424424    DO  i = nxl, nxr
    425        DO  k = nzb+1, nzt-1
     425       DO  k = nzb+1, nzt
    426426          sk_p(k,nys-1,i) = sk_p(k,nyn,i)
    427427          sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
     
    439439    DO  i = nxl, nxr
    440440       DO  j = nys, nyn
    441           DO  k = nzb+1, nzt-1
     441          DO  k = nzb+1, nzt
    442442             zaehler = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
    443443             nenner  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     
    457457!
    458458!-- Allocate temporary arrays
    459     ALLOCATE( a0(nzb+1:nzt-1,nys-1:nyn+1),   a1(nzb+1:nzt-1,nys-1:nyn+1),   &
    460               a2(nzb+1:nzt-1,nys-1:nyn+1),   a12(nzb+1:nzt-1,nys-1:nyn+1),  &
    461               a22(nzb+1:nzt-1,nys-1:nyn+1),  immb(nzb+1:nzt-1,nys-1:nyn+1), &
    462               imme(nzb+1:nzt-1,nys-1:nyn+1), impb(nzb+1:nzt-1,nys-1:nyn+1), &
    463               impe(nzb+1:nzt-1,nys-1:nyn+1), ipmb(nzb+1:nzt-1,nys-1:nyn+1), &
    464               ipme(nzb+1:nzt-1,nys-1:nyn+1), ippb(nzb+1:nzt-1,nys-1:nyn+1), &
    465               ippe(nzb+1:nzt-1,nys-1:nyn+1), m1(nzb+1:nzt-1,nys-2:nyn+2),   &
    466               sw(nzb+1:nzt-1,nys-1:nyn+1)                                   &
     459    ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),   &
     460              a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),  &
     461              a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1), &
     462              imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), &
     463              impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), &
     464              ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), &
     465              ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),   &
     466              sw(nzb+1:nzt,nys-1:nyn+1)                                 &
    467467            )
    468468    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     
    475475!--    Compute polynomial coefficients
    476476       DO  j = nys-1, nyn+1
    477           DO  k = nzb+1, nzt-1
     477          DO  k = nzb+1, nzt
    478478             a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    479479             a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) &
     
    495495!--    *VOCL LOOP,UNROLL(2)
    496496       DO  j = nys, nyn
    497           DO  k = nzb+1, nzt-1
     497          DO  k = nzb+1, nzt
    498498             cip  =  MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    499499             cim  = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
     
    531531!--    Compute monitor function m1
    532532       DO  j = nys-2, nyn+2
    533           DO  k = nzb+1, nzt-1
     533          DO  k = nzb+1, nzt
    534534             m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
    535535             m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     
    549549       sw = 0.0
    550550       DO  j = nys-1, nyn+1
    551           DO  k = nzb+1, nzt-1
     551          DO  k = nzb+1, nzt
    552552             m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
    553553                  MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
     
    575575       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    576576       DO  j = nys, nyn
    577           DO  k = nzb+1, nzt-1
     577          DO  k = nzb+1, nzt
    578578
    579579!--          *VOCL STMT,IF(10)
     
    665665!--    Prognostic equation
    666666       DO  j = nys, nyn
    667           DO  k = nzb+1, nzt-1
     667          DO  k = nzb+1, nzt
    668668             fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    669669                    - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
     
    720720          DO  i = nxl, nxr
    721721             DO  j = nys, nyn
    722                 sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    723                 sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
    724                 sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
     722                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     723                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    725724             ENDDO
    726725          ENDDO
     
    730729!
    731730!--    Temperature boundary condition at the top boundary
    732        IF ( ibc_pt_t == 0 )  THEN
    733 !
    734 !--       Dirichlet
     731       IF ( ibc_pt_t == 0  .OR.  ibc_pt_t == 1 )  THEN
     732!
     733!--       Dirichlet or Neumann (zero gradient)
    735734          DO  i = nxl, nxr
    736735             DO  j = nys, nyn
    737                 sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
    738                 sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
     736                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     737                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
    739738             ENDDO
    740739          ENDDO
    741740
    742        ELSE
    743 !
    744 !--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
     741       ELSEIF ( ibc_pt_t == 2 )  THEN
     742!
     743!--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
    745744          DO  i = nxl, nxr
    746745             DO  j = nys, nyn
    747                 sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_pt_t_val * dzu(nzt)
    748                 sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_pt_t_val * dzu(nzt+1)
    749746                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
     747                sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
    750748             ENDDO
    751749          ENDDO
     
    756754
    757755!
    758 !--    Temperature boundary condition at the bottom boundary
     756!--    Specific humidity boundary condition at the bottom boundary
    759757       IF ( ibc_q_b == 0 )  THEN
    760758!
    761 !--       Dirichlet (fixed surface temperature)
     759!--       Dirichlet (fixed surface humidity)
    762760          DO  i = nxl, nxr
    763761             DO  j = nys, nyn
     
    772770          DO  i = nxl, nxr
    773771             DO  j = nys, nyn
    774                 sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    775                 sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
    776                 sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
     772                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     773                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    777774             ENDDO
    778775          ENDDO
     
    781778
    782779!
    783 !--    Temperature boundary condition at the top boundary
     780!--    Specific humidity boundary condition at the top boundary
    784781       IF ( ibc_q_t == 0 )  THEN
    785782!
     
    787784          DO  i = nxl, nxr
    788785             DO  j = nys, nyn
    789                 sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
    790                 sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
     786                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     787                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
    791788             ENDDO
    792789          ENDDO
     
    794791       ELSE
    795792!
    796 !--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
     793!--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
    797794          DO  i = nxl, nxr
    798795             DO  j = nys, nyn
    799                 sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_q_t_val * dzu(nzt)
    800                 sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_q_t_val * dzu(nzt+1)
    801796                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
     797                sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
    802798             ENDDO
    803799          ENDDO
     
    809805!
    810806!--    TKE boundary condition at bottom and top boundary (generally Neumann)
    811        sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    812        sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
    813        sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
    814        sk_p(nzt,j,i)   = sk_p(nzt-1,j,i)
    815        sk_p(nzt+1,j,i) = sk_p(nzt-1,j,i)
    816        sk_p(nzt+2,j,i) = sk_p(nzt-1,j,i)
     807       sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     808       sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     809       sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
     810       sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
    817811
    818812    ELSE
     
    829823    DO  i = nxl, nxr
    830824       DO  j = nys, nyn
    831           DO  k = nzb, nzt
     825          DO  k = nzb, nzt+1
    832826             zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
    833827             nenner  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
     
    847841!
    848842!-- Allocate temporary arrays
    849     ALLOCATE( a0(nzb:nzt,nys:nyn),       a1(nzb:nzt,nys:nyn),       &
    850               a2(nzb:nzt,nys:nyn),       a12(nzb:nzt,nys:nyn),      &
    851               a22(nzb:nzt,nys:nyn),      immb(nzb+1:nzt-1,nys:nyn), &
    852               imme(nzb+1:nzt-1,nys:nyn), impb(nzb+1:nzt-1,nys:nyn), &
    853               impe(nzb+1:nzt-1,nys:nyn), ipmb(nzb+1:nzt-1,nys:nyn), &
    854               ipme(nzb+1:nzt-1,nys:nyn), ippb(nzb+1:nzt-1,nys:nyn), &
    855               ippe(nzb+1:nzt-1,nys:nyn), m1(nzb-1:nzt+1,nys:nyn),  &
    856               sw(nzb:nzt,nys:nyn)                                   &
     843    ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),   &
     844              a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),  &
     845              a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn), &
     846              imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), &
     847              impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), &
     848              ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), &
     849              ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), &
     850              sw(nzb:nzt+1,nys:nyn)                             &
    857851            )
    858852    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     
    865859!--    Compute polynomial coefficients
    866860       DO  j = nys, nyn
    867           DO  k = nzb, nzt
     861          DO  k = nzb, nzt+1
    868862             a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    869863             a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
     
    885879!--    *VOCL LOOP,UNROLL(2)
    886880       DO  j = nys, nyn
    887           DO  k = nzb+1, nzt-1
     881          DO  k = nzb+1, nzt
    888882             cip  =  MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
    889883             cim  = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
     
    921915!--    Compute monitor function m1
    922916       DO  j = nys, nyn
    923           DO  k = nzb-1, nzt+1
     917          DO  k = nzb-1, nzt+2
    924918             m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
    925919             m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
     
    939933       sw = 0.0
    940934       DO  j = nys, nyn
    941           DO  k = nzb, nzt
     935          DO  k = nzb, nzt+1
    942936             m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
    943937                  MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
     
    965959       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    966960       DO  j = nys, nyn
    967           DO  k = nzb+1, nzt-1
     961          DO  k = nzb+1, nzt
    968962
    969963!--          *VOCL STMT,IF(10)
     
    10551049!--    Prognostic equation
    10561050       DO  j = nys, nyn
    1057           DO  k = nzb+1, nzt-1
     1051          DO  k = nzb+1, nzt
    10581052             fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    10591053                    - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
     
    10821076          DO  sr = 0, statistic_regions
    10831077             DO  j = nys, nyn
    1084                 DO  k = nzb+1, nzt-1
     1078                DO  k = nzb+1, nzt
    10851079                   sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
    10861080                                          m1(k,j) * rmask(j,i,sr)
     
    11031097    DO  i = nxl, nxr
    11041098       DO  j = nys, nyn
    1105           DO  k = nzb+1, nzt-1
     1099          DO  k = nzb+1, nzt
    11061100             tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
    11071101          ENDDO
  • palm/trunk/SOURCE/check_parameters.f90

    r57 r63  
    55! -----------------
    66! "by_user" allowed as initializing action, -data_output_ts,
    7 ! leapfrog with non-flat topography not allowed any more, pt_reference is
    8 ! checked
     7! leapfrog with non-flat topography not allowed any more, loop_optimization
     8! and pt_reference are checked
    99!
    1010! Former revisions:
     
    8282
    8383!
     84!-- Check the general loop optimization method
     85    IF ( loop_optimization == 'default' )  THEN
     86       IF ( host(1:3) == 'nec' )  THEN
     87          loop_optimization = 'vector'
     88       ELSE
     89          loop_optimization = 'cache'
     90       ENDIF
     91    ENDIF
     92    IF ( loop_optimization /= 'noopt'  .AND.  loop_optimization /= 'cache' &
     93         .AND.  loop_optimization /= 'vector' )  THEN
     94       IF ( myid == 0 )  THEN
     95          PRINT*, '+++ check_parameters:'
     96          PRINT*, '    illegal value given for loop_optimization: ', &
     97                  TRIM( loop_optimization )
     98       ENDIF
     99       CALL local_stop
     100    ENDIF
     101
     102!
    84103!-- Check topography setting (check for illegal parameter combinations)
    85104    IF ( topography /= 'flat' )  THEN
     
    257276    END SELECT
    258277
    259     IF ( scalar_advec /= 'pw-scheme'  .AND.  timestep_scheme(1:5) == 'runge' ) &
     278    IF ( scalar_advec == 'ups-scheme'  .AND.  timestep_scheme(1:5) == 'runge' )&
    260279    THEN
    261280       IF ( myid == 0 )  THEN
  • palm/trunk/SOURCE/header.f90

    r60 r63  
    44! Actual revisions:
    55! -----------------
    6 ! Output of netcdf_64bit_3d, particles-package is now part of the default code.
     6! Output of netcdf_64bit_3d, particles-package is now part of the default code,
     7! output of the loop optimization method.
    78!
    89! Former revisions:
     
    189190       WRITE ( io, 118 )
    190191    ENDIF
    191     IF ( momentum_advec /= 'ups-scheme' .AND. scalar_advec /= 'ups-scheme' &
    192          .AND. scalar_advec /= 'bc-scheme'  .AND.  host(1:3) /= 'nec' )  THEN
    193        WRITE ( io, 139 )
    194     ENDIF
     192
     193    WRITE ( io, 139 )  TRIM( loop_optimization )
     194
    195195    IF ( galilei_transformation )  THEN
    196196       IF ( use_ug_for_galilei_tr )  THEN
     
    11821182                  I3,')')
    11831183138 FORMAT ('     Using hybrid version for 1d-domain-decomposition')
    1184 139 FORMAT (' --> Prognostic equations are solved in one single loop (fast', &
    1185                   ' method)')
     1184139 FORMAT (' --> Loop optimization method: ',A)
    11861185140 FORMAT ('     maximum residual allowed:                ',E10.3)
    11871186141 FORMAT ('     fixed number of multigrid cycles:        ',I4)
  • palm/trunk/SOURCE/init_3d_model.f90

    r51 r63  
    756756!
    757757!-- If required, initialize particles
    758     CALL init_particles
     758    IF ( particle_advection )  CALL init_particles
    759759
    760760!
  • palm/trunk/SOURCE/modules.f90

    r57 r63  
    66! -----------------
    77! +array rif_wall
    8 ! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*,
    9 ! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, use_pt_reference,
     8! +loop_optimization, netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*,
     9! id_var_zwwi_*, ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference,
     10! use_pt_reference,
    1011! +age_m in particle_type
    1112! -data_output_ts, dots_n
     
    211212    CHARACTER (LEN=9)   ::  simulated_time_chr
    212213    CHARACTER (LEN=12)  ::  version = 'PALM   3.1c'
    213     CHARACTER (LEN=16)  ::  momentum_advec = 'pw-scheme', &
     214    CHARACTER (LEN=16)  ::  loop_optimization = 'default', &
     215                            momentum_advec = 'pw-scheme', &
    214216                            psolver = 'poisfft', &
    215217                            scalar_advec = 'pw-scheme'
  • palm/trunk/SOURCE/parin.f90

    r61 r63  
    44! Actual revisions:
    55! -----------------
    6 ! +dt_max, netcdf_64bit_3d in d3par, +pt_reference in inipar, -data_output_ts
     6! +dt_max, netcdf_64bit_3d in d3par, +loop_optimization, pt_reference in
     7! inipar, -data_output_ts
    78!
    89! Former revisions:
     
    6061                       inflow_disturbance_end, initializing_actions, &
    6162                       km_constant, km_damp_max, long_filter_factor, &
    62                        mixing_length_1d, moisture, momentum_advec, &
    63                        netcdf_precision, npex, npey, nsor_ini, nx, ny, nz, &
    64                        omega, outflow_damping_width, overshoot_limit_e, &
     63                       loop_optimization, mixing_length_1d, moisture, &
     64                       momentum_advec, netcdf_precision, npex, npey, nsor_ini, &
     65                       nx, ny, nz, omega, outflow_damping_width, &
     66                       overshoot_limit_e, &
    6567                       overshoot_limit_pt, overshoot_limit_u, &
    6668                       overshoot_limit_v, overshoot_limit_w, passive_scalar, &
  • palm/trunk/SOURCE/prognostic_equations.f90

    r57 r63  
    44! Actual revisions:
    55! -----------------
    6 ! z0 removed from arguments in calls of diffusion_u/v/w
     6! z0 removed from arguments in calls of diffusion_u/v/w,
     7! subroutine names changed to .._noopt, .._cache, and .._vector
    78!
    89! Former revisions:
     
    6364
    6465    PRIVATE
    65     PUBLIC prognostic_equations, prognostic_equations_fast, &
    66            prognostic_equations_vec
    67 
    68     INTERFACE prognostic_equations
    69        MODULE PROCEDURE prognostic_equations
    70     END INTERFACE prognostic_equations
    71 
    72     INTERFACE prognostic_equations_fast
    73        MODULE PROCEDURE prognostic_equations_fast
    74     END INTERFACE prognostic_equations_fast
    75 
    76     INTERFACE prognostic_equations_vec
    77        MODULE PROCEDURE prognostic_equations_vec
    78     END INTERFACE prognostic_equations_vec
     66    PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
     67           prognostic_equations_vector
     68
     69    INTERFACE prognostic_equations_noopt
     70       MODULE PROCEDURE prognostic_equations_noopt
     71    END INTERFACE prognostic_equations_noopt
     72
     73    INTERFACE prognostic_equations_cache
     74       MODULE PROCEDURE prognostic_equations_cache
     75    END INTERFACE prognostic_equations_cache
     76
     77    INTERFACE prognostic_equations_vector
     78       MODULE PROCEDURE prognostic_equations_vector
     79    END INTERFACE prognostic_equations_vector
    7980
    8081
     
    8283
    8384
    84  SUBROUTINE prognostic_equations
     85 SUBROUTINE prognostic_equations_noopt
    8586
    8687!------------------------------------------------------------------------------!
     
    615616
    616617
    617  END SUBROUTINE prognostic_equations
    618 
    619 
    620  SUBROUTINE prognostic_equations_fast
     618 END SUBROUTINE prognostic_equations_noopt
     619
     620
     621 SUBROUTINE prognostic_equations_cache
    621622
    622623!------------------------------------------------------------------------------!
     
    996997
    997998
    998  END SUBROUTINE prognostic_equations_fast
    999 
    1000 
    1001  SUBROUTINE prognostic_equations_vec
     999 END SUBROUTINE prognostic_equations_cache
     1000
     1001
     1002 SUBROUTINE prognostic_equations_vector
    10021003
    10031004!------------------------------------------------------------------------------!
     
    12421243!
    12431244!-- pt-tendency terms with communication
     1245    sat = tsc(1)
     1246    sbt = tsc(2)
    12441247    IF ( scalar_advec == 'bc-scheme' )  THEN
    1245 !
    1246 !--    Bott-Chlond scheme always uses Euler time step. Thus:
    1247        sat = 1.0
    1248        sbt = 1.0
     1248
     1249       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1250!
     1251!--       Bott-Chlond scheme always uses Euler time step when leapfrog is
     1252!--       switched on. Thus:
     1253          sat = 1.0
     1254          sbt = 1.0
     1255       ENDIF
    12491256       tend = 0.0
    12501257       CALL advec_s_bc( pt, 'pt' )
    12511258    ELSE
    1252        sat = tsc(1)
    1253        sbt = tsc(2)
    12541259       IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
    12551260          tend = 0.0
     
    13401345!
    13411346!--    Scalar/q-tendency terms with communication
     1347       sat = tsc(1)
     1348       sbt = tsc(2)
    13421349       IF ( scalar_advec == 'bc-scheme' )  THEN
    1343 !
    1344 !--       Bott-Chlond scheme always uses Euler time step. Thus:
    1345           sat = 1.0
    1346           sbt = 1.0
     1350
     1351          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1352!
     1353!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
     1354!--          switched on. Thus:
     1355             sat = 1.0
     1356             sbt = 1.0
     1357          ENDIF
    13471358          tend = 0.0
    13481359          CALL advec_s_bc( q, 'q' )
    13491360       ELSE
    1350           sat = tsc(1)
    1351           sbt = tsc(2)
    13521361          IF ( tsc(2) /= 2.0 )  THEN
    13531362             IF ( scalar_advec == 'ups-scheme' )  THEN
     
    14381447!--    TKE-tendency terms with communication
    14391448       CALL production_e_init
     1449
     1450       sat = tsc(1)
     1451       sbt = tsc(2)
    14401452       IF ( .NOT. use_upstream_for_tke )  THEN
    14411453          IF ( scalar_advec == 'bc-scheme' )  THEN
    1442 !
    1443 !--          Bott-Chlond scheme always uses Euler time step. Thus:
    1444              sat = 1.0
    1445              sbt = 1.0
     1454
     1455             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1456!
     1457!--             Bott-Chlond scheme always uses Euler time step when leapfrog is
     1458!--             switched on. Thus:
     1459                sat = 1.0
     1460                sbt = 1.0
     1461             ENDIF
    14461462             tend = 0.0
    14471463             CALL advec_s_bc( e, 'e' )
    14481464          ELSE
    1449              sat = tsc(1)
    1450              sbt = tsc(2)
    14511465             IF ( tsc(2) /= 2.0 )  THEN
    14521466                IF ( scalar_advec == 'ups-scheme' )  THEN
     
    15501564
    15511565
    1552  END SUBROUTINE prognostic_equations_vec
     1566 END SUBROUTINE prognostic_equations_vector
    15531567
    15541568
  • palm/trunk/SOURCE/read_var_list.f90

    r57 r63  
    44! Actual revisions:
    55! -----------------
    6 ! +pt_reference
     6! +loop_optimization, pt_reference
    77!
    88! Former revisions:
     
    208208          CASE ( 'long_filter_factor' )
    209209             READ ( 13 )  long_filter_factor
     210          CASE ( 'loop_optimization' )
     211             READ ( 13 )  loop_optimization
    210212          CASE ( 'mixing_length_1d' )
    211213             READ ( 13 )  mixing_length_1d
  • palm/trunk/SOURCE/time_integration.f90

    r48 r63  
    55! -----------------
    66! Move call of user_actions( 'after_integration' ) below increment of times
    7 ! and counters
     7! and counters,
     8! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
     9! .._vector, these calls are now controlled by switch loop_optimization
    810!
    911! Former revisions:
     
    98100!--       in the other versions a good vectorization is prohibited due to
    99101!--       inlining problems.
    100           IF ( host(1:3) == 'nec' )  THEN
    101              CALL prognostic_equations_vec
     102          IF ( loop_optimization == 'vector' )  THEN
     103             CALL prognostic_equations_vector
    102104          ELSE
    103105             IF ( momentum_advec == 'ups-scheme'  .OR.  &
     
    105107                  scalar_advec == 'bc-scheme' )        &
    106108             THEN
    107                 CALL prognostic_equations
     109                CALL prognostic_equations_noopt
    108110             ELSE
    109                 CALL prognostic_equations_fast
     111                CALL prognostic_equations_cache
    110112             ENDIF
    111113          ENDIF
     
    114116!--       Particle advection (only once during intermediate steps, because
    115117!--       it uses an Euler-step)
    116           IF ( simulated_time >= particle_advection_start  .AND. &
     118          IF ( particle_advection  .AND.                         &
     119               simulated_time >= particle_advection_start  .AND. &
    117120               intermediate_timestep_count == 1 )  THEN
    118121             CALL advec_particles
  • palm/trunk/SOURCE/write_var_list.f90

    r57 r63  
    44! Actual revisions:
    55! -----------------
    6 ! +pt_refrence
     6! +loop_optimization, pt_refrence
    77!
    88! Former revisions:
     
    174174    WRITE ( 14 )  'long_filter_factor            '
    175175    WRITE ( 14 )  long_filter_factor
     176    WRITE ( 14 )  'loop_optimization             '
     177    WRITE ( 14 )  loop_optimization
    176178    WRITE ( 14 )  'mixing_length_1d              '
    177179    WRITE ( 14 )  mixing_length_1d
Note: See TracChangeset for help on using the changeset viewer.