program qcom

!     pgf90 -c print.f
!     pgf90 -o qcom QCOM.f90 print.o
!     ./qcom

!     v. 2 (Sept 2000) --
!     Array subscript ranges now start with 0 instead of 1.	

!     v. 3 (Feb 2009) --
!     All comments now start with !
!     Added f90 declarations and procedures.

      parameter (jt = 20, kt = 20)
      parameter (jv = jt, kv = kt)
      real, dimension (0:jv+1, 0:kv+1) :: v
      real, dimension (1:jv, 1:kv, 2) :: fv

!     Assign parameters for run

!     Initialize all values of arrays for v, w, theta, pi, fv, fw, ftheta, fpi,
!     including values for boundary points for v, w, theta, pi.

      CALL INIT  ! initialize all variables

      ITT = 1 ! itt is time step index

!     USE FORWARD SCHEME TO do first step

      A = 1.
      B = 0.
      N1 = MOD ( ITT    , 2 ) + 1
      N2 = MOD ( ITT - 1, 2 ) + 1

      CALL STEP ( N1, N2, A, B ) ! do first time step

!     ADAMS - BASHFORTH TWO - LEVEL SCHEME

      A =   3. / 2. 
      B = - 1. / 2. 

      ITTNOW = ITT  + 1

      DO ITT = ITTNOW, ITTMAX

      N1 = MOD ( ITT    , 2 ) + 1
      N2 = MOD ( ITT - 1, 2 ) + 1

      CALL STEP ( N1, N2, A, B ) ! do subsequent time steps

      end do
	
!     END-OF-RUN OUTPUT ROUTINES GO HERE

contains

      SUBROUTINE STEP ( N1, N2, A, B )

!     This is the entire subroutine.

      CALL RCALC ( N2 ) ! calculate forcing terms from variables at current time
      CALL AB ( N1, N2, A, B ) ! update variables using a time scheme
      CALL BOUND ! apply boundary conditions to variables
!
      END SUBROUTINE STEP

      SUBROUTINE RCALC (  N2 )

!     CALCULATES FORCING TERMS FOR V(J,K), ETC.; STORES THEM  IN FV(J,K,N2), ETC.

      DO K = 1, KT
      DO J = 1, JT
      FV(J,K,N2)  = [fv for v(j,k)]
      END DO
      END DO

!     ETC (forcing for w, theta, and pi)

      END SUBROUTINE RCALC

      SUBROUTINE AB ( N1, N2, A, B )

!     THE FOLLOWING LOOP UPDATES V USING EITHER THE FORWARD OR THE ADAMS-BASHFORTH 
!     SCHEME DEPENDING ON THE VALUES OF A, B.  
!     SUBSCRIPT N2 OF FV ALWAYS REFERS TO THE MOST RECENTLY CALCULATED VALUES FOR FV.

      DO K = 1, KT
      DO J = 1, JT
      V(J,K)  = V(J,K)  + DT * ( A * FV(J,K,N2)  + B * FV(J,K,N1) )
      END DO
      END DO

!     ETC (update w, theta, and pi)

      END SUBROUTINE AB

      SUBROUTINE BOUND

!     apply boundary conditions for v

      do j = 1, jt
      v(j,0) = v(j,1) ! free slip b.c.
      v(j,kt+1) = v(j,kt) ! free slip b.c.
      end do

      do k = 0, kt+1
      v(0,k) = v(jt,k) ! periodic b.c.
      v(jt+1,k) = v(1,k) ! periodic b.c.
      end do

!     ETC (apply b.c. for w, theta, and pi)

      END SUBROUTINE BOUND
      
      SUBROUTINE INIT
		
!     initialize all variables 
		
      CALL BOUND
      
      END SUBROUTINE INIT
      
end program qcom