text:
*> \brief \b SCHKGL
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SCHKGL( NIN, NOUT )
*
*       .. Scalar Arguments ..
*       INTEGER            NIN, NOUT
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SCHKGL tests SGGBAL, a routine for balancing a matrix pair (A, B).
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] NIN
*> \verbatim
*>          NIN is INTEGER
*>          The logical unit number for input.  NIN > 0.
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*>          NOUT is INTEGER
*>          The logical unit number for output.  NOUT > 0.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup single_eig
*
*  =====================================================================
      SUBROUTINE SCHKGL( NIN, NOUT )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            NIN, NOUT
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            LDA, LDB, LWORK
      PARAMETER          ( LDA = 20, LDB = 20, LWORK = 6*LDA )
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IHI, IHIIN, ILO, ILOIN, INFO, J, KNT, N,
     $                   NINFO
      REAL               ANORM, BNORM, EPS, RMAX, VMAX
*     ..
*     .. Local Arrays ..
      INTEGER            LMAX( 5 )
      REAL               A( LDA, LDA ), AIN( LDA, LDA ), B( LDB, LDB ),
     $                   BIN( LDB, LDB ), LSCALE( LDA ), LSCLIN( LDA ),
     $                   RSCALE( LDA ), RSCLIN( LDA ), WORK( LWORK )
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLANGE
      EXTERNAL           SLAMCH, SLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGGBAL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     ..
*     .. Executable Statements ..
*
      LMAX( 1 ) = 0
      LMAX( 2 ) = 0
      LMAX( 3 ) = 0
      NINFO = 0
      KNT = 0
      RMAX = ZERO
*
      EPS = SLAMCH( 'Precision' )
*
   10 CONTINUE
*
      READ( NIN, FMT = * )N
      IF( N.EQ.0 )
     $   GO TO 90
      DO 20 I = 1, N
         READ( NIN, FMT = * )( A( I, J ), J = 1, N )
   20 CONTINUE
*
      DO 30 I = 1, N
         READ( NIN, FMT = * )( B( I, J ), J = 1, N )
   30 CONTINUE
*
      READ( NIN, FMT = * )ILOIN, IHIIN
      DO 40 I = 1, N
         READ( NIN, FMT = * )( AIN( I, J ), J = 1, N )
   40 CONTINUE
      DO 50 I = 1, N
         READ( NIN, FMT = * )( BIN( I, J ), J = 1, N )
   50 CONTINUE
*
      READ( NIN, FMT = * )( LSCLIN( I ), I = 1, N )
      READ( NIN, FMT = * )( RSCLIN( I ), I = 1, N )
*
      ANORM = SLANGE( 'M', N, N, A, LDA, WORK )
      BNORM = SLANGE( 'M', N, N, B, LDB, WORK )
*
      KNT = KNT + 1
*
      CALL SGGBAL( 'B', N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
     $             WORK, INFO )
*
      IF( INFO.NE.0 ) THEN
         NINFO = NINFO + 1
         LMAX( 1 ) = KNT
      END IF
*
      IF( ILO.NE.ILOIN .OR. IHI.NE.IHIIN ) THEN
         NINFO = NINFO + 1
         LMAX( 2 ) = KNT
      END IF
*
      VMAX = ZERO
      DO 70 I = 1, N
         DO 60 J = 1, N
            VMAX = MAX( VMAX, ABS( A( I, J )-AIN( I, J ) ) )
            VMAX = MAX( VMAX, ABS( B( I, J )-BIN( I, J ) ) )
   60    CONTINUE
   70 CONTINUE
*
      DO 80 I = 1, N
         VMAX = MAX( VMAX, ABS( LSCALE( I )-LSCLIN( I ) ) )
         VMAX = MAX( VMAX, ABS( RSCALE( I )-RSCLIN( I ) ) )
   80 CONTINUE
*
      VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
*
      IF( VMAX.GT.RMAX ) THEN
         LMAX( 3 ) = KNT
         RMAX = VMAX
      END IF
*
      GO TO 10
*
   90 CONTINUE
*
      WRITE( NOUT, FMT = 9999 )
 9999 FORMAT( 1X, '.. test output of SGGBAL .. ' )
*
      WRITE( NOUT, FMT = 9998 )RMAX
 9998 FORMAT( 1X, 'value of largest test error            = ', E12.3 )
      WRITE( NOUT, FMT = 9997 )LMAX( 1 )
 9997 FORMAT( 1X, 'example number where info is not zero  = ', I4 )
      WRITE( NOUT, FMT = 9996 )LMAX( 2 )
 9996 FORMAT( 1X, 'example number where ILO or IHI wrong  = ', I4 )
      WRITE( NOUT, FMT = 9995 )LMAX( 3 )
 9995 FORMAT( 1X, 'example number having largest error    = ', I4 )
      WRITE( NOUT, FMT = 9994 )NINFO
 9994 FORMAT( 1X, 'number of examples where info is not 0 = ', I4 )
      WRITE( NOUT, FMT = 9993 )KNT
 9993 FORMAT( 1X, 'total number of examples tested        = ', I4 )
*
      RETURN
*
*     End of SCHKGL
*
      END

max_stars_repo_name: sydha/Lapack

max_issues_repo_name: sydha/Lapack

max_forks_repo_name: sydha/Lapack


text:
!BSD 3-Clause License
!
!Copyright (c) 2021, Kiyle Winborne
!All rights reserved.
!
!Redistribution and use in source and binary forms, with or without
!modification, are permitted provided that the following conditions are met:
!
!1. Redistributions of source code must retain the above copyright notice, this
!   list of conditions and the following disclaimer.
!
!2. Redistributions in binary form must reproduce the above copyright notice,
!   this list of conditions and the following disclaimer in the documentation
!   and/or other materials provided with the distribution.
!
!3. Neither the name of the copyright holder nor the names of its
!   contributors may be used to endorse or promote products derived from
!   this software without specific prior written permission.
!
!THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
!DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
!FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
!DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
!SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
!CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
!OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
!OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! Unit name: List_Module
! Author: Kiyle Winborne
! Date: 3/2/2021
!
! Dependencies: List_Node_Module
!
! Purpose: Defines a singly linked list data structure.
! Procedure, Type-bound: Append - Add a node to the end of the list.
!       Argument types: List, List_Node
!
! Procedure, Type-bound: get_next - Returns the next node in the list, starting from the head node.
!       Argument types: List, Integer
!
! Procedure, Type-bound: free_list - Removes each node in the list and deallocates the memory.
!       Argument types: List
!
! Procedure, Type-bound: Remove - Remove a node at any position in the list.
!       Argument types: List, Integer

module List_Module
    use List_Node_Module
    implicit none

    !========================================================================
    !Definition of a singly linked list data structure.
    !Each list consists of a head node, defined in List_Node_Module
    !and a tail node. Structure also contains a list length for processing.
    !========================================================================
    type, public :: List
        integer :: length_of = 0
        type(List_Node), pointer :: head
        type(List_Node), pointer :: tail
    contains
        procedure, pass(this) :: append
        procedure, pass(this) :: remove_head
        procedure, pass(this) :: remove
        procedure, pass(this) :: remove_tail
        procedure, pass(this) :: free_list
        procedure, pass(this) :: get_next
    end type List

    contains
    !=============================
    !  Append a node to the list
    !=============================
    subroutine append (this, element)
        class(List), intent(in out) :: this
        class(*), intent(in) :: element
        if(associated (this%tail)) then
            allocate(this%tail%next, source = new_node(element, null()))
            this%tail => this%tail%next
            this%tail%List_Position = this%length_of+1
        else
            allocate(this%head, source = new_node(element, null()))
            this%tail => this%head
            this%head%list_position = this%length_of +1
        end if
        this%length_of = this%length_of + 1
    end subroutine append
    !================================
    ! Remove all nodes from the list
    !================================
    subroutine free_list (this)
        class(List), intent(inout) :: this
        do while(this%length_of /= 0)
            call this%remove_tail()
        end do
    end subroutine free_list
    !===================================
    ! Return the next node in the list
    !===================================
    function get_next(this, node_position) result (node)
        class(List), intent(in) :: this
        integer, intent(in) :: node_position
        type(List_Node) :: node
        type(List_Node), pointer :: current_node
        integer :: i
        current_node => this%head
        do i = 2, node_position
        current_node => current_node%next
            end do
            node = current_node
    end function get_next

    !===================================
    ! Remove a single node from the list
    !===================================
    subroutine remove (this, node_position)
        class(List), intent(inout) :: this
        integer, intent(in) :: node_position
        type(List_Node) :: node
        type(List_Node), pointer :: current_node
        type(List_Node), pointer :: temp_node
        integer :: i
        current_node => this%head
        if(node_position == 1) then
            call this%remove_head()
            return
        else if(node_position == this%length_of) then
            call this%remove_tail()
            return
        end if
        do i = 2, node_position
            if(i == node_position) then
                temp_node => current_node
            end if
        current_node => current_node%next
        end do
        temp_node%next => current_node%next
        deallocate(current_node%element)
        current_node%next => null()
        this%length_of = this%length_of - 1
    end subroutine remove

        subroutine remove_head (this)
        class(list), intent(inout) :: this
        type(List_Node) :: node_to_remove
        node_to_remove = this%head
        this%head = this%head%next
        deallocate(node_to_remove%element)
        node_to_remove%next => null()
        this%length_of = this%length_of -1
    end subroutine remove_head

    subroutine remove_tail (this)
        class(list), intent(inout) :: this
        type(List_Node) :: node_to_remove
        node_to_remove = this%tail
        deallocate(node_to_remove%element)
        this%length_of = this%length_of - 1
    end subroutine remove_tail

end module List_Module

max_stars_repo_path: Singly Linked List/src/list.f90

max_stars_repo_name: kwinborne/Fortran-Structures

max_issues_repo_path: Singly Linked List/src/list.f90

max_issues_repo_name: kwinborne/Fortran-Structures

max_forks_repo_path: Singly Linked List/src/list.f90

max_forks_repo_name: kwinborne/Fortran-Structures


text:
! @@name:	affinity.3f
! @@type:	F-fixed
! @@compilable:	yes
! @@linkable:	yes
! @@expect:	success
      PROGRAM EXAMPLE
!$OMP PARALLEL PROC_BIND(CLOSE) NUM_THREADS(4)
      CALL WORK()
!$OMP END PARALLEL
      END PROGRAM EXAMPLE

ext: f

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
1 FORMAT (1PE12.4, I10)
2 FORMAT (I12, /, ' Dates: ', 2 (2I3, I5))
3 FORMAT (' Enter the node number: ', $ )
end

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
MODULE functions

  ABSTRACT INTERFACE
      SUBROUTINE pr (a)
         INTEGER, INTENT(IN) :: a
      END SUBROUTINE pr
  END INTERFACE

  !Use a type to hold the function pointer so we can have an array of them
  TYPE f_ptr
    PROCEDURE (pr), pointer, nopass :: ptr
  END TYPE f_ptr

  CONTAINS

  SUBROUTINE print_i(a)
    INTEGER, INTENT(IN) :: a
    PRINT*, a
  END SUBROUTINE

  SUBROUTINE print_fizz(a)
    INTEGER, INTENT(IN) :: a
    PRINT*, "fizz"
  END SUBROUTINE

  SUBROUTINE print_buzz(a)
    INTEGER, INTENT(IN) :: a
    PRINT*, "buzz"
  END SUBROUTINE

  SUBROUTINE print_fizzbuzz(a)
    INTEGER, INTENT(IN) :: a
    PRINT*, "fizzbuzz"
  END SUBROUTINE

END MODULE

PROGRAM main

  USE functions
  IMPLICIT NONE

  TYPE(f_ptr), DIMENSION(4) :: printers
  TYPE(f_ptr) :: printer
  INTEGER :: i, offset
  INTEGER, PARAMETER :: en = 100, fizz = 3, buzz = 5

  printers(1)%ptr => print_i
  printers(2)%ptr => print_fizz
  printers(3)%ptr => print_buzz
  printers(4)%ptr => print_fizzbuzz

  DO i = 1, en

    !Work out the offset we need
    ! We want '+1' if MOD(i, fizz) == 0
    ! and then '+2' if MOD(i, buzz) == 0
    ! In Fortran we can't manipulate logicals, so we use 'MIN'
    ! This looks horrible, and it is. This is bad cleverness.

    offset =  MAX(0, 1+MIN(0, -MOD(i, fizz))) + &
        2 * MAX(0, 1+MIN(0, -MOD(i, buzz)))

    printer%ptr => printers(1+offset)%ptr
    CALL printer%ptr(i)

  END DO

END PROGRAM

max_stars_repo_name: WarwickRSE/SOUP

max_issues_repo_name: WarwickRSE/SOUP

max_forks_repo_name: WarwickRSE/SOUP


text:
      SUBROUTINE CALQQ2T (ISTL_)
C
C CHANGE RECORD
C  FIXED DYNAMIC TIME STEPPING
C **  SUBROUTINE CALQQ CALCULATES THE TURBULENT INTENSITY SQUARED AT
C **  TIME LEVEL (N+1).  THE VALUE OF ISTL INDICATES THE NUMBER OF
C **  TIME LEVELS INVOLVED
C
      USE GLOBAL
!$      USE OMP_LIB
      REAL,PARAMETER::BETAVEG_P=1.0, BETAVEG_D=5.1,CE4VEG=0.9 !from Katul et al. 2003
      REAL,PARAMETER::BETASUP_P=1.0, BETASUP_D=5.1,CE4SUP=0.9 !from Katul et al. 2003
      REAL,SAVE,ALLOCATABLE::PQQVEGI(:,:),PQQVEGE(:,:)
      REAL,SAVE,ALLOCATABLE::PQQMHKI(:,:),PQQMHKE(:,:)
      REAL,SAVE,ALLOCATABLE::PQQSUPI(:,:),PQQSUPE(:,:)
      REAL,SAVE,ALLOCATABLE::PQQDRPI(:,:),PQQDRPE(:,:)
      REAL::TMPQQI,TMPQQE
      IF(.NOT.ALLOCATED(PQQVEGI)) THEN
        ALLOCATE(PQQVEGI(LCM,KCM))
        PQQVEGI=0.0
        ALLOCATE(PQQVEGE(LCM,KCM))
        PQQVEGE=0.0
        ALLOCATE(PQQMHKI(LCM,KCM))
        PQQMHKI=0.0
        ALLOCATE(PQQMHKE(LCM,KCM))
        PQQMHKE=0.0
        ALLOCATE(PQQSUPI(LCM,KCM))
        PQQSUPI=0.0
        ALLOCATE(PQQSUPE(LCM,KCM))
        PQQSUPE=0.0
        ALLOCATE(PQQDRPI(LCM,KCM))
        PQQDRPI=0.0
        ALLOCATE(PQQDRPE(LCM,KCM))
        PQQDRPE=0.0
      END IF
      IF(ISDYNSTP.EQ.0)THEN
        DELT=DT
        DELTD2=0.5*DT
        DELTI=1./DELT
      ELSE
        DELT=DTDYN
        DELTD2=0.5*DTDYN
        DELTI=1./DELT
      END IF
      S2TL=0.0
      BSMALL=1.E-12
      DO K=1,KS
        DO L=2,LA
          QQ2(L,K)=QQ(L,K)+QQ(L,K)
          QQL2(L,K)=QQL(L,K)+QQL(L,K)
        ENDDO
      ENDDO
C

!$OMP PARALLEL PRIVATE(LS,LW,WB,UHUW,VHVW)
      DO K = 1,KC
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
        DO L=1,LC
          UUU(L,K)=0.
          VVV(L,K)=0.
          FUHU(L,K)=0.
          FUHV(L,K)=0.
          FVHU(L,K)=0.
          FUHV(L,K)=0.
          PQQVEGI(L,K)=0.
          PQQVEGE(L,K)=0.
          PQQMHKI(L,K)=0.
          PQQMHKE(L,K)=0.
          PQQSUPI(L,K)=0.
          PQQSUPE(L,K)=0.
          PQQDRPI(L,K)=0.
          PQQDRPE(L,K)=0.
        END DO
!$OMP END DO NOWAIT
      ENDDO

C
C **  CALCULATE ADVECTIVE FLUXES BY UPWIND DIFFERENCE WITH TRANSPORT
C **  AVERAGED BETWEEN (N) AND (N+1) AND TRANSPORTED FIELD AT (N) OR
C **  TRANSPORT BETWEEN (N-1) AND (N+1) AND TRANSPORTED FIELD AT (N-1)
C **  FOR ISTL EQUAL TO 2 AND 3 RESPECTIVELY
C **  FUHQQ=FUHU, FVHQQ=FVHU, FUHQQL=FUHV, FVHQQL=FVHV
C
      IF(IDRYTBP.EQ.0)THEN
        DO K=1,KC
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA
            WB=0.5*DXYP(L)*(W2(L,K-1)+W2(L,K))
            FWQQ(L,K)=MAX(WB,0.)*QQ(L,K-1)
     &          +MIN(WB,0.)*QQ(L,K)
            FWQQL(L,K)=MAX(WB,0.)*QQL(L,K-1)*H1P(L)
     &          +MIN(WB,0.)*QQL(L,K)*H1P(L)
          ENDDO
!$OM END DO NOWAIT
        ENDDO
      ELSE
        DO K=1,KC
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              WB=0.5*DXYP(L)*(W2(L,K-1)+W2(L,K))
              FWQQ(L,K)=MAX(WB,0.)*QQ(L,K-1)
     &            +MIN(WB,0.)*QQ(L,K)
              FWQQL(L,K)=MAX(WB,0.)*QQL(L,K-1)*H1P(L)
     &            +MIN(WB,0.)*QQL(L,K)*H1P(L)
            ELSE
              FWQQ(L,K)=0.
              FWQQL(L,K)=0.
            ENDIF
          ENDDO
!$OMP END DO NOWAIT
        ENDDO
      ENDIF

C
      IF(IDRYTBP.EQ.0)THEN
        DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA
            LS=LSC(L)
            LW=LWEST(L)
            UHUW=0.5*(UHDY2(L,K)+UHDY2(L,K+1))
            VHVW=0.5*(VHDX2(L,K)+VHDX2(L,K+1))
            FUHU(L,K)=MAX(UHUW,0.)*QQ(LW,K)
     &               +MIN(UHUW,0.)*QQ(L ,K)
            FVHU(L,K)=MAX(VHVW,0.)*QQ(LS,K)
     &               +MIN(VHVW,0.)*QQ(L ,K)
            FUHV(L,K)=MAX(UHUW,0.)*QQL(LW,K)*H1P(LW)
     &               +MIN(UHUW,0.)*QQL(L ,K)*H1P(L )
            FVHV(L,K)=MAX(VHVW,0.)*QQL(LS,K)*H1P(LS)
     &               +MIN(VHVW,0.)*QQL(L ,K)*H1P(L )
          ENDDO
!$OMP END DO NOWAIT
        ENDDO
      ELSE
        DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              LS=LSC(L)
              LW=LWEST(L)
              UHUW=0.5*(UHDY2(L,K)+UHDY2(L,K+1))
              VHVW=0.5*(VHDX2(L,K)+VHDX2(L,K+1))
              FUHU(L,K)=MAX(UHUW,0.)*QQ(LW,K)
     &                 +MIN(UHUW,0.)*QQ(L ,K)
              FVHU(L,K)=MAX(VHVW,0.)*QQ(LS,K)
     &                 +MIN(VHVW,0.)*QQ(L ,K)
              FUHV(L,K)=MAX(UHUW,0.)*QQL(LW,K)*H1P(LW)
     &                 +MIN(UHUW,0.)*QQL(L ,K)*H1P(L )
              FVHV(L,K)=MAX(VHVW,0.)*QQL(LS,K)*H1P(LS)
     &                 +MIN(VHVW,0.)*QQL(L ,K)*H1P(L )
            ELSE
              FUHU(L,K)=0.
              FUHV(L,K)=0.
              FVHU(L,K)=0.
              FUHV(L,K)=0.
            ENDIF
          ENDDO
!$OMP END DO NOWAIT
        ENDDO
      ENDIF
!$OMP END PARALLEL

C
C **  CALCULATE PRODUCTION, LOAD BOUNDARY CONDITIONS AND SOLVE
C **  TRANSPORT EQUATIONS
C **  FUHQQ=FUHU, FVHQQ=FVHU, FUHQQL=FUHV, FVHQQL=FVHV
C **  CU1=CUQ, CU2=CUQL, UUU=QQH, VVV=QQLH
C
      DO K=1,KC
        DO LL=1,NCBS
          L=LCBS(LL)
          LN=LNC(L)
          IF(FVHU(LN,K).GT.0)THEN
            FVHU(LN,K)=0.0
            FVHV(LN,K)=0.0
          ENDIF
        ENDDO
        DO LL=1,NCBW
          L=LCBW(LL)
          LE=LEAST(L)
          IF(FUHU(LE,K).GT.0)THEN
            FUHU(LE,K)=0.0
            FUHV(LE,K)=0.0
          ENDIF
        ENDDO
        DO LL=1,NCBE
          L=LCBE(LL)
          IF(FUHU(L,K).LT.0.)THEN
            FUHU(L,K)=0.0
            FUHV(L,K)=0.0
          ENDIF
        ENDDO
        DO LL=1,NCBN
          L=LCBN(LL)
          IF(FVHU(L,K).LT.0.)THEN
            FVHU(L,K)=0.0
            FVHV(L,K)=0.0
          ENDIF
        ENDDO
      ENDDO
!      WHERE(FVHU(LCBN(1:NCBN),1:KC)<0.0)
!        FVHU(LCBN(1:NCBN),1:KC)=0.0
!        FVHV(LCBN(1:NCBN),1:KC)=0.0
!      ENDWHERE
      FVHU(LCBN(1:NCMB),1:KC)=MAX(FVHU(LCBN(1:NCMB),1:KC),0.0)
      FVHV(LCBN(1:NCMB),1:KC)=MAX(FVHV(LCBN(1:NCMB),1:KC),0.0)
      IF(ISVEG.GT.0)THEN !SCJ vegetative/MHK impact on K-epsilon
        DO K=1,KS
          DO L=2,LA
            LE=LEAST(L)
            LN=LNC(L)
            TMPQQI=0.25*BETAVEG_P
            TMPQQE=0.25*BETAVEG_D
            PQQVEGI(L,K)=TMPQQI*(FXVEG(L ,K  )+FXVEG(L ,K+1)
     &                          +FXVEG(LE,K  )+FXVEG(LE,K+1)
     &                          +FYVEG(L ,K  )+FYVEG(L ,K+1)
     &                          +FYVEG(LN,K  )+FYVEG(LN,K+1))
            PQQVEGE(L,K)=TMPQQE*(FXVEG(L ,K  )*U(L ,K  )*U(L ,K  )
     &                          +FXVEG(L ,K+1)*U(L ,K+1)*U(L ,K+1)
     &                          +FXVEG(LE,K  )*U(LE,K  )*U(LE,K  )
     &                          +FXVEG(LE,K+1)*U(LE,K+1)*U(LE,K+1)
     &                          +FYVEG(L ,K  )*V(L ,K  )*V(L ,K  )
     &                          +FYVEG(L ,K+1)*V(L ,K+1)*V(L ,K+1)
     &                          +FYVEG(LN,K  )*V(LN,K  )*V(LN,K  )
     &                          +FYVEG(LN,K+1)*V(LN,K+1)*V(LN,K+1))
            IF(MVEGL(L)>90)THEN
              TMPQQI=0.25*BETAMHK_P
              TMPQQE=0.25*BETAMHK_D
              PQQMHKI(L,K)=TMPQQI*(FXMHK(L ,K  )+FXMHK(L ,K+1)
     &                            +FXMHK(LE,K  )+FXMHK(LE,K+1)
     &                            +FYMHK(L ,K  )+FYMHK(L ,K+1)
     &                            +FYMHK(LN,K  )+FYMHK(LN,K+1))
              PQQMHKE(L,K)=TMPQQE*(FXMHK(L ,K  )*U(L ,K  )*U(L ,K  )
     &                            +FXMHK(L ,K+1)*U(L ,K+1)*U(L ,K+1)
     &                            +FXMHK(LE,K  )*U(LE,K  )*U(LE,K  )
     &                            +FXMHK(LE,K+1)*U(LE,K+1)*U(LE,K+1)
     &                            +FYMHK(L ,K  )*V(L ,K  )*V(L ,K  )
     &                            +FYMHK(L ,K+1)*V(L ,K+1)*V(L ,K+1)
     &                            +FYMHK(LN,K  )*V(LN,K  )*V(LN,K  )
     &                            +FYMHK(LN,K+1)*V(LN,K+1)*V(LN,K+1))
              TMPQQI=0.25*BETASUP_P
              TMPQQE=0.25*BETASUP_D
              PQQSUPI(L,K)=TMPQQI*(FXSUP(L ,K  )+FXSUP(L ,K+1)
     &                            +FXSUP(LE,K  )+FXSUP(LE,K+1)
     &                            +FYSUP(L ,K  )+FYSUP(L ,K+1)
     &                            +FYSUP(LN,K  )+FYSUP(LN,K+1))
              PQQSUPE(L,K)=TMPQQE*(FXSUP(L ,K  )*U(L ,K  )*U(L ,K  )
     &                            +FXSUP(L ,K+1)*U(L ,K+1)*U(L ,K+1)
     &                            +FXSUP(LE,K  )*U(LE,K  )*U(LE,K  )
     &                            +FXSUP(LE,K+1)*U(LE,K+1)*U(LE,K+1)
     &                            +FYSUP(L ,K  )*V(L ,K  )*V(L ,K  )
     &                            +FYSUP(L ,K+1)*V(L ,K+1)*V(L ,K+1)
     &                            +FYSUP(LN,K  )*V(LN,K  )*V(LN,K  )
     &                            +FYSUP(LN,K+1)*V(LN,K+1)*V(LN,K+1))
            ENDIF
          ENDDO
        ENDDO
      ENDIF
      IF(AQUAINCL==1 .AND. AQUATURB==1)THEN !FOD suspended canopy impact on K-epsilon; based on SCJ version
        DO K=1,KS
          DO L=2,LA
            LE=L+1
            LN=LNC(L)
            TMPQQI=0.25*BETADRP_P
            TMPQQE=0.25*BETADRP_D
            PQQDRPI(L,K)=TMPQQI*(FXDRPR(L ,K  )+FXDRPR(L ,K+1)
     &                          +FXDRPR(LE,K  )+FXDRPR(LE,K+1)
     &                          +FYDRPR(L ,K  )+FYDRPR(L ,K+1)
     &                          +FYDRPR(LN,K  )+FYDRPR(LN,K+1))
            PQQDRPE(L,K)=TMPQQE*(FXDRPR(L ,K  )*U(L ,K  )*U(L ,K  )
     &                          +FXDRPR(L ,K+1)*U(L ,K+1)*U(L ,K+1)
     &                          +FXDRPR(LE,K  )*U(LE,K  )*U(LE,K  )
     &                          +FXDRPR(LE,K+1)*U(LE,K+1)*U(LE,K+1)
     &                          +FYDRPR(L ,K  )*V(L ,K  )*V(L ,K  )
     &                          +FYDRPR(L ,K+1)*V(L ,K+1)*V(L ,K+1)
     &                          +FYDRPR(LN,K  )*V(LN,K  )*V(LN,K  )
     &                          +FYDRPR(LN,K+1)*V(LN,K+1)*V(LN,K+1))
          ENDDO
        ENDDO
      ENDIF

      IF(ISWAVE.LE.1.OR.ISWAVE.EQ.3)THEN
        IF(IDRYTBP.EQ.0)THEN
!$OMP PARALLEL PRIVATE(LN,LS,LE,
!$OMP& PQQB,PQQU,PQQ,PQQL)
          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              LN=LNC(L)
              LS=LSC(L)
              LE=LEAST(L)
              UUU(L,K)=QQ(L,K)*H1P(L)
     &            +DELT*(FUHU(L,K)-FUHU(LE,K)+FVHU(L,K)-FVHU(LN,K)
     &            +(FWQQ(L,K)-FWQQ(L,K+1))*DZIG(K))*DXYIP(L)

              VVV(L,K)=QQL(L,K)*H1P(L)*H1P(L)
     &            +DELT*(FUHV(L,K)-FUHV(LE,K)+FVHV(L,K)-FVHV(LN,K)
     &            +(FWQQL(L,K)-FWQQL(L,K+1))*DZIG(K))*DXYIP(L)
            END DO
!$OMP END DO NOWAIT
          END DO

          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              UUU(L,K)=MAX(UUU(L,K),0.)
              VVV(L,K)=MAX(VVV(L,K),0.)
              LN=LNC(L)
              LE=LEAST(L)
              PQQB=AB(L,K)*GP*HP(L)*DZIG(K)*(B(L,K+1)-B(L,K))
              PQQU=AV(L,K)*DZIGSD4(K)*(U(LE,K+1)-U(LE,K)+U(L,K+1)-
     &            U(L,K))**2+AV(L,K)*DZIGSD4(K)*(V(LN,K+1)-V(LN,K)+V(
     &            L,K+1)-V(L,K))**2

              PQQ=DELT*(PQQB+PQQU+
     &                  PQQVEGE(L,K)+PQQMHKE(L,K)+PQQSUPE(L,K)+PQQDRPE(L,K))
              UUU(L,K)=UUU(L,K)+2.*PQQ
              PQQL=DELT*H1P(L)*(CTE3*PQQB+CTE1*PQQU
     &                         +CE4VEG*PQQVEGE(L,K)
     &                         +CE4MHK*PQQMHKE(L,K)
     &                         +CE4SUP*PQQSUPE(L,K)
     &                         +CE4DRP*PQQDRPE(L,K))
              VVV(L,K)=VVV(L,K)+DML(L,K)*PQQL
            ENDDO
!$OMP END DO NOWAIT
          ENDDO
!$OMP END PARALLEL
        ELSE
!$OMP PARALLEL PRIVATE(LN,LS,LE,
!$OMP& PQQB,PQQU,PQQ,PQQL)
          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                LN=LNC(L)
                LE=LEAST(L)
                UUU(L,K)=QQ(L,K)*H1P(L)
     &              +DELT*(FUHU(L,K)-FUHU(LE,K)+FVHU(L,K)-FVHU(LN,K)
     &              +(FWQQ(L,K)-FWQQ(L,K+1))*DZIG(K))*DXYIP(L)
                VVV(L,K)=QQL(L,K)*H1P(L)*H1P(L)
     &              +DELT*(FUHV(L,K)-FUHV(LE,K)+FVHV(L,K)-FVHV(LN,K)
     &              +(FWQQL(L,K)-FWQQL(L,K+1))*DZIG(K))*DXYIP(L)

                UUU(L,K)=MAX(UUU(L,K),0.)
                VVV(L,K)=MAX(VVV(L,K),0.)
              ELSE
                UUU(L,K)=0.0
                VVV(L,K)=0.0
              ENDIF
            ENDDO
!$OMP END DO NOWAIT
          ENDDO

          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                LN=LNC(L)
                LS=LSC(L)
                LE=LEAST(L)
                PQQB=AB(L,K)*GP*HP(L)*DZIG(K)*(B(L,K+1)-B(L,K))
                PQQU=AV(L,K)*DZIGSD4(K)*(U(LE,K+1)-U(LE,K)+U(L,K+1)
     &              -U(L,K))**2+AV(L,K)*DZIGSD4(K)*(V(LN,K+1)-V(LN,K)
     &              +V(L,K+1)-V(L,K))**2
              PQQ=DELT*(PQQB+PQQU+PQQVEGE(L,K)
     &                 +PQQMHKE(L,K)+PQQSUPE(L,K)+PQQDRPE(L,K))
              UUU(L,K)=UUU(L,K)+2.*PQQ
              PQQL=DELT*H1P(L)*(CTE3*PQQB+CTE1*PQQU
     &                         +CE4VEG*PQQVEGE(L,K)
     &                         +CE4MHK*PQQMHKE(L,K)
     &                         +CE4SUP*PQQSUPE(L,K)
     &                         +CE4DRP*PQQDRPE(L,K))
                VVV(L,K)=VVV(L,K)+DML(L,K)*PQQL
              ENDIF
            ENDDO
!$OMP END DO NOWAIT
          ENDDO
!$OMP END PARALLEL
        ENDIF
      ENDIF
C
C *** DSLLC BEGIN BLOCK
C
      IF(ISWAVE.EQ.2)THEN
        IF(N.LT.NTSWV)THEN
          TMPVAL=FLOAT(N)/FLOAT(NTSWV)
          WVFACT=0.5-0.5*COS(PI*TMPVAL)
        ELSE
          WVFACT=1.0
        ENDIF
        DO K=1,KS
          DO L=2,LA
            TVAR1W(L,K)=WVDTKEM(K)*WVDISP(L,K)+WVDTKEP(K)*WVDISP(L,K+1)
          ENDDO
        ENDDO
        DO K=1,KS
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              LN=LNC(L)
              LS=LSC(L)
              LE=LEAST(L)
              PQQB=AB(L,K)*GP*HP(L)*DZIG(K)*(B(L,K+1)-B(L,K))
              PQQU=AV(L,K)*DZIGSD4(K)*
     &          (U(LE,K+1)-U(LE,K)+U(L,K+1)-U(L,K))**2
              PQQV=AV(L,K)*DZIGSD4(K)*
     &          (V(LN,K+1)-V(LN,K)+V(L,K+1)-V(L,K))**2
              PQQW= WVFACT*TVAR1W(L,K)
              PQQ=DELT*(PQQU+PQQV+PQQB+PQQW
     &                 +PQQVEGE(L,K)+PQQMHKE(L,K)+PQQSUPE(L,K)+PQQDRPE(L,K))
              FFTMP=MAX(FUHU(L,K)-FUHU(LE,K)+FVHU(L,K)-FVHU(LN,K) +
     &          (FWQQ(L,K)-FWQQ(L,K+1))*DZIG(K),0.)
              UUU(L,K)=QQ(L,K)*H1P(L)+DELT*FFTMP*DXYIP(L) + 2.*PQQ
              PQQL=DELT*H1P(L)*(CTE3*PQQB+CTE1*(PQQU+PQQV+PQQW)
     &                         +CE4VEG*PQQVEGE(L,K)
     &                         +CE4MHK*PQQMHKE(L,K)
     &                         +CE4SUP*PQQSUPE(L,K)
     &                         +CE4DRP*PQQDRPE(L,K))
              FFTMP=MAX(FUHV(L,K)-FUHV(LE,K)+FVHV(L,K)-FVHV(LN,K) +
     &          (FWQQL(L,K)-FWQQL(L,K+1))*DZIG(K),0.)
              VVV(L,K)=QQL(L,K)*H1P(L)*H1P(L)+DELT*FFTMP*DXYIP(L) +
     &          DML(L,K)*PQQL
            ENDIF
          ENDDO
        ENDDO
      ENDIF
C
C *** DSLLC END BLOCK
C
      IF(KC.LE.2)THEN
        IF(IDRYTBP.EQ.0)THEN
          DO L=2,LA
            CLQTMP=-DELT*CDZKK(1)*AQ(L,1)*HPI(L)
            CUQTMP=-DELT*CDZKKP(1)*AQ(L,2)*HPI(L)
            CMQTMP=1.-CLQTMP-CUQTMP
     &         +2.*DELT*QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L))
     &         +0.5*DELT*HPI(L)*(PQQVEGI(L,1)+PQQMHKI(L,1)+PQQSUPI(L,1)+PQQDRPI(L,1))
            CMQLTMP=1.-CLQTMP-CUQTMP
     &          +DELT*(QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L)))*(1.
     &          +CTE2*DML(L,1)*DML(L,1)*FPROX(1))
     &          +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,1)
     &                       +CE4MHK*PQQMHKI(L,1)
     &                       +CE4SUP*PQQSUPI(L,1)
     &                       +CE4DRP*PQQDRPI(L,1))

            EQ=1./CMQTMP
            EQL=1./CMQLTMP
            CU1(L,1)=CUQTMP*EQ
            CU2(L,1)=CUQTMP*EQL
            UUU(L,1)=(UUU(L,1)-CLQTMP*HP(L)*QQ(L,0)
     &          -CUQTMP*HP(L)*QQ(L,KC))*EQ
            VVV(L,1)=VVV(L,1)*EQL
          ENDDO
        ELSE
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              CLQTMP=-DELT*CDZKK(1)*AQ(L,1)*HPI(L)
              CUQTMP=-DELT*CDZKKP(1)*AQ(L,2)*HPI(L)
              CMQTMP=1.-CLQTMP-CUQTMP
     &            +2.*DELT*QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L))
     &            +DELT*HPI(L)*(PQQVEGI(L,1)+PQQMHKI(L,1)+PQQSUPI(L,1)+PQQDRPI(L,1))
              CMQLTMP=1.-CLQTMP-CUQTMP
     &            +DELT*(QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L)))*(1.
     &            +CTE2*DML(L,1)*DML(L,1)*FPROX(1))
     &            +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,1)
     &                         +CE4MHK*PQQMHKI(L,1)
     &                         +CE4SUP*PQQSUPI(L,1)
     &                         +CE4DRP*PQQDRPI(L,1))
              EQ=1./CMQTMP
              EQL=1./CMQLTMP
              CU1(L,1)=CUQTMP*EQ
              CU2(L,1)=CUQTMP*EQL
              UUU(L,1)=(UUU(L,1)-CLQTMP*HP(L)*QQ(L,0)
     &            -CUQTMP*HP(L)*QQ(L,KC))*EQ
              VVV(L,1)=VVV(L,1)*EQL
            ENDIF
          ENDDO
        ENDIF
      ENDIF
      IF(KC.GT.2)THEN
        IF(IDRYTBP.EQ.0)THEN
!$OMP PARALLEL PRIVATE(CLQTMP,CUQTMP,CMQTMP,
!$OMP& CMQLTMP,EQ,EQL)
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA
            CLQTMP=-DELT*CDZKK(1)*AQ(L,1)*HPI(L)
            CUQTMP=-DELT*CDZKKP(1)*AQ(L,2)*HPI(L)
            CMQTMP=1.-CLQTMP-CUQTMP
     &          +2.*DELT*QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L))
     &          +DELT*HPI(L)*(PQQVEGI(L,1)+PQQMHKI(L,1)+PQQSUPI(L,1)+PQQDRPI(L,1))
            CMQLTMP=1.-CLQTMP-CUQTMP
     &          +DELT*(QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L)))*(1.
     &          +CTE2*DML(L,1)*DML(L,1)*FPROX(1))
     &          +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,1)
     &                       +CE4MHK*PQQMHKI(L,1)
     &                       +CE4SUP*PQQSUPI(L,1)
     &                       +CE4DRP*PQQDRPI(L,1))
            EQ=1./CMQTMP
            EQL=1./CMQLTMP
            CU1(L,1)=CUQTMP*EQ
            CU2(L,1)=CUQTMP*EQL
            UUU(L,1)=(UUU(L,1)-CLQTMP*HP(L)*QQ(L,0))*EQ
            VVV(L,1)=VVV(L,1)*EQL
            CUQTMP=-DELT*CDZKKP(KS)*AQ(L,KC)*HPI(L)
            UUU(L,KS)=UUU(L,KS)-CUQTMP*HP(L)*QQ(L,KC)
          ENDDO
!$OMP END DO NOWAIT
          DO K=2,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              CLQTMP=-DELT*CDZKK(K)*AQ(L,K)*HPI(L)
              CUQTMP=-DELT*CDZKKP(K)*AQ(L,K+1)*HPI(L)
              CMQTMP=1.-CLQTMP-CUQTMP
     &            +2.*DELT*QQSQR(L,K)/(CTURBB1(L,K)*DML(L,K)*HP(L))
     &            +DELT*HPI(L)*(PQQVEGI(L,K)+PQQMHKI(L,K)+PQQSUPI(L,K)+PQQDRPI(L,K))
              CMQLTMP=1.-CLQTMP-CUQTMP
     &            +DELT*(QQSQR(L,K)/(CTURBB1(L,K)*DML(L,K)*HP(L)))*(1.
     &            +CTE2*DML(L,K)*DML(L,K)*FPROX(K))
     &            +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,K)
     &                         +CE4MHK*PQQMHKI(L,K)
     &                         +CE4SUP*PQQSUPI(L,K)
     &                         +CE4DRP*PQQDRPI(L,K))
              EQ=1./(CMQTMP-CLQTMP*CU1(L,K-1))
              EQL=1./(CMQLTMP-CLQTMP*CU2(L,K-1))
              CU1(L,K)=CUQTMP*EQ
              CU2(L,K)=CUQTMP*EQL
              UUU(L,K)=(UUU(L,K)-CLQTMP*UUU(L,K-1))*EQ
              VVV(L,K)=(VVV(L,K)-CLQTMP*VVV(L,K-1))*EQL
            ENDDO
!$OMP END DO NOWAIT
          ENDDO

          DO K=KS-1,1,-1
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              UUU(L,K)=UUU(L,K)-CU1(L,K)*UUU(L,K+1)
              VVV(L,K)=VVV(L,K)-CU2(L,K)*VVV(L,K+1)
            ENDDO
!$OMP END DO NOWAIT
          ENDDO
!$OMP END PARALLEL

        ELSE
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              CLQTMP=-DELT*CDZKK(1)*AQ(L,1)*HPI(L)
              CUQTMP=-DELT*CDZKKP(1)*AQ(L,2)*HPI(L)
              CMQTMP=1.-CLQTMP-CUQTMP
     &            +2.*DELT*QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L))
     &            +DELT*HPI(L)*(PQQVEGI(L,1)+PQQMHKI(L,1)+PQQSUPI(L,1)+PQQDRPI(L,1))
              CMQLTMP=1.-CLQTMP-CUQTMP
     &            +DELT*(QQSQR(L,1)/(CTURBB1(L,1)*DML(L,1)*HP(L)))*(1.
     &            +CTE2*DML(L,1)*DML(L,1)*FPROX(1))
     &            +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,1)
     &                         +CE4MHK*PQQMHKI(L,1)
     &                         +CE4SUP*PQQSUPI(L,1)
     &                         +CE4DRP*PQQDRPI(L,1))
              EQ=1./CMQTMP
              EQL=1./CMQLTMP
              CU1(L,1)=CUQTMP*EQ
              CU2(L,1)=CUQTMP*EQL
              UUU(L,1)=(UUU(L,1)-CLQTMP*HP(L)*QQ(L,0))*EQ
              VVV(L,1)=VVV(L,1)*EQL
              CUQTMP=-DELT*CDZKKP(KS)*AQ(L,KC)*HPI(L)
              UUU(L,KS)=UUU(L,KS)-CUQTMP*HP(L)*QQ(L,KC)
            ENDIF
          ENDDO
          DO K=2,KS
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                CLQTMP=-DELT*CDZKK(K)*AQ(L,K)*HPI(L)
                CUQTMP=-DELT*CDZKKP(K)*AQ(L,K+1)*HPI(L)
                CMQTMP=1.-CLQTMP-CUQTMP
     &             +2.*DELT*QQSQR(L,K)/(CTURBB1(L,K)*DML(L,K)*HP(L))
     &             +DELT*HPI(L)*(PQQVEGI(L,K)+PQQMHKI(L,K)+PQQSUPI(L,K)+PQQDRPI(L,K))
               CMQLTMP=1.-CLQTMP-CUQTMP
     &              +DELT*(QQSQR(L,K)/(CTURBB1(L,K)*DML(L,K)*HP(L)))
     &              *(1.+CTE2*DML(L,K)*DML(L,K)*FPROX(K))
     &              +DELT*HPI(L)*(CE4VEG*PQQVEGI(L,K)
     &                           +CE4MHK*PQQMHKI(L,K)
     &                           +CE4SUP*PQQSUPI(L,K)
     &                           +CE4DRP*PQQDRPI(L,K))
                EQ=1./(CMQTMP-CLQTMP*CU1(L,K-1))
                EQL=1./(CMQLTMP-CLQTMP*CU2(L,K-1))
                CU1(L,K)=CUQTMP*EQ
                CU2(L,K)=CUQTMP*EQL
                UUU(L,K)=(UUU(L,K)-CLQTMP*UUU(L,K-1))*EQ
                VVV(L,K)=(VVV(L,K)-CLQTMP*VVV(L,K-1))*EQL
              ENDIF
            ENDDO
          ENDDO
          DO K=KS-1,1,-1
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                UUU(L,K)=UUU(L,K)-CU1(L,K)*UUU(L,K+1)
                VVV(L,K)=VVV(L,K)-CU2(L,K)*VVV(L,K+1)
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ENDIF
      IF(IDRYTBP.EQ.0)THEN
C
C **  ORIGINAL FORM MODIFED FOR DIMENSIONAL LENGTH SCALE TRANSPORT
C
        IF(ISTOPT(0).EQ.1)THEN

!$OMP PARALLEL PRIVATE(QQHDH,DMLTMP,DELB,DMLMAX)
          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              QQ1(L,K)=QQ(L,K)
              QQHDH=UUU(L,K)*HPI(L)
              QQ(L,K)=MAX(QQHDH,QQMIN)
              QQL1(L,K)=QQL(L,K)
              QQHDH=VVV(L,K)*HPI(L)*HPI(L)
              QQL(L,K)=MAX(QQHDH,QQLMIN)
              DMLTMP=QQL(L,K)/QQ(L,K)
              DMLTMP=MAX(DMLTMP,DMLMIN)
              DELB=B(L,K)-B(L,K+1)
              IF(DELB.GT.0.0)THEN
                DMLMAX=0.53*SQRT(QQ(L,K)/(G*HP(L)*DZIG(K)*DELB))
                DML(L,K)=MIN(DMLMAX,DMLTMP)
              ELSE
                DML(L,K)=DMLTMP
              ENDIF
            ENDDO
          ENDDO
!$OMP END PARALLEL
        ENDIF
C
C **  BUCHARD'S MODIFED CLOSURE FOR DIMENSIONAL LENGTH SCALE TRANSPORT
C
        IF(ISTOPT(0).EQ.2)THEN
!$OMP PARALLEL PRIVATE(QQHDH,DMLTMP)
          DO K=1,KS
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              QQ1(L,K)=QQ(L,K)
              QQHDH=UUU(L,K)*HPI(L)
              QQ(L,K)=MAX(QQHDH,QQMIN)
              QQL1(L,K)=QQL(L,K)
              QQHDH=VVV(L,K)*HPI(L)*HPI(L)
              QQL(L,K)=MAX(QQHDH,QQLMIN)
              DMLTMP=QQL(L,K)/QQ(L,K)
              DML(L,K)=MAX(DMLTMP,DMLMIN)
            ENDDO
          ENDDO
!$OMP END PARALLEL
        ENDIF
      ELSE     ! else on if LMASKDRY(L)
        DO K=1,KS
          DO L=2,LA
            IF(LMASKDRY(L))THEN
              QQ1(L,K)=QQ(L,K)
              QQHDH=UUU(L,K)*HPI(L)
              QQ(L,K)=MAX(QQHDH,QQMIN)
            ENDIF
          ENDDO
        ENDDO
C
C **  ORIGINAL FORM MODIFED FOR DIMENSIONAL LENGTH SCALE TRANSPORT
C
        IF(ISTOPT(0).EQ.1)THEN
          DO K=1,KS
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                QQL1(L,K)=QQL(L,K)
                QQHDH=VVV(L,K)*HPI(L)*HPI(L)
                QQL(L,K)=MAX(QQHDH,QQLMIN)
                DMLTMP=QQL(L,K)/QQ(L,K)
                DMLTMP=MAX(DMLTMP,DMLMIN)
                DELB=B(L,K)-B(L,K+1)
                IF(DELB.GT.0.0)THEN
                  DMLMAX=0.53*SQRT(QQ(L,K)/(G*HP(L)*DZIG(K)*DELB))
                  DML(L,K)=MIN(DMLMAX,DMLTMP)
                ELSE
                  DML(L,K)=DMLTMP
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
C
C **  BUCHARD'S MODIFED CLOSURE FOR DIMENSIONAL LENGTH SCALE TRANSPORT
C
        IF(ISTOPT(0).EQ.2)THEN
          DO K=1,KS
            DO L=2,LA
              IF(LMASKDRY(L))THEN
                QQL1(L,K)=QQL(L,K)
                QQHDH=VVV(L,K)*HPI(L)*HPI(L)
                QQL(L,K)=MAX(QQHDH,QQLMIN)
                DMLTMP=QQL(L,K)/QQ(L,K)
                DML(L,K)=MAX(DMLTMP,DMLMIN)
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ENDIF
C
C      QQMXSV=-1.E+12
C      QQMNSV=1.E+12
C      QQLMXSV=-1.E+12
C      QQLMNSV=1.E+12
C
      DO K=1,KS
        DO LL=1,NCBS
          L=LCBS(LL)
          LN=LNC(L)
          QQ(L,K)=QQ(LN,K)
          QQL(L,K)=QQL(LN,K)
          DML(L,K)=DML(LN,K)
        ENDDO
      ENDDO
      DO K=1,KS
        DO LL=1,NCBW
          L=LCBW(LL)
          LE=LEAST(L)
          QQ(L,K)=QQ(LE,K)
          QQL(L,K)=QQL(LE,K)
          DML(L,K)=DML(LE,K)
        ENDDO
      ENDDO
      DO K=1,KS
        DO LL=1,NCBE
          L=LCBE(LL)
          LW=LWEST(L)
          QQ(L,K)=QQ(LW,K)
          QQL(L,K)=QQL(LW,K)
          DML(L,K)=DML(LW,K)
        ENDDO
      ENDDO
      DO K=1,KS
        DO LL=1,NCBN
          L=LCBN(LL)
          LS=LSC(L)
          QQ(L,K)=QQ(LS,K)
          QQL(L,K)=QQL(LS,K)
          DML(L,K)=DML(LS,K)
        ENDDO
      ENDDO
C *** DSLLC BEGIN BLOCK
      DO K=1,KS
        DO L=1,LC
          QQSQR(L,K)=SQRT(QQ(L,K))  ! *** DSLLC
        ENDDO
      ENDDO
C *** DSLLC END BLOCK
  110 FORMAT('    I    J   QQ BOT        QQ MID        QQ SURF',
     &    '       PROD+ADV      1./DIAGON')
  111 FORMAT(2I5,5E14.5)
  600 FORMAT('N,QX,QN,QLX,QLN,CX,CN=',I5,6E12.4)
  601 FORMAT('NEG QQ I,J,K,QQ=',3I5,E13.5)
      RETURN
      END


program generate_testfile
    implicit none
    integer :: unit
    real(4) :: x4(2)
    real(8) :: x8(2)
    integer(4) :: i4(2), ishort
    integer(8) :: i8(2), ilong

    open(file="testfile.bin", newunit=unit, form="unformatted")

    x4 = (/1e0, 2e0/)
    x8 = (/3d0, 4d0/)
    i4 = (/5, 6/)
    i8 = (/7, 8/)
    ishort = 1234
    ilong = 123456789
    ! Write scalars
    write(unit) ishort
    write(unit) ilong
    write(unit) 1e0
    write(unit) 2d0

    ! Write vectors
    write(unit) i4
    write(unit) i8
    write(unit) x4
    write(unit) x8
    close(unit)

    open(file="testfile_attrs.bin", newunit=unit, form="unformatted")
    ! Write scalars
    write(unit) 1
    write(unit) 1, 2, 3, 4
    write(unit) 1d0, 2d0, 3d0, 4d0
    write(unit) "test"
    close(unit)
end program


text:
!=============================================================================!
        subroutine rpenta1p( nsys, np, nblk, mat, rhs, per, per2, code )
!=============================================================================!
!
!  Solves a block pentadiagonal system of equations without pivoting.
!  Vectorized over the first index.
!
!  This routine solves a PERIODIC block pentadiagonal system.
!
!  WARNING: Since pivoting is not done, this algorithm could fail
!
!  input:
!    nsys       : number of systems to vectorize over
!    n          : number of rows in the systems
!    nblk       : dimension of the blocks in each row
!    code       : 0 = LU factorization only
!               : 1 = Forward and back solves given LU
!               : 2 = LU and solve
!
!  inout:
!    mat        : the block pentadiagonal matrix to be solved
!    rhs        : the right hand side
!=============================================================================!
        implicit none

        integer np, nblk, nsys, code
        real mat(5,nblk,nblk,nsys,np), rhs(nblk,nsys,np)
        real per(2,nblk,nblk,nsys,np), per2(np-2:np-1,nblk,nblk,nsys,np)
        real mult, fact

!.... Useful constants

        real zero, one
        parameter (zero = 0.0, one = 1.0)

!.... local variables

        integer nc, i, j, l, m, n, p, q, r, iv
!=============================================================================!

!.... don't solve for the redundant nodes

        n = np - 1

!=============================================================================!
!                       L U   F A C T O R I Z A T I O N
!=============================================================================!
        if (code .eq. 0 .or. code .eq. 2) then

!$omp parallel do private(nc,i,l,m,p,q,r,iv,fact,mult)
        do iv = 1, nsys

        do i = 1, np
           do q = 1, nblk
              do p = 1, nblk
                 per(1,p,q,iv,i) = zero
                 per(2,p,q,iv,i) = zero
                 per2(n-1,p,q,iv,i) = zero
                 per2(n,p,q,iv,i) = zero
              end do
           end do
        end do

!.... setup the periodic fillin vector

        do q = 1, nblk
          do p = 1, nblk
            per(1,p,q,iv,1) = mat(1,p,q,iv,1)
            per(2,p,q,iv,1) = mat(2,p,q,iv,1)
            per(1,p,q,iv,2) = zero
            per(2,p,q,iv,2) = mat(1,p,q,iv,2)

            per2(n-1,p,q,iv,1) = mat(5,p,q,iv,n-1)
            per2(n-1,p,q,iv,2) = zero
            per2(n,p,q,iv,1)   = mat(4,p,q,iv,n)
            per2(n,p,q,iv,2)   = mat(5,p,q,iv,n)
          end do
        end do

        do i = 1, n - 6

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
                mat(5,l,q,iv,i) = mat(5,l,q,iv,i) + mult * mat(5,nc,q,iv,i)
                per(1,l,q,iv,i) = per(1,l,q,iv,i) + mult * per(1,nc,q,iv,i)
                per(2,l,q,iv,i) = per(2,l,q,iv,i) + mult * per(2,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          nc = 1

          p = i + 1
          r = i + 2

          fact = one / mat(3,nc,nc,iv,i)
          do m = 1, nblk
            mult = -mat(2,m,nc,iv,p) * fact
            mat(2,m,nc,iv,p) = mult
            do q = nc + 1, nblk
              mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
              per(1,m,q,iv,p) = per(1,m,q,iv,p) + mult * per(1,nc,q,iv,i)
              per(2,m,q,iv,p) = per(2,m,q,iv,p) + mult * per(2,nc,q,iv,i)
            end do

            mult = -mat(1,m,nc,iv,r) * fact
            mat(1,m,nc,iv,r) = mult
            do q = nc + 1, nblk
              mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
              mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
              per(1,m,q,iv,r) = mult * per(1,nc,q,iv,i)
              per(2,m,q,iv,r) = mult * per(2,nc,q,iv,i)
            end do

!.... the periodic stuff on the next-to-the-last row

            mult = -per2(n-1,m,nc,iv,i) * fact
            per2(n-1,m,nc,iv,i) = mult
            do q = nc + 1, nblk
              per2(n-1,m,q,iv,i) = per2(n-1,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              per2(n-1,m,q,iv,p) = per2(n-1,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              per2(n-1,m,q,iv,r) = mult * mat(5,nc,q,iv,i)
              mat(3,m,q,iv,n-1) = mat(3,m,q,iv,n-1) + mult * per(1,nc,q,iv,i)
              mat(4,m,q,iv,n-1) = mat(4,m,q,iv,n-1) + mult * per(2,nc,q,iv,i)
            end do

!.... the periodic stuff on the last row

            mult = -per2(n,m,nc,iv,i) * fact
            per2(n,m,nc,iv,i) = mult
            do q = nc + 1, nblk
              per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              per2(n,m,q,iv,p) = per2(n,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              per2(n,m,q,iv,r) = mult * mat(5,nc,q,iv,i)
              mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * per(1,nc,q,iv,i)
              mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
            end do

          end do        ! end of loop on m

          do nc = 2, nblk

            p = i + 1
            r = i + 2

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
                per(1,m,q,iv,p) = per(1,m,q,iv,p) + mult * per(1,nc,q,iv,i)
                per(2,m,q,iv,p) = per(2,m,q,iv,p) + mult * per(2,nc,q,iv,i)
              end do

              mult = -mat(1,m,nc,iv,r) * fact
              mat(1,m,nc,iv,r) = mult
              do q = nc + 1, nblk
                mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
                mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                per(1,m,q,iv,r) = per(1,m,q,iv,r) + mult * per(1,nc,q,iv,i)
                per(2,m,q,iv,r) = per(2,m,q,iv,r) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the next-to-the-last row

              mult = -per2(n-1,m,nc,iv,i) * fact
              per2(n-1,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n-1,m,q,iv,i) = per2(n-1,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(n-1,m,q,iv,p) = per2(n-1,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                per2(n-1,m,q,iv,r) = per2(n-1,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(3,m,q,iv,n-1) = mat(3,m,q,iv,n-1) + mult * per(1,nc,q,iv,i)
                mat(4,m,q,iv,n-1) = mat(4,m,q,iv,n-1) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the last row

              mult = -per2(n,m,nc,iv,i) * fact
              per2(n,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(n,m,q,iv,p) = per2(n,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                per2(n,m,q,iv,r) = per2(n,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * per(1,nc,q,iv,i)
                mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
              end do

            end do      ! end of loop on m

          end do        ! end of loop on nc

        end do          ! end of loop on i

!====================================================================================
        i = n - 5

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
                mat(5,l,q,iv,i) = mat(5,l,q,iv,i) + mult * mat(5,nc,q,iv,i)
                per(1,l,q,iv,i) = per(1,l,q,iv,i) + mult * per(1,nc,q,iv,i)
                per(2,l,q,iv,i) = per(2,l,q,iv,i) + mult * per(2,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          nc = 1

          p = i + 1
          r = i + 2

          fact = one / mat(3,nc,nc,iv,i)
          do m = 1, nblk
            mult = -mat(2,m,nc,iv,p) * fact
            mat(2,m,nc,iv,p) = mult
            do q = nc + 1, nblk
              mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
              per(1,m,q,iv,p) = per(1,m,q,iv,p) + mult * per(1,nc,q,iv,i)
              per(2,m,q,iv,p) = per(2,m,q,iv,p) + mult * per(2,nc,q,iv,i)
            end do

            mult = -mat(1,m,nc,iv,r) * fact
            mat(1,m,nc,iv,r) = mult
            do q = nc + 1, nblk
              mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
              mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
              mat(4,m,q,iv,r) = mat(4,m,q,iv,r) + mult * per(1,nc,q,iv,i)
              per(2,m,q,iv,r) = mult * per(2,nc,q,iv,i)
            end do

!.... the periodic stuff on the next-to-the-last row

            mult = -per2(n-1,m,nc,iv,i) * fact
            per2(n-1,m,nc,iv,i) = mult
            do q = nc + 1, nblk
              per2(n-1,m,q,iv,i) = per2(n-1,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              per2(n-1,m,q,iv,p) = per2(n-1,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              mat(1,m,q,iv,n-1) = mat(1,m,q,iv,n-1) + mult * mat(5,nc,q,iv,i)
              mat(3,m,q,iv,n-1) = mat(3,m,q,iv,n-1) + mult * per(1,nc,q,iv,i)
              mat(4,m,q,iv,n-1) = mat(4,m,q,iv,n-1) + mult * per(2,nc,q,iv,i)
            end do

!.... the periodic stuff on the last row

            mult = -per2(n,m,nc,iv,i) * fact
            per2(n,m,nc,iv,i) = mult
            do q = nc + 1, nblk
              per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
            end do
            do q = 1, nblk
              per2(n,m,q,iv,p) = per2(n,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              per2(n,m,q,iv,r) = mult * mat(5,nc,q,iv,i)
              mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * per(1,nc,q,iv,i)
              mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
            end do

          end do        ! end of loop on m

          do nc = 2, nblk

            p = i + 1
            r = i + 2

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
                per(1,m,q,iv,p) = per(1,m,q,iv,p) + mult * per(1,nc,q,iv,i)
                per(2,m,q,iv,p) = per(2,m,q,iv,p) + mult * per(2,nc,q,iv,i)
              end do

              mult = -mat(1,m,nc,iv,r) * fact
              mat(1,m,nc,iv,r) = mult
              do q = nc + 1, nblk
                mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
                mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(4,m,q,iv,r) = mat(4,m,q,iv,r) + mult * per(1,nc,q,iv,i)
                per(2,m,q,iv,r) = per(2,m,q,iv,r) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the next-to-the-last row

              mult = -per2(n-1,m,nc,iv,i) * fact
              per2(n-1,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n-1,m,q,iv,i) = per2(n-1,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(n-1,m,q,iv,p) = per2(n-1,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(1,m,q,iv,n-1) = mat(1,m,q,iv,n-1) + mult * mat(5,nc,q,iv,i)
                mat(3,m,q,iv,n-1) = mat(3,m,q,iv,n-1) + mult * per(1,nc,q,iv,i)
                mat(4,m,q,iv,n-1) = mat(4,m,q,iv,n-1) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the last row

              mult = -per2(n,m,nc,iv,i) * fact
              per2(n,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(n,m,q,iv,p) = per2(n,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                per2(n,m,q,iv,r) = per2(n,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * per(1,nc,q,iv,i)
                mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
              end do

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 4

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
                mat(5,l,q,iv,i) = mat(5,l,q,iv,i) + mult * mat(5,nc,q,iv,i)
                per(1,l,q,iv,i) = per(1,l,q,iv,i) + mult * per(1,nc,q,iv,i)
                per(2,l,q,iv,i) = per(2,l,q,iv,i) + mult * per(2,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
                mat(5,m,q,iv,p) = mat(5,m,q,iv,p) + mult * per(1,nc,q,iv,i)
                per(2,m,q,iv,p) = per(2,m,q,iv,p) + mult * per(2,nc,q,iv,i)
              end do

              mult = -mat(1,m,nc,iv,r) * fact
              mat(1,m,nc,iv,r) = mult
              do q = nc + 1, nblk
                mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
                mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(4,m,q,iv,r) = mat(4,m,q,iv,r) + mult * per(1,nc,q,iv,i)
                mat(5,m,q,iv,r) = mat(5,m,q,iv,r) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the next-to-the-last row

              mult = -per2(n-1,m,nc,iv,i) * fact
              per2(n-1,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n-1,m,q,iv,i) = per2(n-1,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(1,m,q,iv,n-1) = mat(1,m,q,iv,n-1) + mult * mat(4,nc,q,iv,i)
                mat(2,m,q,iv,n-1) = mat(2,m,q,iv,n-1) + mult * mat(5,nc,q,iv,i)
                mat(3,m,q,iv,n-1) = mat(3,m,q,iv,n-1) + mult * per(1,nc,q,iv,i)
                mat(4,m,q,iv,n-1) = mat(4,m,q,iv,n-1) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the last row

              mult = -per2(n,m,nc,iv,i) * fact
              per2(n,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(n,m,q,iv,p) = per2(n,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(1,m,q,iv,n) = mat(1,m,q,iv,n) + mult * mat(5,nc,q,iv,i)
                mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * per(1,nc,q,iv,i)
                mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
              end do

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 3

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
                mat(5,l,q,iv,i) = mat(5,l,q,iv,i) + mult * mat(5,nc,q,iv,i)
                per(2,l,q,iv,i) = per(2,l,q,iv,i) + mult * per(2,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
                mat(5,m,q,iv,p) = mat(5,m,q,iv,p) + mult * per(2,nc,q,iv,i)
              end do

              mult = -mat(1,m,nc,iv,r) * fact
              mat(1,m,nc,iv,r) = mult
              do q = nc + 1, nblk
                mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
                mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
                mat(4,m,q,iv,r) = mat(4,m,q,iv,r) + mult * per(2,nc,q,iv,i)
              end do

!.... the periodic stuff on the last row

              mult = -per2(n,m,nc,iv,i) * fact
              per2(n,m,nc,iv,i) = mult
              do q = nc + 1, nblk
                per2(n,m,q,iv,i) = per2(n,m,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(1,m,q,iv,n) = mat(1,m,q,iv,n) + mult * mat(4,nc,q,iv,i)
                mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult * mat(5,nc,q,iv,i)
                mat(3,m,q,iv,n) = mat(3,m,q,iv,n) + mult * per(2,nc,q,iv,i)
              end do

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 2

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
                mat(5,l,q,iv,i) = mat(5,l,q,iv,i) + mult * mat(5,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
                mat(4,m,q,iv,p) = mat(4,m,q,iv,p) + mult * mat(5,nc,q,iv,i)
              end do

              mult = -mat(1,m,nc,iv,r) * fact
              mat(1,m,nc,iv,r) = mult
              do q = nc + 1, nblk
                mat(1,m,q,iv,r) = mat(1,m,q,iv,r) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,r) = mat(2,m,q,iv,r) + mult * mat(4,nc,q,iv,i)
                mat(3,m,q,iv,r) = mat(3,m,q,iv,r) + mult * mat(5,nc,q,iv,i)
              end do

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 1

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
              mat(3,l,nc,iv,i) = mult
              do q = nc + 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(4,l,q,iv,i) = mat(4,l,q,iv,i) + mult * mat(4,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1

            fact = one / mat(3,nc,nc,iv,i)
            do m = 1, nblk
              mult = -mat(2,m,nc,iv,p) * fact
              mat(2,m,nc,iv,p) = mult
              do q = nc + 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult * mat(3,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult * mat(4,nc,q,iv,i)
              end do
            end do

          end do        ! end of loop on nc

!.... do the last row

        i = n

        do nc = 1, nblk - 1

          do m = nc, nblk - 1
            l = m + 1
            mult = -mat(3,l,nc,iv,i) / mat(3,nc,nc,iv,i)
            mat(3,l,nc,iv,i) = mult
            do q = nc + 1, nblk
              mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult * mat(3,nc,q,iv,i)
            end do
          end do

        end do  ! end of loop on nc

        end do  ! loop over nsys

        end if
!=============================================================================!
!                     F O R W A R D   E L I M I N A T I O N
!=============================================================================!
        if (code.eq.1 .or. code.eq.2) then

!$omp parallel do private(nc,i,l,m,p,q,r,iv,fact,mult)
        do iv = 1, nsys

        do i = 1, n - 6

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          nc = 1

          p = i + 1
          r = i + 2

          do m = 1, nblk
            mult = mat(2,m,nc,iv,p)
            rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

            mult = mat(1,m,nc,iv,r)
            rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the next-to-the-last row

            mult = per2(n-1,m,nc,iv,i)
            rhs(m,iv,n-1) = rhs(m,iv,n-1) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

            mult = per2(n,m,nc,iv,i)
            rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)

          end do        ! end of loop on m

          do nc = 2, nblk

            p = i + 1
            r = i + 2

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

              mult = mat(1,m,nc,iv,r)
              rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the next-to-the-last row

              mult = per2(n-1,m,nc,iv,i)
              rhs(m,iv,n-1) = rhs(m,iv,n-1) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

              mult = per2(n,m,nc,iv,i)
              rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)

            end do      ! end of loop on m

          end do        ! end of loop on nc

        end do          ! end of loop on i

!====================================================================================
        i = n - 5

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          nc = 1

          p = i + 1
          r = i + 2

          do m = 1, nblk
            mult = mat(2,m,nc,iv,p)
            rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

            mult = mat(1,m,nc,iv,r)
            rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the next-to-the-last row

            mult = per2(n-1,m,nc,iv,i)
            rhs(m,iv,n-1) = rhs(m,iv,n-1) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

            mult = per2(n,m,nc,iv,i)
            rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)
          end do        ! end of loop on m

          do nc = 2, nblk

            p = i + 1
            r = i + 2

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

              mult = mat(1,m,nc,iv,r)
              rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the next-to-the-last row

              mult = per2(n-1,m,nc,iv,i)
              rhs(m,iv,n-1) = rhs(m,iv,n-1) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

              mult = per2(n,m,nc,iv,i)
              rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 4

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

              mult = mat(1,m,nc,iv,r)
              rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the next-to-the-last row

              mult = per2(n-1,m,nc,iv,i)
              rhs(m,iv,n-1) = rhs(m,iv,n-1) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

              mult = per2(n,m,nc,iv,i)
              rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 3

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

              mult = mat(1,m,nc,iv,r)
              rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)

!.... the periodic stuff on the last row

              mult = per2(n,m,nc,iv,i)
              rhs(m,iv,n) = rhs(m,iv,n) + mult * rhs(nc,iv,i)

            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 2

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1
            r = i + 2

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)

              mult = mat(1,m,nc,iv,r)
              rhs(m,iv,r) = rhs(m,iv,r) + mult * rhs(nc,iv,i)
            end do      ! end of loop on m

          end do        ! end of loop on nc

!====================================================================================
        i = n - 1

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult = mat(3,l,nc,iv,i)
              rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p = i + 1

            do m = 1, nblk
              mult = mat(2,m,nc,iv,p)
              rhs(m,iv,p) = rhs(m,iv,p) + mult * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... do the last row

        i = n

        do nc = 1, nblk - 1

          do m = nc, nblk - 1
            l = m + 1
            mult = mat(3,l,nc,iv,i)
            rhs(l,iv,i) = rhs(l,iv,i) + mult * rhs(nc,iv,i)
          end do

        end do  ! end of loop on nc

!=============================================================================!
!                   B A C K W A R D   S U B S T I T U T I O N
!=============================================================================!

!.... do the last row first

        i = n

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(3,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(3,m,m,iv,i)
        end do

!.... do the next-to-the-last row

        i = n - 1
        p = i + 1

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(3,m,q,iv,i)
          end do
          do q = 1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(4,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(3,m,m,iv,i)
        end do

        i = n - 2
        p = i + 1
        r = i + 2

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(3,m,q,iv,i)
          end do
          do q = 1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(4,m,q,iv,i)
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,r) * mat(5,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(3,m,m,iv,i)
        end do

        i = n - 3
        p = i + 1
        r = i + 2

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(3,m,q,iv,i)
          end do
          do q = 1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(4,m,q,iv,i)
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,r) * mat(5,m,q,iv,i)
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,n) * per(2,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(3,m,m,iv,i)
        end do

!.... now do the rest of the rows in reverse order

        do i = n - 4, 1, -1

          p = i + 1
          r = i + 2

          do m = nblk, 1, -1
            do q = m+1, nblk
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(3,m,q,iv,i)
            end do
            do q = 1, nblk
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(4,m,q,iv,i)
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,r) * mat(5,m,q,iv,i)
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,n-1) * per(1,m,q,iv,i)
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,n  ) * per(2,m,q,iv,i)
            end do
            rhs(m,iv,i) = rhs(m,iv,i) / mat(3,m,m,iv,i)
          end do

        end do

!.... account for periodicity

        do q = 1, nblk
          rhs(q,iv,np) = rhs(q,iv,1)
        end do

        end do  ! loop over nsys

        end if
!=============================================================================!
        return
        end