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

preliminary changes concerning update of BC-scheme

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.