! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE MATMULT_SSS_NTERM ( MAT_A_NAME, NROW_A, NTERM_A, SYM_A, I_A, J_A,                                                 &
                                     MAT_B_NAME, NCOL_B, NTERM_B, SYM_B, J_B, I_B, AROW_MAX_TERMS, MAT_C_NAME, NTERM_C )

! Setup routine for subr MATMULT_SSS (which performs the matrix operation C = cons*A*B With A and B in sparse format
! Matrices A and C are stored in compressed row storage (CRS) format. In addition, if A is symmetric it can be stored  with only the
! terms on, and above, the diagonal. Matrix B must be stored in compressed col storage (CCS) format.

! Input  matrix A is stored in compressed row storage (CRS) format using arrays I_A(NROW_A+1), J_A(NTERM_A) A(NTERM_A) where NROW_A
! is the number of rows in matrix A and NTERM_A are the number of nonzero terms in matrix A:

!      I_A is an array of NROW_A+1 integers that is used to specify the number of nonzero terms in rows of matrix A. That is:
!          I_A(I+1) - I_A(I) are the number of nonzero terms in row I of matrix A

!      J_A is an integer array giving the col numbers of the NTERM_A nonzero terms in matrix A

!        A is a real array of the nonzero terms in matrix A. If SYM_A='Y' then only the terms on, and above, the diag are stored.
!          array A is not used in this subr - it is used in the actual matrix multiplication routine.

! Input  matrix B is stored in compressed col storage (CCS) format using arrays J_B(NCOL_B+1), I_B(NTERM_B), B(NTERM_B) where NCOL_B
! is the number of columns in matrix B and NTERM_B are the number of nonzero terms in matrix B:

!      J_B is an array of NCOL_B+1 integers that is used to specify the number of nonzero terms in rows of matrix A. That is:
!          J_B(I+1) - J_B(I) are the number of nonzero terms in col I of matrix B

!      I_B is an integer array giving the row numbers of the NTERM_B nonzero terms in matrix B

!        B is a real array of the nonzero terms in matrix B. All of the nonzero terms of B must be included (i.e., no SYM option)
!          array B is not used in this subr - it is used in the actual matrix multiplication routine.

! This subr determines the storage required (NTERM_C) for array J_C (and C) from the data in arrays I_A, J_A, I_B and J_B. This
! information is used to allocate memory for arrays J_C and C prior to calling subr MATMULT_SSS so that it can it can do the sparse
! matrix multiply

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  MATMULT_SSS_NTERM_BEGEND
      USE SPARSE_ALG_ARRAYS, ONLY     :  J_AROW
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG

      USE MATMULT_SSS_NTERM_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'MATMULT_SSS_NTERM'
      CHARACTER(LEN=*), INTENT(IN)    :: MAT_A_NAME             ! Name of matrix A
      CHARACTER(LEN=*), INTENT(IN)    :: MAT_B_NAME             ! Name of matrix B
      CHARACTER(LEN=*), INTENT(IN)    :: MAT_C_NAME             ! Name of matrix C
      CHARACTER(LEN=*), INTENT(IN)    :: SYM_A                  ! ='Y' if matrix A is input sym (terms on and above diag only)
      CHARACTER(LEN=*), INTENT(IN)    :: SYM_B                  ! ='Y' if matrix B is input sym (terms on and above diag only)

      INTEGER(LONG), INTENT(IN )      :: NCOL_B                 ! Number of cols in input matrix B
      INTEGER(LONG), INTENT(IN )      :: NROW_A                 ! Number of rows in input matrix A
      INTEGER(LONG), INTENT(IN )      :: NTERM_A                ! Number of nonzero terms in input matrix A
      INTEGER(LONG), INTENT(IN )      :: NTERM_B                ! Number of nonzero terms in input matrix B
      INTEGER(LONG), INTENT(IN )      :: I_A(NROW_A+1)          ! I_A(I+1) - I_A(I) = no. terms in row I of matrix A
      INTEGER(LONG), INTENT(IN )      :: J_A(NTERM_A)           ! Col no's for nonzero terms in matrix A
      INTEGER(LONG), INTENT(IN )      :: J_B(NCOL_B+1)          ! J_B(I+1) - J_B(I) = no. terms in row I of matrix B
      INTEGER(LONG), INTENT(IN )      :: I_B(NTERM_B)           ! Row no's for nonzero terms in matrix B
      INTEGER(LONG), INTENT(OUT)      :: AROW_MAX_TERMS         ! Max number of terms in any row of A
      INTEGER(LONG), INTENT(OUT)      :: NTERM_C                ! Number of nonzero terms in output matrix C
      INTEGER(LONG)                   :: A_ROW_COLJ_BEG(NROW_A) ! jth term is row number in A where col j nonzeros begin
      INTEGER(LONG)                   :: A_ROW_COLJ_END(NROW_A) ! jth term is row number in A where col j nonzeros end
      INTEGER(LONG)                   :: A_NTERM_ROW_I          ! Number of terms in row I of matrix A
      INTEGER(LONG)                   :: A_COL_NUM              ! A col number in matrix A
      INTEGER(LONG)                   :: A_ROW_BEG              ! Index into array I_A where a row of matrix A begins
      INTEGER(LONG)                   :: A_ROW_END              ! Index into array I_A where a row of matrix A ends
      INTEGER(LONG)                   :: B_NTERM_COL_J          ! Number of terms in col J of matrix B
      INTEGER(LONG)                   :: B_ROW_NUM              ! A row number in matrix B
      INTEGER(LONG)                   :: B_COL_BEG              ! Index into array J_B where a col of matrix B begins
      INTEGER(LONG)                   :: B_COL_END              ! Index into array J_B where a col of matrix B ends
      INTEGER(LONG)                   :: DELTA_NTERM_C = 0      ! Incr in NTERM_C (0,1) resulting from mult row of A times col of B
      INTEGER(LONG)                   :: I,J,K,L,II             ! DO loop indices
      INTEGER(LONG)                   :: I1,I2                  ! DO loop range
      INTEGER(LONG)                   :: NHITS                  ! Num of "hits" of terms in a col of A existing where terms in a
!                                                                 row of B exist when a row of A is multiplied by a col of B
      INTEGER(LONG)                   :: NHITS_TOT_FOR_ROW_OF_A ! Num of "hits" of terms in a col of A existing where terms in any
!                                                                 row of B exist when a row of A is multiplied by al cols of B
      INTEGER(LONG)                   :: NTERM_AROW             ! Max number of nonzero terms in one row of A
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = MATMULT_SSS_NTERM_BEGEND

      INTRINSIC                       :: MAX

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Make sure B is stored as nonsym (all terms)

      IF (SYM_B /= 'N') THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,948) SUBR_NAME, MAT_B_NAME, SYM_B
         WRITE(F06,948) SUBR_NAME, MAT_B_NAME, SYM_B
         CALL OUTA_HERE ( 'Y' )
      ENDIF

! Initialize outputs

      AROW_MAX_TERMS = 0
      NTERM_C        = 0

! First, set up array J_AROW so it can handle the row from A that has the most terms in it

      A_ROW_BEG         = 1
      DO I=1,NROW_A
         A_NTERM_ROW_I = I_A(I+1) - I_A(I)
         A_ROW_END = A_ROW_BEG + A_NTERM_ROW_I - 1
         NTERM_AROW = 0
         IF (SYM_A == 'Y') THEN
            DO K=1,A_ROW_BEG-1
               IF (J_A(K) == I) THEN
                  NTERM_AROW = NTERM_AROW + 1
               ENDIF
            ENDDO
         ENDIF
         DO K=A_ROW_BEG,A_ROW_END
            NTERM_AROW = NTERM_AROW + 1
         ENDDO
         IF (NTERM_AROW > AROW_MAX_TERMS) THEN
            AROW_MAX_TERMS = NTERM_AROW
         ENDIF
         A_ROW_BEG = A_ROW_END + 1
      ENDDO
      NTERM_AROW = 0
                                                           ! Allocate integer vector to hold the column no's of the terms in the row
      CALL ALLOCATE_SPARSE_ALG ( 'J_AROW', AROW_MAX_TERMS, 0, SUBR_NAME)

      IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '1', '   ' )

! Create arrays that give the row numbers at which col j begins and ends. When SYM_A = 'Y'', we have to get
! terms for the A matrix that are not explicitly in A. This is done by getting A terms in the column (above the
! diagonal of a row) as well as the explicit terms from A that are there from the diagonal out to the end of the row. The two
! arrays: A_ROW_COLJ_BEG and A_ROW_COLJ_END are used to aid in getting the terms in the column above the diagonal.
! A_ROW_COLJ_BEG is an array that gives, for each col of A, the starting row number of nonzero terms in that column.
! A_ROW_COLJ_END is an array that gives, for each col of A, the ending   row number of nonzero terms in that column.
! The span: A_ROW_COLJ_BEG to A_ROW_COLJ_END is used when we search for terms in the columns.
! We only need A_ROW_COLJ_BEG and A_ROW_COLJ_END when A is input symmetric and MATOUT is not to be output as symmetric

      IF (SYM_A == 'Y') THEN

         CALL ROW_AT_COLJ_BEGEND ( MAT_A_NAME, NROW_A, NROW_A, NTERM_A, I_A, J_A, A_ROW_COLJ_BEG, A_ROW_COLJ_END )

      ENDIF

! Now count the terms that go into C, using values put into AROW and J_AROW for each row of A. This is done to facilitate the
! SYM option for matrix A.

      NHITS   = 0                                          ! Initialize variables used in the matrix multiplication in the DO I loop
      NHITS_TOT_FOR_ROW_OF_A = 0
      A_ROW_BEG  = 1
i_do: DO I=1,NROW_A                                        ! Matrix multiply loop. Range over the rows in A

         A_NTERM_ROW_I = I_A(I+1) - I_A(I)                 ! Number of terms in matrix A in row I
         IF (A_NTERM_ROW_I == 0) CYCLE i_do
         A_ROW_END = A_ROW_BEG + A_NTERM_ROW_I - 1         ! A_ROW_BEG to A_ROW_END is the range of indices of terms in A for row I
         IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '2', '   ' )

         DO K=1,AROW_MAX_TERMS                             ! Null J_AROW and AROW each time we begin a new row of A
            J_AROW(K) = 0
         ENDDO
                                                           ! Formulate J_AROW
         NTERM_AROW = 0                                    ! 1st, look for terms that would be in this row, but are not, due to SYM
         IF (SYM_A == 'Y') THEN
            DO K=1,A_ROW_BEG-1
               IF (J_A(K) == I) THEN
                  NTERM_AROW = NTERM_AROW + 1
                  I1 = A_ROW_COLJ_BEG(I)
                  I2 = A_ROW_COLJ_END(I)
                  DO II=I1,I2
                     IF ((K >= I_A(II)) .AND. (K < I_A(II+1))) THEN
                        J_AROW(NTERM_AROW) = II
                        IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '5', ' #2' )
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
         DO K=A_ROW_BEG,A_ROW_END                          ! 2nd, get terms from this row of A from the diagonal out
            NTERM_AROW = NTERM_AROW + 1
            J_AROW(NTERM_AROW) = J_A(K)
            IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '6', ' #1' )
         ENDDO

         B_COL_BEG = 1
j_do:    DO J=1,NCOL_B                                     ! J loops over the number of columns in B
            B_NTERM_COL_J = J_B(J+1) - J_B(J)
            IF (B_NTERM_COL_J == 0) CYCLE j_do
            B_COL_END = B_COL_BEG + B_NTERM_COL_J - 1

k_do:       DO K=1,NTERM_AROW                              ! The following 2 loops det  the ij-th term of C
               A_COL_NUM = J_AROW(K)
l_do:          DO L=B_COL_BEG,B_COL_END
                  B_ROW_NUM = I_B(L)
                  IF (A_COL_NUM == B_ROW_NUM) THEN
                     NHITS = NHITS + 1
                     NHITS_TOT_FOR_ROW_OF_A = NHITS_TOT_FOR_ROW_OF_A + 1
                     DELTA_NTERM_C = 1
                  ENDIF
               ENDDO l_do
            ENDDO k_do
            B_COL_BEG = B_COL_END + 1

            NTERM_C = NTERM_C + DELTA_NTERM_C
            IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '7', '   ' )
            DELTA_NTERM_C = 0
            NHITS = 0

         ENDDO j_do

         IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) THEN
            WRITE(F06,*)
         ENDIF

         IF (NHITS_TOT_FOR_ROW_OF_A == 0) THEN
             IF ((DEBUG(83) == 1) .OR. (DEBUG(83) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '8', '   ' )
         ENDIF
         NHITS_TOT_FOR_ROW_OF_A = 0
         A_ROW_BEG = A_ROW_END + 1

      ENDDO i_do

      IF ((DEBUG(84) == 1) .OR. (DEBUG(84) == 3)) CALL MATMULT_SSS_NTERM_DEB ( '9', '   ' )

      CALL DEALLOCATE_SPARSE_ALG ( 'J_AROW' )

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  948 FORMAT(' *ERROR   948: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' MATRIX B WHICH IS "',A,'" MUST BE STORED AS SYM_B = "N" FOR THIS SUBR TO WORK BUT IS STORED AS',      &
                           ' SYM_B = "',A,'"')

! ##################################################################################################################################

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE MATMULT_SSS_NTERM_DEB ( WHICH, ALG )

      CHARACTER( 3*BYTE), INTENT(IN)  :: ALG                    ! Which algorithm is used (#1 for terms above diag when SYM_A='Y'
!                                                                 or #2 for terms in row from diag out)
      CHARACTER( 1*BYTE)              :: WHICH                  ! Decides what to print out for this call to this subr

! **********************************************************************************************************************************
      IF      (WHICH == '1') THEN

         WRITE(F06,*)
         WRITE(F06,2011)
         WRITE(F06,2012)
         WRITE(F06,2013)
         WRITE(F06,2014) MAT_A_NAME, MAT_B_NAME, MAT_C_NAME
         WRITE(F06,2016) NROW_A, NTERM_A, NCOL_B, NTERM_B
         IF (SYM_A == 'Y') THEN
            WRITE(F06,2017)
         ELSE
            WRITE(F06,2018)
         ENDIF
         WRITE(F06,2019)
         WRITE(F06,*)

      ELSE IF (WHICH == '2') THEN

         WRITE(F06,2021)
         WRITE(F06,2022) I
         WRITE(F06,2023) I, I, A_ROW_BEG, A_ROW_END
         WRITE(F06,2024)
         WRITE(F06,*)

      ELSE IF (WHICH == '3') THEN

      ELSE IF (WHICH == '4') THEN

      ELSE IF (WHICH == '5') THEN

         WRITE(F06,2051) ALG,               K , II, J_A(K) , NTERM_AROW, J_AROW(NTERM_AROW)

      ELSE IF (WHICH == '6') THEN

         WRITE(F06,2061) ALG, K, I, J_A(K),                  NTERM_AROW, J_AROW(NTERM_AROW)

      ELSE IF (WHICH == '7') THEN

         IF (NHITS > 0) THEN
            IF (J == 1) THEN
               WRITE(F06,*)
            ENDIF
            WRITE(F06,2071) I, J, NHITS, NTERM_C
         ENDIF

      ELSE IF (WHICH == '8') THEN

         WRITE(F06,2081) I
         WRITE(F06,*)

      ELSE IF (WHICH == '9') THEN

         WRITE(F06,*)
         WRITE(F06,2021)
         WRITE(F06,2091) MAT_C_NAME, NTERM_C
         WRITE(F06,9100)
         WRITE(F06,*)

      ENDIF

! **********************************************************************************************************************************
 2011 FORMAT(' __________________________________________________________________________________________________________________',&
             '_________________'                                                                                ,//,&
             ' ::::::::::::::::::::::::::::::::::::START DEBUG(84) OUTPUT FROM SUBROUTINE MATMULT_SSS_NTERM::::::::::::::::::::::',&
              ':::::::::::::::::',/)

 2012 FORMAT(' SSS SETUP FOR SPARSE MATRIX MULTIPLY ROUTINE: Determine memory required for sparse Compressed Row Storage (CRS)'   ,&
' formatted matrix C',/,' --------------------------------------------',/,' resulting from the multiplication of sparse CRS'      ,&
' matrix A and sparse Compressed Column Storage (CCS) formatted matrix B',/)

 2013 FORMAT(' A may be stored as symmetric (only terms on and above the diagonal) or with all nonzero terms included.'         ,/,&
' Matrix B must be stored with all nonzero terms.',/)

 2014 FORMAT(40X,' The name of CRS formatted matrix A is: ',A                                                                   ,/,&
             40x,' The name of CRS formatted matrix B is: ',A                                                                   ,/,&
             40x,' The name of CRS formatted matrix C is: ',A,/)

 2016 FORMAT(36X,' Matrix A has ',I8,' rows and '  ,I12,' nonzero terms'                                                        ,/,&
             36X,' Matrix B has ',I8,' cols and '  ,I12,' nonzero terms'                                                        ,/,&
             36X,' Matrix C will have                       20 nonzero terms*'                                                  ,/,&
             22X,' *(as detrmined by subr MATMULT_SSS_NTERM which had to have been run prior to this subr',/)

 2017 FORMAT(' Matrix A was flagged as being a symmetric CRS array (only nonzero terms on and above the diagonal)',/)

 2018 FORMAT(' Matrix A was flagged as a CRS array that contains all nonzero terms',/)

 2019 FORMAT(                                                                                &
' In order to handle symmetric A matrices, which do not have all terms in a row (only those from the diag out), array J_AROW is',/,&
' used. J_AROW is a 1D array that contains the column numbers for all nonzero terms in one row of matrix A including those that',/,&
' may not actually be in array A due to symmetry storage. This array is used in a simulated multiply of one row of matrix A',/    ,&
' times columns of matrix B to determine whether there are any "hits" of terms in the row of A with terms in a column of B. In' ,/,&
' this manner, the number of nonzero terms (NTERM_C) that will be in array C can be estimated. The estimate is exact except for',/,&
' the case where a row of A multiplies a column of B and has "hits" that turn out to sum to zero when subr MATMULT_SSS is run' ,//,&
' Alg #1 (below) gets data for array J_AROW directly from the Compressed Row Storage (CRS) format of array A'                   ,/,&
' Alg #2 (below) is only needed if matrix A is input as symmetric (only terms on and above the diagonal) and gets terms for'    ,/,&
'         J_AROW from column I of matrix A while working on row I of matrix A. These are the terms that would be below the diag',/,&
'         in matrix A but are not explicitly in the array due to symmetry storage',//,                                             &
' For each row of matrix A, the following shows the development of array AROW and the result of multiplying AROW times matrix B',/,&
' to get one row of result matrix C. Output is only given for non null rows of matrix A and non null cols of matrix B',/)

 2021 FORMAT(' ******************************************************************************************************************',&
              '*****************')

 2022 FORMAT(30X,' W O R K I N G   O N   R O W ',I8,'   O F   O U T P U T   M A T R I X   C',/)

 2023 FORMAT(' Simulate the multiplication of row ',I8,' of matrix A times all cols of matrix B to get row ',I8,' of matrix C'  ,/,&
             ' This row of A begins in array A(K) at index K  = ',I8,' and ends at index K = ',I8,//)

 2024 FORMAT(11X,'Data from input array A           Data below diag of matrix A not in array A             Data for array AROW' ,/,&
 9X,'--------------------------          ------------------------------------------             -------------------'            ,/,&
 ' Alg     Index       Row     Col                      Index       Row     Col                         Index    Col No'           &
,/,13X,'K         I    J_A(K)                        K        II    J_A(K)                           M  J_AROW(M)')

 2051 FORMAT(1X,A3,10X,10X,10X,I25,I10,I10,I28,I11)

 2061 FORMAT(1X,A3,I10,I10,I10,25X,10X,10X,I28,I11)

 2071 FORMAT(' Row ',I8,' of A times col ',I8,' of B results in ',I8,' hits and accum number of nonzero terms in C = ',I9)

 2081 FORMAT(' There were no "hits" of any terms in row ',I8,' of matrix A multiplying terms from any col in matrix B')

 2091 FORMAT(' Total number of nonzero terms needed for allocating arrays J_C and C for matrix C = ',A,' is NTERM_C = ',I8,//)

 9100 FORMAT(' :::::::::::::::::::::::::::::::::::::END DEBUG(84) OUTPUT FROM SUBROUTINE MATMULT_SSS_NTERM:::::::::::::::::::::::',&
              ':::::::::::::::::'                                                                                ,/,&
             ' __________________________________________________________________________________________________________________',&
             '_________________',/)

! **********************************************************************************************************************************

      END SUBROUTINE MATMULT_SSS_NTERM_DEB

      END SUBROUTINE MATMULT_SSS_NTERM


text:
program square_ex

  use mathiness, only: square
  implicit none

  integer, parameter :: x = 5
  integer :: r

  r = square(x)

  write(*,*) "Square example"
  write(*,*) "A number: ", x
  write(*,*) "Its square: ", r

end program square_ex

max_stars_repo_name: mdpiper/fortran-explorations

max_issues_repo_name: mdpiper/fortran-explorations

max_forks_repo_name: mdpiper/fortran-explorations


text:
      PROGRAM CKDRIV
      implicit none
      CHARACTER*72 CHEMKIN_input, THERMO_input, TRANLIB_input, log_file,
     &     CHEMKIN_linking_file, TRANLIB_c_file, file_list
C
      integer LIN, LINCK, LOUT, LTHRM, LINMC, LTRAN, LNAME
      PARAMETER (LIN = 34, LINCK = 35, LOUT = 36, LTHRM = 37,
     &     LINMC = 38, LTRAN = 39, LNAME=40)
C
      integer IDIM, KDIM, MDIM, MAXORD, MAXSP, MAXTB, MAXTP
      PARAMETER
     +  (IDIM = 1000,
     +   KDIM = 200,
     +   MDIM = 20,
     +   MAXORD = 10,
     +   MAXSP = 12,
     +   MAXTB = 20,
     +   MAXTP = 3)
C
      integer LCWORK, LIWORK, LLWORK, LRWORK
      PARAMETER
     +  (LCWORK = KDIM + MDIM,
     +   LIWORK = IDIM * (27 + MAXORD + 2 * MAXSP + MAXTB)
     +      + KDIM * (4 + MDIM),
     +   LLWORK = KDIM,
     +   LRWORK = IDIM * (33 + MAXORD + MAXSP + MAXTB)
     +      + 2 * KDIM * (4 * MAXTP - 3) + MDIM)
C
      integer LCMCWK, LIMCWK, LLMCWK, LRMCWK
      PARAMETER
     +  (LCMCWK = LCWORK + KDIM,
     +   LIMCWK = LIWORK + 3 * KDIM,
     +   LLMCWK = LLWORK,
     +   LRMCWK = LRWORK + 300 + KDIM * (48 + 4 * KDIM + 8 * MAXTP))

      CHARACTER*16 C(LCMCWK)
      INTEGER I(LIMCWK)
      LOGICAL L(LLMCWK)
      DOUBLE PRECISION R(LRMCWK)
C
      namelist / files / CHEMKIN_input, THERMO_input, TRANLIB_input, log_file,
     &     CHEMKIN_linking_file, TRANLIB_c_file

      call getarg(1,file_list)
      OPEN (LNAME, STATUS = 'OLD', FORM = 'FORMATTED', FILE = trim(file_list))
      read(LNAME,files)
      close(LNAME)

      OPEN (LIN, STATUS = 'OLD', FORM = 'FORMATTED',
     +   FILE = trim(CHEMKIN_input))
      OPEN (LTHRM, STATUS = 'OLD', FORM = 'FORMATTED',
     +   FILE = trim(THERMO_input))
      OPEN (LTRAN, STATUS = 'OLD', FORM = 'FORMATTED',
     +   FILE = trim(TRANLIB_input))
      OPEN (LINCK, STATUS = 'UNKNOWN', FORM = 'FORMATTED',
     +   FILE = trim(CHEMKIN_linking_file))
      OPEN (LINMC, STATUS = 'UNKNOWN', FORM = 'FORMATTED',
     +   FILE = trim(TRANLIB_c_file))
      OPEN (LOUT, STATUS = 'UNKNOWN', FORM = 'FORMATTED',
     +   FILE = trim(log_file))
C
      CALL CKINTP
     +  (MDIM, KDIM, IDIM, MAXTP, MAXSP, MAXTB, LIN, LOUT, LTHRM, LINCK,
     +   MAXORD, LIWORK, I, LRWORK, R, LCWORK, C, L)

      CALL TRANFT
     +  (LINMC, LINCK, LTRAN, MAXTP, LIMCWK, LRMCWK, LCMCWK, I, R,
     +   C)

      CLOSE (LIN)
      CLOSE (LTRAN)
      CLOSE (LTHRM)
      CLOSE (LINCK)
      CLOSE (LINMC)
      CLOSE (LOUT)

      END

ext: f

max_stars_repo_name: marient/PelePhysics

max_issues_repo_name: marient/PelePhysics

max_forks_repo_name: marient/PelePhysics


text:
This page is experiencing technical difficulties. Go take a long walk and come back later.

20050906 09:57:13 nbsp Cool!  I like your webpage, by the way.  And the story you linked about the analog PDA is hilarious. Users/CarlMcCabe

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      recursive subroutine STCO( vertex )
!----------------------------------------------------------------
!	... Permutate sparse matrix in order to find the
!           the strongly connected blocks in a lower triangular
!           block form (LTBF)
!----------------------------------------------------------------

      use SP_MODS

      implicit none

!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: vertex

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i, j, vn
      integer :: adj_vertex

      nb = nb + 1
      lowlink(vertex) = nb
      number(vertex) = nb
      sp = sp + 1
      vstack(sp) = vertex
      do vn = rp(vertex),rp(vertex+1)-1
	 adj_vertex = ci(vn)
	 if( number(adj_vertex) == 0 ) then
	    call STCO( adj_vertex )
	    lowlink(vertex) = MIN( lowlink(vertex),lowlink(adj_vertex) )
	 else if( number(adj_vertex) < number(vertex) ) then
            do j = 1,sp
	       if( vstack(j) == adj_vertex ) then
		  lowlink(vertex) = MIN( lowlink(vertex),number(adj_vertex) )
		  exit
	       end if
	    end do
	 end if
      end do

      if( lowlink(vertex) == number(vertex) ) then
	 blkcnt = blkcnt + 1
	 stcoblk(blkcnt) = pp + 1
	 do
	    if( sp == 0 ) then
	       exit
	    end if
	    vn = vstack(sp)
	    if( number(vn) >= number(vertex) ) then
	       sp = sp - 1
	       pp = pp + 1
	       perm(pp) = vn
	       blkmemcnt(blkcnt) = blkmemcnt(blkcnt) + 1
	    else
	       exit
	    end if
	 end do
      end if

      end subroutine STCO

      recursive subroutine STCO_PAT( vertex, order )
!----------------------------------------------------------------
!	... Permutate sparse matrix in order to find the
!           the strongly connected blocks in a lower triangular
!           block form (LTBF)
!----------------------------------------------------------------

      use SP_MODS

      implicit none

!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: vertex, order

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i, j, vn
      integer :: adj_vertex

      nb = nb + 1
      lowlink(vertex) = nb
      number(vertex) = nb
      sp = sp + 1
      vstack(sp) = vertex
      do vn = 1,order
	 if( matrix(vertex,vn) ) then
	    adj_vertex = vn
	    if( number(adj_vertex) == 0 ) then
	       call STCO_PAT( adj_vertex, order )
	       lowlink(vertex) = MIN( lowlink(vertex),lowlink(adj_vertex) )
	    else if( number(adj_vertex) < number(vertex) ) then
               do j = 1,sp
	          if( vstack(j) == adj_vertex ) then
		     lowlink(vertex) = MIN( lowlink(vertex),number(adj_vertex) )
		     exit
	          end if
	       end do
	    end if
	 end if
      end do

      if( lowlink(vertex) == number(vertex) ) then
	 blkcnt = blkcnt + 1
	 stcoblk(blkcnt) = pp + 1
	 do
	    if( sp == 0 ) then
	       exit
	    end if
	    vn = vstack(sp)
	    if( number(vn) >= number(vertex) ) then
	       sp = sp - 1
	       pp = pp + 1
	       perm(pp) = vn
	       blkmemcnt(blkcnt) = blkmemcnt(blkcnt) + 1
	    else
	       exit
	    end if
	 end do
      end if

      end subroutine STCO_PAT

      subroutine DIAG_MARK( order, matrix, perm )
!---------------------------------------------------------------------------
!	... Find permuatation via diagonal markowitz to produce
!           near optimal LU fillin
!---------------------------------------------------------------------------

      implicit none

!---------------------------------------------------------------------------
!	... Dummy args
!---------------------------------------------------------------------------
      integer, intent(in)  :: order
      integer, intent(out) :: perm(order)
      logical, intent(in)  :: matrix(order,order)

!---------------------------------------------------------------------------
!	... Local variables
!---------------------------------------------------------------------------
      integer :: row, rowp1, col, maxrow
      integer :: beta, cnt
      integer :: i, j
      logical :: holder(order)
      logical :: pattern(order,order)

      cnt = COUNT( matrix )
      perm(:order) = (/ (row,row=1,order) /)        ! no permutations
      pattern = matrix
      do row = 1,order-1
	 rowp1 = row + 1
	 beta = (order - 1)**2
	 maxrow = row
	 do col = row,order
	    cnt = (COUNT(pattern(col,row:order)) - 1) * &
	          (COUNT(pattern(row:order,col)) - 1)
	    if( cnt < beta ) then
	       beta = cnt
	       maxrow = col
	    end if
	 end do
	 if( maxrow /= row ) then
!----------------------------------------------------------------
!	... Row and column permuatation
!----------------------------------------------------------------
	    holder(row:order) = pattern(row,row:order)
	    pattern(row,row:order) = pattern(maxrow,row:order)
	    pattern(maxrow,row:order) = holder(row:order)
	    holder(row:order) = pattern(row:order,row)
	    pattern(row:order,row) = pattern(row:order,maxrow)
	    pattern(row:order,maxrow) = holder(row:order)
	    beta = perm(row)
	    perm(row) = perm(maxrow)
	    perm(maxrow) = beta
	 end if
!----------------------------------------------------------------
!	... Now do "symbolic" decomposition on sub-matrix
!----------------------------------------------------------------
	 do col = rowp1,order
	    if( pattern(row,col) ) then
	       pattern(rowp1:order,col) = pattern(rowp1:order,row) .or. &
					  pattern(rowp1:order,col)
	    end if
	 end do
      end do

      end subroutine DIAG_MARK

      subroutine SYM_FAC( order, matrix, adds, mults )
!----------------------------------------------------------------
!	... Do a symbolic LU decomposition on sparse matrix
!----------------------------------------------------------------

      implicit none

!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: order                     ! matrix order
      logical, intent(inout) :: matrix(order,order)    ! sparse matrix
      integer, intent(out) :: adds(2), mults(2)        ! operation counts

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i, im1, j, rcnt, ccnt

      adds = 0; mults = 0;
      do i = 2,order
	 im1 = i - 1
	 ccnt = COUNT( matrix(i:order,im1) )
	 rcnt = COUNT( matrix(im1,i:order) )
	 mults(1) = mults(1) + ccnt * (1 + rcnt)
	 mults(2) = mults(2) + ccnt
	 do j = i,order
	    if( matrix(im1,j) ) then
	       adds(1) = adds(1) + COUNT( matrix(i:order,j) .and. matrix(i:order,im1) )
	       matrix(i:order,j) = matrix(i:order,j) .or. matrix(i:order,im1)
	    end if
	 end do
      end do

      do i = order,2,-1
	 ccnt = COUNT( matrix(1:i-1,i) )
	 mults(2) = mults(2) + ccnt
      end do
      adds(2) = mults(2)

      end subroutine SYM_FAC

      subroutine GEN_PAT( order, pattern, rp, ci )
!----------------------------------------------------------------
!	... Generate sparsity pattern from row storage
!----------------------------------------------------------------

      implicit none
!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: order
      integer, intent(in) :: rp(order+1)
      integer, intent(in) :: ci(*)
      logical, intent(out) :: pattern(order,order)

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i, row

      pattern = .false.
      do i = 1,order
	 do row = rp(i),rp(i+1)-1
	    pattern(i,ci(row)) = .true.
	 end do
      end do

      end subroutine GEN_PAT

      subroutine PERM_MAT( order, matrix, perm)
!----------------------------------------------------------------
!	... Matrix row and column permutation
!----------------------------------------------------------------

      implicit none

!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: order
      integer, intent(in) :: perm(order)
      logical, intent(inout) :: matrix(order,order)

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i
      logical :: copy(order,order)

      copy = matrix
!----------------------------------------------------------------
!	... Row permutation
!----------------------------------------------------------------
      do i = 1,order
	 if( perm(i) /= i ) then
	    copy(i,:order) = matrix(perm(i),:order)
	 end if
      end do
      matrix = copy
!----------------------------------------------------------------
!	... Col permutation
!----------------------------------------------------------------
      do i = 1,order
	 if( perm(i) /= i ) then
	    matrix(:order,i) = copy(:order,perm(i))
	 end if
      end do

      end subroutine PERM_MAT

      subroutine DRAW_MAT( order, matrix )
!----------------------------------------------------------------
!	... Draw matrix sparsity
!----------------------------------------------------------------

      implicit none

!----------------------------------------------------------------
!	... Dummy args
!----------------------------------------------------------------
      integer, intent(in) :: order
      logical, intent(in) :: matrix(order,order)

!----------------------------------------------------------------
!	... Local variables
!----------------------------------------------------------------
      integer :: i, row
      character(len=1)    :: matrow(order)
      character(len=32)   :: frmt
      character(len=132)  :: line

      write(*,*) ' '
      line = ' '
      do i = 10,order,10
	 write(line(5+i*2:5+i*2),'(i1)') i/10
      end do
      write(*,'(a)') line(:LEN_TRIM(line))
      line = ' '
      do i = 1,order
	 write(line(5+i*2:5+i*2),'(i1)') MOD( i,10 )
      end do
      write(*,'(a)') line(:LEN_TRIM(line))
      write(*,*) ' '
      frmt = '(1x,i2,2x,'
      i = LEN_TRIM( frmt ) + 1
      if( order < 10 ) then
	 write(frmt(i:i),'(i1)') order
      else
	 write(frmt(i:i+1),'(i2)') order
      end if
      frmt(LEN_TRIM(frmt)+1:) = '(1x,a))'
      do row = 1,order
	 matrow = ' '
	 where( matrix(row,:order) )
	    matrow(:order) = 'X'
	 endwhere
	 write(*,frmt) row, matrow
      end do

      end subroutine DRAW_MAT

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
ccoorree -- System Administration

Format of a core-dump file
#iinncclluuddee <ssyyss/ccoorree.hh>

When  a   process  terminates  abnormally   because  a  it   encounters  an
unrecoverable  error  or  receives  an  asynchronous  signal  from  another
process, COHERENT tries to write a  copy of its image in memory into a file
called ccoorree. You can examine this file with the debugger ddbb and other tools
to try to determine what went wrong.

The structure cchh_iinnffoo appears at the  head of a ccoorree file.  The header file
ccoorree.hh defines it as follows:

    struct  ch_info {
        unsigned short ch_magic;
        unsigned int ch_info_len;
    };

Field cchh_mmaaggiicc  is always  set to  the constant CCOORREE_MMAAGGIICC.  This ``magic''
value signifies  to COHERENT that  this is a core  file.  Field cchh_iinnffoo_lleenn
gives a count of information bytes  in the core file, including the cchh_iinnffoo
structure itself.

If the value of cchh_iinnffoo_lleenn exceeds the size of the cchh_iinnffoo structure, this
indicates that  data follow the  cchh_iinnffoo structure.  These  data follow the
cchh_iinnffoo structure,  and are in  the form of a  ccoorree_pprroocc structure.  Header
file <ssyyss/ccoorree.hh> defines this structure as follows:

    struct core_proc {
        gregset_t cp_registers;
        int cp_signal_number;
        struct _fpstate cp_floating_point;
        dregset_t cp_debug_registers;
    };

This substitutes  for a dump of  the uu area, whose  information is reserved
for the kernel alone.

This is  followed by an image  of each process segment.   The data for each
segment consists of  the following: a header, which is  a structure of type
ccoorree_sseegg; ccss_ppaatthhlleenn bytes of data that give the path name of the file from
the segment data originated; and ccss_dduummppeedd bytes of core-image data.

ccoorree.hh defines strucutre ccoorree_sseegg as follows:

     struct core_seg {
          size_t cs_pathlen;            /* length of path name */
          off_t cs_dumped;              /* dumped size in bytes */
          caddr_t cs_base;              /* virtual base address */
          off_t cs_size;                /* full size in bytes */
          unsigned long cs_reserved[8];
     };

The order  of the segments is  the text segment first (if  it is present --
usually it  is omitted), followed by  the data segment, and  then the stack
segment.  The  contents of the text segment can  usually be identified from
the program being debugged.  The patchable kernel variable DDUUMMPP_TTEEXXTT allows
the  COHERENT kernel  to  dump text  segments  as well  as  data and  stack
segments.

Patchable  kernel variable  DDUUMMPP_LLIIMM  sets the  maximum size  of a  segment
within a  ccoorree file.  The  system uses this  limit to keep  keep core files
from getting out of hand.

_S_e_e _A_l_s_o
AAddmmiinniisstteerriinngg CCOOHHEERREENNTT, ccoorree.hh, ssiiggnnaall(), wwaaiitt()

_D_i_a_g_n_o_s_t_i_c_s
COHERENT will not write ccoorree if  that file already exists as a non-ordinary
file or if  there is more than one link to  it.  The 0200 bit in the status
returned to the parent process by wwaaiitt() indicates a successful dump.

For  a list  of signals  that automatically  trigger a  core dump,  see the
Lexicon entry for ssiiggnnaall().

max_stars_repo_name: gspu/Coherent

max_issues_repo_name: gspu/Coherent

max_forks_repo_name: gspu/Coherent


text:
MODULE GALACTIC
USE dflib

IMPLICIT NONE

 real(8), parameter, private :: Leo = 282.85948083 ! °
 real(8), parameter, private :: AL  = Leo-90.0_8 ! °
 real(8), parameter, private :: L0  =  32.931918056 ! °
 real(8), parameter, private :: si  = 0.88998807641_8 ! sin 62.871748611°
 real(8), parameter, private :: ci  = 0.45598379779_8 ! cos 62.871748611°

 real(8), parameter :: Pi  = 3.1415926535897932384626433832795_8

 real(8), parameter, private :: Leo1 = 282.25 ! °
 real(8), parameter, private :: L01  =  33 ! °
 real(8), parameter, private :: si1  = 0.88781538514_8 ! sin 62.6°
 real(8), parameter, private :: ci1  = 0.46019978478_8 ! cos 62.6°

 ! Прецессионные константы
REAL(8), PARAMETER :: M_1975 = 3.07374_8*15.0_8/3600.0_8, N_1975 = 20.0405_8/3600.0_8
REAL(8), PARAMETER :: M_1950 = 3.07327_8*15.0_8/3600.0_8, N_1950 = 20.0426_8/3600.0_8

CONTAINS

INTEGER(4) FUNCTION Round(x)
REAL(8), INTENT(IN) :: x

 Round=INT(ANINT(x))

END FUNCTION

 Subroutine Galaxy(a,d,l,b)
 real(8),intent(in) :: a,d
 real(8),intent(out) :: l,b
 real(8) :: al,sa,ca,sd,cd

  al=a-Leo
  sa=dsind(al);  ca=dcosd(al)
  sd=dsind(d);   cd=dcosd(d)
  b=dasind(sd*ci-cd*si*sa)
  l=datan2d(sd*si+cd*ci*sa,cd*ca)+L0
  if (l<0) then
   l=l+360.0_8
  endif
 end subroutine

Subroutine Galaxy1950(a,d,l,b)
 real(8),intent(in) :: a,d
 real(8),intent(out) :: l,b
 real(8) :: al,sa,ca,sd,cd

  al=a-Leo1
  sa=dsind(al);  ca=dcosd(al)
  sd=dsind(d);   cd=dcosd(d)
  b=dasind(sd*ci1-cd*si1*sa)
  l=datan2d(sd*si1+cd*ci1*sa,cd*ca)+L01
  if (l<0) then
   l=l+360.0_8
  endif
 end subroutine

 Subroutine Equatorial(l,b,a,d)
 real(8),intent(in) :: l,b
 real(8),intent(out) :: a,d
 real(8) cb,sb,l1,cl,sl
 l1=l-L0
 cb=dcosd(b);  sb=dsind(b)
 cl=dcosd(l1); sl=dsind(l1)

 d=dasind(cb*si*sl+sb*ci)
 a=datan2d(cb*cl,sb*si-cb*ci*sl)+AL
 if (a>360.0_8) then
   a=a-360.0_8
 end if
 end subroutine

 Subroutine GalaxMu(mua,mud,l,b,d, mul,mub)
 real(8),intent(in) :: mua,mud,l,b,d
 real(8),intent(out) :: mul,mub
 real(8) :: cd,sfi,cfi
   cd=dcosd(d);
   sfi=si*dcosd(l-L0)/cd
   cfi=(dcosd(b)*ci-dsind(b)*si*dsind(l-L0))/cd
   mul= cfi*mua+sfi*mud
   mub=-sfi*mua+cfi*mud
 end subroutine

 REAL(4) function Ranorm(s)
 ! s - дисперсия
 real(4), intent(in) :: s
 real(4) :: a
 integer, save :: z=37843827
 integer(4) ::   i;
  a=0.0
  DO i=1,12
   a=a+ran(z)
  END DO
  Ranorm=(a-6.0)*s
 end FUNCTION RANORM

 REAL(8) FUNCTION PRECRA(RA,DE, Dt)
 ! Прецессия по прямому восхождению в градусах
 ! RA, DE - в градусах
 ! Dt - в годах
 REAL(8), INTENT(IN) :: RA,DE,Dt
  PRECRA=(M_1975+N_1975*DSIND(RA)*DTAND(DE))*Dt
 END FUNCTION

 REAL(8) FUNCTION PRECDE(RA,DE, Dt)
 ! Прецессия по склонению восхождению в градусах
 REAL(8), INTENT(IN) :: RA,DE,Dt
  PRECDE=N_1975*DCOSD(RA)*Dt
 END FUNCTION

 REAL(8) FUNCTION RHO(L1,B1,L2,B2)
 ! Расстояние между объектами в градусах
 REAL(8), INTENT(IN) :: L1,B1,L2,B2
  RHO=DACOSD( DSIND(B1)*DSIND(B2)+ DCOSD(B1)*DCOSD(B2)*DCOSD(L1-L2) )
 END FUNCTION

END MODULE

max_stars_repo_name: movsesyanpv/gaiaread

max_issues_repo_name: movsesyanpv/gaiaread

max_forks_repo_name: movsesyanpv/gaiaread


text:
! Ode solver
!
! Adapted from Numerical Recipes 77

submodule (solver) solver_sm

    real(8), parameter :: pi = 4*atan(1.0d0)
    real(8), parameter :: two_pi = 2*pi
    real(8), allocatable, target :: xp_data(:), yp_data(:,:)
    real(8), allocatable, target :: x_error(:,:)

contains
!
!   Cash-Karp - an adaptive Runge-Kutta method
!   see https://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method
    module function rk_cash_karp(h,t,x,fp) result(xout)
        real(8), intent(in)             :: h, t
        real(8), intent(in), contiguous :: x(:,:)
        procedure(f_p), pointer         :: fp
        real(8)                         :: xout(size(x,1),size(x,2))
        real(8), dimension(size(x,1),size(x,2)) :: k1, k2, k3, k4, k5, k6
        real(8) :: a2,a3,a4,a5,a6,b21,b31,b32,b41,b42,b43,b51,b52,b53
        real(8) :: b54,b61,b62,b63,b64,b65,c1,c3,c4,c6,dc1,dc3,dc4,dc5,dc6

        parameter (a2=1/5.0d0, a3=3/10.0d0, a4=3/5.0d0, a5=1.0d0 ,a6=7/8.0d0,                               &
                   b21=1/5.0d0,                                                                             &
                   b31=3/40.0d0, b32=9/40.0d0,                                                                  &
                   b41=.3, b42=-.9, b43=1.2,                                                                &
                   b51=-11./54., b52=2.5, b53=-70./27., b54=35./27.,                                        &
                   b61=1631./55296.,b62=175./512., b63=575./13824.,b64=44275./110592.,b65=253./4096.,       &
                   c1=37./378., c3=250./621., c4=125./594., c6=512./1771.,                                  &
                   dc1=c1-2825./27648., dc3=c3-18575./48384., dc4=c4-13525./55296., dc5=-277./14336., dc6=c6-.25)

        k1 = h*fp(t         , x   )
        k2 = h*fp(t + a2*h  , x + b21*k1  )
        k3 = h*fp(t + a3*h  , x + b31*k1 + b32*k2  )
        k4 = h*fp(t + a4*h  , x + b41*k1 + b42*k2 + b43*k3   )
        k5 = h*fp(t + a5*h  , x + b51*k1 + b52*k2 + b53*k3 + b54*k4   )
        k6 = h*fp(t + a6*h  , x + b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5   )
        xout = x + c1*k1 + c3*k3 + c4*k4 + c6*k6
        x_error = dc1*k1 + dc3*k3 + dc4*k4 + dc5*k5 + dc6*k6
    end function rk_cash_karp

    module function rk_get_error() result(r)
        real(8), dimension(:,:), pointer :: r
        if (allocated(x_error)) then
            r => x_error
        else
            nullify(r)
        end if
    end function rk_get_error

    ! Explicit 4th order Runge-Kutta method
    module function rk_exp4(h,t,x,fp) result(xout)
        real(8), intent(in)             :: h, t
        real(8), intent(in), contiguous :: x(:,:)
        procedure(f_p), pointer :: fp
        real(8)                         :: xout(size(x,1),size(x,2))
        real(8), dimension(size(x,1),size(x,2)) :: k1, k2, k3, k4
        k1   = h*fp(t    , x     )
        k2   = h*fp(t+h/2, x+k1/2)
        k3   = h*fp(t+h/2, x+k2/2)
        k4   = h*fp(t+h  , x+k3  )
        xout = x + k1/6 + k2/3 + k3/3 + k4/6
    end function rk_exp4

    ! Explicit 4th order Runge-Kutta method - the 3/8 rule
    module function rk_exp38(h,t,x,fp) result(xout)
        real(8), intent(in)             :: h, t
        real(8), intent(in), contiguous :: x(:,:)
        procedure(f_p), pointer :: fp
        real(8)                         :: xout(size(x,1),size(x,2))
        real(8), dimension(size(x,1),size(x,2)) :: k1, k2, k3, k4
        k1   = h*fp(t      , x     )
        k2   = h*fp(t+h/3  , x + k1/3)
        k3   = h*fp(t+2*h/3, x - k1/3 + k2)
        k4   = h*fp(t+h    , x + k1 - k2 + k3  )
        xout = x + k1/8 + 3*k2/8 + 3*k3/8 + k4/8
    end function rk_exp38

    module subroutine four1(data,nn,isign)
        integer, intent(in) :: isign, nn
        real(8) :: data(2*nn)
        integer :: i,istep,j,m,mmax,n
        real(8) :: temp(2)
        real(8) :: theta,wi,wpi,wpr,wr,wtemp
        n = 2*nn
        j = 1
        do i=1,n,2
            if(j > i) then
              temp = data(j:j+1)
              data(j:j+1) = data(i:i+1)
              data(i:i+1) = temp
            endif
            m = n/2
            do
                if((m < 2) .or. (j <= m)) exit
                j = j-m
                m = m/2
            end do
            j = j+m
        end do

        mmax = 2
        do
            if (n <= mmax) exit
            istep = 2*mmax
            theta = two_pi/(isign*mmax)
            wpr = -2*sin(0.5d0*theta)**2
            wpi = sin(theta)
            wr = 1
            wi = 0
            do m=1,mmax,2
                do i=m,n,istep
                    j=i+mmax
                    temp = [wr*data(j)-wi*data(j+1), wr*data(j+1)+wi*data(j)]
                    data(j:j+1) = data(i:i+1) - temp
                    data(i:i+1) = data(i:i+1) + temp
                end do
                wtemp = wr
                wr = wr*wpr - wi*wpi+wr
                wi = wi*wpr + wtemp*wpi+wi
            end do
            mmax = istep
        end do
        return
    end subroutine four1
!
    module subroutine realft(data, isign)
        integer, intent(in)     :: isign
        real(8), intent(inout) :: data(:)

        integer  :: i,i1,i2,i3,i4,n2p3, n
        real(8) :: c1,c2,h1i,h1r,h2i,h2r
        real(8) :: theta,wi,wpi,wpr,wr,wtemp

        n = size(data)
        theta = 2*pi/n
        c1 = 0.5
        if (isign == 1) then
            c2 = -0.5d0
            call four1(data,n/2,+1)
        else
            c2 = 0.5d0
            theta = -theta
        end if
        wpr = -2*sin(0.5d0*theta)**2
        wpi = sin(theta)
        wr = 1 + wpr
        wi = wpi
        n2p3 = n+3
        do i=2,n/4
            i1 = 2*i-1
            i2 = i1+1
            i3 = n2p3-i2
            i4 = i3+1
            h1r = c1*(data(i1)+data(i3))
            h1i = c1*(data(i2)-data(i4))
            h2r = -c2*(data(i2)+data(i4))
            h2i = c2*(data(i1)-data(i3))
            data(i1) = h1r+wr*h2r-wi*h2i
            data(i2) = h1i+wr*h2i+wi*h2r
            data(i3) = h1r-wr*h2r+wi*h2i
            data(i4) = -h1i+wr*h2i+wi*h2r
            wtemp = wr
            wr = wr*wpr-wi*wpi+wr
            wi = wi*wpr+wtemp*wpi+wi
        end do
        if (isign == 1) then
            h1r = data(1)
            data(1) = h1r+data(2)
            data(2) = h1r-data(2)
        else
            h1r = data(1)
            data(1) = c1*(h1r+data(2))
            data(2) = c1*(h1r-data(2))
            call four1(data,n/2,-1)
        end if
        return
    end subroutine realft

end submodule solver_sm

max_stars_repo_name: sgeard/code

max_issues_repo_name: sgeard/code

max_forks_repo_name: sgeard/code


!----------------------------------------------------------------------------
! This file is part of UCLALES.
!
! UCLALES is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! UCLALES is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1999-2007, Bjorn B. Stevens, Dep't Atmos and Ocean Sci, UCLA
!----------------------------------------------------------------------------
!
! PRSS: Pressure Solver:  Solves the anelastic or bousinessq system
! for the pressure using a fractional step method, which is implemented
! using fft's and a tri-diagonal solver in the vertical
!
MODULE prss

   USE defs, ONLY : pi
   IMPLICIT NONE

CONTAINS
   !
   !----------------------------------------------------------------------
   ! Subroutine poisson: called by timesteping driver to invert the
   ! poisson equation for pressure and apply the velocity tendencies.
   !
   SUBROUTINE poisson

      USE grid, ONLY : nxp, nyp, nzp, dtlt, dxi, dyi, dzm, dzt, a_up,                     &
                       a_uc, a_ut, a_vp, a_vc, a_vt, a_wp, a_wc, a_wt, a_press, a_pexnr,  &
                       th00, dn0, wsavex, wsavey
      USE stat, ONLY : fill_scalar, sflg
      USE util, ONLY : ae1mm

      COMPLEX, ALLOCATABLE :: s1(:,:,:)
      REAL    :: mxdiv, awpbar(nzp)
      INTEGER :: ix,iy

      ix = max(1,nxp-4)
      iy = max(1,nyp-4)

      ALLOCATE (s1(ix,iy,nzp))
      s1 = 0.0
      !
      ! -------
      ! Do first step of asselin filter, first saving corrlations of
      ! tendencies for TKE budget on statistics timestep
      !

      CALL asselin(1)
      CALL apl_tnd(nzp,nxp,nyp,a_up,a_vp,a_wp,a_ut,a_vt,a_wt,dtlt)
      !
      ! ------
      ! Pressure Solve
      !
      CALL poiss(nzp,nxp,nyp,ix,iy,a_up,a_vp,a_wp,a_pexnr,a_press,dn0,th00,dzt, &
                 dzm,dxi,dyi,dtlt,s1,wsavex,wsavey)
      CALL ae1mm(nzp,nxp,nyp,a_wp,awpbar)
      !
      ! -------
      ! Do second step of asselin filter, first saving corrlations of
      ! tendencies for TKE budget on statistics timestep, note that the
      ! old centered velocity resides in up,vp,wp after the second step
      ! of the asselin filter, hence the pressure correlation terms in the
      ! tke budget include the effects of time filtering
      !
      CALL asselin(2)
      CALL velocity_bcs
      CALL get_diverg(nzp,nxp,nyp,ix,iy,s1,a_up,a_vp,a_wp,dn0,dzt,dxi,dyi,  &
                      dtlt,mxdiv)

      IF (sflg) THEN
         CALL get_diverg(nzp,nxp,nyp,ix,iy,s1,a_up,a_vp,a_wp,dn0,dzt,dxi,dyi,  &
                         dtlt,mxdiv)
         CALL fill_scalar(2,mxdiv)
         CALL prs_cor(nzp,nxp,nyp,a_pexnr,a_up,a_vp,a_wp,dzm,dxi,dyi,th00)
         CALL chk_tsplt(nzp,nxp,nyp,a_up,a_vp,a_wp,a_uc,a_vc,a_wc)
      END IF
      DEALLOCATE (s1)

   END SUBROUTINE poisson

   !
   ! --------------------------------------------------------------------
   ! Subroutine apl_tnd: applies tendencies to velocity field
   !
   SUBROUTINE apl_tnd(n1,n2,n3,u,v,w,ut,vt,wt,dtlt)

      USE util, ONLY : velset

      INTEGER :: n1,n2,n3,i,k,j
      REAL    :: u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)
      REAL    :: ut(n1,n2,n3),vt(n1,n2,n3),wt(n1,n2,n3),dtlt,dt

      dt = dtlt*2.

      DO j = 1, n3
         DO i = 1, n2
            DO k = 2, n1-1
               w(k,i,j) = w(k,i,j)+wt(k,i,j)*dt
               u(k,i,j) = u(k,i,j)+ut(k,i,j)*dt
               v(k,i,j) = v(k,i,j)+vt(k,i,j)*dt
            END DO
         END DO
      END DO

      CALL velset(n1,n2,n3,u,v,w)

   END SUBROUTINE apl_tnd
   !
   ! --------------------------------------------------------------------
   ! Subroutine poiss: called each timestep to evaluate the pressure
   ! in accordance with the anelastic continuity equation, and then apply
   ! the pressure to the velocity terms for three dimensional flow,
   ! cyclic in x and y.  pp and pc are used as scratch arrays in the
   ! call to trdprs.  pp is filled with its diagnostic value in fll_prs
   !
   SUBROUTINE poiss(n1,n2,n3,ix,iy,u,v,w,pp,pc,dn0,th00,dzt,dzm,dx,dy,  &
                    dtlt,s1,wsvx,wsvy)

      USE util, ONLY  : velset, get_fft_twodim

      INTEGER :: n1,n2,n3,ix,iy
      REAL    :: pp(n1,n2,n3),pc(n1,n2,n3),dmy
      REAL    :: u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)
      REAL    :: wsvx(1:),wsvy(1:),dn0(n1),dzt(n1),dzm(n1),dx,dy,dtlt,th00
      COMPLEX :: s1(ix,iy,n1)

      CALL get_diverg(n1,n2,n3,ix,iy,s1,u,v,w,dn0,dzt,dx,dy,dtlt,dmy)
      CALL get_fft_twodim(ix,iy,n1,s1,wsvx,wsvy,-1)
      CALL trdprs(n1,ix,iy,s1,dn0,dzt,dzm,dx,dy)
      CALL get_fft_twodim(ix,iy,n1,s1,wsvx,wsvy,+1)

      CALL fll_prs(n1,n2,n3,ix,iy,pp,s1)
      CALL prs_grd(n1,n2,n3,pp,u,v,w,dzm,dx,dy,dtlt)

      CALL velset(n1,n2,n3,u,v,w)

      pp(:,:,:) = pp(:,:,:)/th00
      pc(:,:,:) = pp(:,:,:)

   END SUBROUTINE poiss
   !
   ! --------------------------------------------------------------------
   ! Subroutine get_diverg: gets velocity divergence and puts it into
   ! a complexvalue array for use in pressure calculation
   !
   SUBROUTINE get_diverg(n1,n2,n3,ix,iy,s1,u,v,w,dn0,dz,dx,dy,dt,mxdiv)

      INTEGER, INTENT (in)  :: n1,n2,n3,ix,iy
      REAL, INTENT (in)     :: dz(n1),dn0(n1),dx,dy,dt
      REAL, INTENT (in)     :: u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)
      REAL, INTENT (out)    :: mxdiv
      COMPLEX, INTENT (out) :: s1(ix,iy,n1)

      INTEGER :: k,i,j,l,m
      REAL    :: xf,yf,zf,wf1,wf2,dtv,dti

      s1(:,:,:) = (0.0,0.0)
      dtv = dt*2.
      dti = 1./dtv
      m = 0
      DO j = 3, n3-2
         m = m+1
         l = 0
         DO i = 3, n2-2
            l = l+1
            DO k = 2, n1-1
               wf1 = 0.5*(dn0(k+1)+dn0(k))
               wf2 = 0.5*(dn0(k)+dn0(k-1))
               IF (k == 2 )   wf2 = 0.
               IF (k == n1-1) wf1 = 0.
               xf = dn0(k)*dx*dti
               yf = dn0(k)*dy*dti
               zf = dti*dz(k)
               s1(l,m,k) = (wf1*w(k,i,j)-wf2*w(k-1,i,j))*zf                  &
                           +(v(k,i,j)-v(k,i,j-1))*yf+(u(k,i,j)-u(k,i-1,j))*xf
            END DO
         END DO
      END DO
      !
      ! save mxdiv to a statistics array, no reduction necessary as this is done
      ! in post processing
      !
      mxdiv = maxval(REAL(s1))

   END SUBROUTINE get_diverg
   !
   !----------------------------------------------------------------------
   ! Subroutine fll_prs: writes the pressure to the appropriate array
   !
   SUBROUTINE fll_prs(n1,n2,n3,ix,iy,pp,s1)

      USE mpi_interface, ONLY : cyclics,cyclicc

      INTEGER :: n1,n2,n3,ix,iy,k,i,j,l,m,req(16)
      REAL    :: pp(n1,n2,n3)
      COMPLEX :: s1(ix,iy,n1)

      pp(:,:,:) = 0.0
      DO k = 2, n1-1
         l = 0
         DO i = 3, n2-2
            l = l+1
            m = 0
            DO j = 3, n3-2
               m = m+1
               pp(k,i,j) = REAL(s1(l,m,k))
            END DO
         END DO
      END DO
      CALL cyclics(n1,n2,n3,pp,req)
      CALL cyclicc(n1,n2,n3,pp,req)

   END SUBROUTINE fll_prs
   !
   !---------------------------------------------------------------------
   ! TRDPRS: solves for the wave number (l,m) component of
   ! pressure in a vertical column using a tri-diagonal solver.
   !
   SUBROUTINE trdprs(n1,ix,iy,s1,dn0,dzt,dzm,dx,dy)

      USE mpi_interface, ONLY : yoffset, nypg, xoffset, wrxid, wryid, nxpg
      USE util, ONLY          : tridiff

      INTEGER, INTENT (in)    :: n1,ix,iy
      REAL, INTENT (in)       :: dn0(n1),dzt(n1),dzm(n1),dx,dy
      COMPLEX, INTENT (inout) :: s1(ix,iy,n1)

      REAL    :: ak(ix,n1),dk(ix,n1),bk(ix,n1),ck(ix,n1)
      REAL    :: xk(ix,n1),yk(ix,n1),wv(ix,iy)

      INTEGER :: k,l,m
      REAL    :: fctl,fctm,xl,xm,af,cf
      fctl = 2.*pi/float(nxpg-4)
      fctm = 2.*pi/float(nypg-4)

      DO l = 1, ix
         IF(l+xoffset(wrxid) <= (nxpg-4)/2+1) THEN
            xl = float(l-1+xoffset(wrxid))
         ELSE
            xl = float(l-(nxpg-4)-1+xoffset(wrxid))
         END IF

         DO m = 1, iy
            IF(m+yoffset(wryid) <= (nypg-4)/2+1) THEN
               xm = float(m-1+yoffset(wryid))
            ELSE
               xm = float(m-(nypg-4)-1+yoffset(wryid))
            END IF
            wv(l,m) = 2.*((cos(fctl*xl)-1.)*dx*dx + (cos(fctm*xm)-1.)*dy*dy)
         END DO
      END DO

      IF(wrxid == 0 .AND. wryid == 0 ) THEN
         wv(1,1) = 0.
      END IF
      !
      ! configure vectors for tri-diagonal solver
      !
      DO m = 1, iy
         DO k = 2, n1-1
            af = (dn0(k)+dn0(k-1))*.5
            cf = (dn0(k+1)+dn0(k))*.5
            IF (k == 2   ) af = 0.
            IF (k == n1-1) cf =0.
            DO l = 1, ix
               ak(l,k) = dzt(k)*dzm(k-1)*af
               bk(l,k) = s1(l,m,k)
               ck(l,k) = dzt(k)*dzm(k)*cf
               dk(l,k) = dn0(k)*wv(l,m)-(ak(l,k)+ck(l,k))
            END DO
         END DO
         !
         ! solve for fourier components, x_k, given a tri-diagonal matrix of the
         ! form a_k x_k-1 + d_k x_k + c_k x_k+1 = b_k.  y_k is a scratch array.
         !
         CALL tridiff(ix,n1-1,ix,ak,dk,ck,bk,xk,yk)

         DO k = 2, n1-1
            DO l = 1, ix
               IF (m+yoffset(wryid)+l+xoffset(wrxid) > 2) bk(l,k) = aimag(s1(l,m,k))
               IF (m+yoffset(wryid)+l+xoffset(wrxid) > 2) s1(l,m,k) = xk(l,k)
            END DO
         END DO

         CALL tridiff(ix,n1-1,ix,ak,dk,ck,bk,xk,yk)

         DO k = 2, n1-1
            DO l = 1, ix
               IF (m+yoffset(wryid)+l+xoffset(wrxid) > 2)                &
                  s1(l,m,k) = cmplx(REAL(s1(l,m,k)),xk(l,k))
            END DO
         END DO

      END DO

   END SUBROUTINE trdprs
   !
   !---------------------------------------------------------------------
   ! Subroutine Prs_grd: apply the pressure gradient term
   !
   SUBROUTINE prs_grd(n1,n2,n3,p,u,v,w,dz,dx,dy,dtlt)

      INTEGER, INTENT (in) :: n1,n2,n3
      REAL, INTENT (in)    :: p(n1,n2,n3),dz(n1),dx,dy,dtlt
      REAL, INTENT (inout) :: u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)

      INTEGER :: i,j,k

      DO j = 1, n3-1
         DO i = 1, n2-1
            DO k = 2, n1-1
               IF(k /= n1-1)w(k,i,j) = w(k,i,j)-dz(k)*2.*dtlt*(p(k+1,i,j)-p(k,i,j))
               u(k,i,j) = u(k,i,j)-dx*2.*dtlt*(p(k,i+1,j)-p(k,i,j))
               v(k,i,j) = v(k,i,j)-dy*2.*dtlt*(p(k,i,j+1)-p(k,i,j))
            END DO
         END DO
      END DO

   END SUBROUTINE prs_grd
   !
   !---------------------------------------------------------------------
   ! Subroutine Prs_cor: correlate the pressure tendency with velocity
   ! field for TKE budget
   !
   SUBROUTINE prs_cor(n1,n2,n3,p,u,v,w,dz,dx,dy,th00)

      USE stat, ONLY : updtst
      USE util, ONLY : get_cor

      INTEGER, INTENT (in) :: n1,n2,n3
      REAL, INTENT (in)    :: p(n1,n2,n3),dz(n1),dx,dy,th00
      REAL, INTENT (in)    :: u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3)

      REAL, DIMENSION (n2,n3) :: pgx, pgy, pgz, ufld, vfld, wfld
      REAL                    :: fx, fy, fz, v1da(n1), v1db(n1), v1dc(n1)
      INTEGER                 :: i,j,k,ip1,jp1

      v1da = 0.0
      v1db = 0.0
      v1dc = 0.0

      DO k = 2, n1-1
         fx = dx*th00
         fy = dy*th00
         fz = dz(k)*th00
         DO j = 1, n3
            DO i = 1, n2
               ip1 = min(n2,i+1)
               jp1 = min(n3,j+1)
               pgx(i,j) = -fx*(p(k,ip1,j)-p(k,i,j))
               pgy(i,j) = -fy*(p(k,i,jp1)-p(k,i,j))
               pgz(i,j) = -fz*(p(k+1,i,j)-p(k,i,j))
               ufld(i,j) = u(k,i,j)
               vfld(i,j) = v(k,i,j)
               wfld(i,j) = w(k,i,j)
            END DO
         END DO
         v1da(k) = get_cor(1,n2,n3,1,ufld,pgx)
         v1db(k) = get_cor(1,n2,n3,1,vfld,pgy)
         v1dc(k) = get_cor(1,n2,n3,1,wfld,pgz)
      END DO
      CALL updtst(n1,'prs',1,v1da,1)
      CALL updtst(n1,'prs',2,v1db,1)
      CALL updtst(n1,'prs',3,v1dc,1)

   END SUBROUTINE prs_cor
   !
   ! --------------------------------------------------------------------
   ! Subroutine chk_tsplt
   !
   SUBROUTINE chk_tsplt(n1,n2,n3,up,vp,wp,uc,vc,wc)

      INTEGER, INTENT (in) :: n1,n2,n3
      REAL, INTENT (in)    :: up(n1,n2,n3),vp(n1,n2,n3),wp(n1,n2,n3)
      REAL, INTENT (in)    :: uc(n1,n2,n3),vc(n1,n2,n3),wc(n1,n2,n3)

      REAL :: wmx,umx,vmx

      wmx = maxval((wp-wc)/(wp+wc+1.e-5))
      umx = maxval((up-uc)/(up+uc+1.e-5))
      vmx = maxval((vp-vc)/(vp+vc+1.e-5))

   END SUBROUTINE chk_tsplt
   !
   !----------------------------------------------------------------------
   ! Subroutine Asselin:  Applies the asselin filter in two stages
   ! depending on the value of iac
   !
   SUBROUTINE asselin(iac)

      USE grid, ONLY : a_up,a_vp,a_wp,a_uc,a_vc,a_wc, nxyzp, runtype

      INTEGER :: iac
      INTEGER, SAVE :: ncall = 0

      IF (runtype == 'HISTORY') ncall = 1

      CALL predict(nxyzp,a_uc,a_up,iac,ncall)
      CALL predict(nxyzp,a_vc,a_vp,iac,ncall)
      CALL predict(nxyzp,a_wc,a_wp,iac,ncall)
      IF (iac == 2) ncall = ncall+1

   END SUBROUTINE asselin
   !
   !----------------------------------------------------------------------
   ! Subroutine predict:  This subroutine advances the leapfrog terms
   ! in two stages.  It applies the filter equation:
   !
   !         a(n) = a(n) + eps * (a(n-1) - 2*a(n) + a(n+1))
   !
   ! the first stage of the filter applies all but the a(n+2) term.
   ! the second stage renames the variables and applies this term, i.e.,
   ! a(n+1) -> a(n), a(n+2) -> a(n+1).  Note that for iac=2 ap=a(n+2)
   ! because the tendencies have been updated in pressure solver.  Durran,
   ! in his text cites values of eps of 0.2 for convective cloud models,
   ! we seem to get by with eps=0.1, perhaps because of the coupling
   ! provided by the staggered forward step.
   !
   SUBROUTINE predict(npts,ac,ap,iac,iflag)

      INTEGER :: m,npts,iac,iflag
      REAL    :: ac(npts),ap(npts),af(npts),epsu
      REAL, PARAMETER :: eps = 0.1

      epsu = eps
      IF (iflag == 0) epsu = 0.5
      IF (iac == 1) THEN
         DO m = 1, npts
            ac(m)=ac(m)+epsu*(ap(m)-2.*ac(m))
         END DO
      ELSE IF (iac == 2) THEN
         DO m = 1, npts
            af(m) = ap(m)
            ap(m) = ac(m)+epsu*af(m)
            ac(m) = af(m)
         END DO
      END IF

   END SUBROUTINE predict
   !
   !----------------------------------------------------------------------
   ! Subroutine Velocity_bcs: Applies boundary conditions on veolicities
   !
   SUBROUTINE velocity_bcs

      USE grid, ONLY : a_up,a_vp,a_wp,a_uc,a_vc,a_wc,a_pexnr,a_press,           &
                       nxp, nyp, nzp, dzm
      USE util, ONLY : velset, sclrset

      CALL velset(nzp,nxp,nyp,a_up,a_vp,a_wp)
      CALL velset(nzp,nxp,nyp,a_uc,a_vc,a_wc)
      CALL sclrset('grad',nzp,nxp,nyp,a_pexnr,dzm)
      CALL sclrset('grad',nzp,nxp,nyp,a_press,dzm)

   END SUBROUTINE velocity_bcs

END MODULE prss


  program pertenergy

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !PROGRAM: pertenergy - POstprocess perturbation vector(s) to output energy field
!                        as a dynamic vector (in *.nc4 file).  If multiple perturbation
!                        vectors are processed, then result is the average energy field.
!
! !USES:

  use precision

  use m_mpif90,only : MP_init
  use m_mpif90,only : MP_finalize
  use m_mpif90,only : MP_comm_rank
  use m_mpif90,only : MP_comm_world

  use m_pertenergy, only : pertenergy_Set
  use m_pertenergy, only : pertenergy_Process
  use m_pertenergy, only : pertenergy_clean

  use m_pertutil, only : pertutil_setparam

  use m_die, only : MP_die, die

  implicit none

! !DESCRIPTION: This routine makes the actual perturbation vector
!
! To do:
!  \begin{itemize}
!      \item MPI part
!      \item double check this is really correct
!  \end{itemize}
!
! !REVISION HISTORY:
!
!  20Nov2002  Todling    Initial code: projection matrix assumed diagonal.
!  20Dec2002  Gelaro     Properly defined projection operator.
!  01Oct2004  Winslow    Hacked to work with adjoint initialization
!  16Aug2005  RG/EN      Adopted postStats.f90 initial code to create
!                           a stand-alone post-processing utility
!  03May2007  Gelaro     Added rcfile capability to process multiple perturbations
!                        a stand-alone post-processing utility
!  02Feb2007  Todling    Placed mpi hooks (only runs in 1 pe)
!  10Jan2007  Todling    - Remove energy-type opt(generalized)
!                        - Add end date/time to command line arg list
!  09Aug2017  Todling    Opt not to time-tag output file
!
!EOP
!-------------------------------------------------------------------------

  character(len=*), parameter :: myname= 'pertenergy'

  integer, parameter :: mfiles = 3        ! max no. of dyn-vector and txt files
  integer, parameter :: ROOT = 0 ! MPI ROOT PE

  character(len=255) :: dynfiles(mfiles)  ! input dyn filename(s)
  character(len=255) :: expid             ! name of experiment
  character(len=255) :: rcfile            ! input resource file

  integer :: knmax        ! total no. of dyn-like arrays
  logical :: pick
  logical :: notag
  integer :: myID, ier
  integer :: nymd         ! beginning (or only) date (YYYYMMDD) to process
  integer :: nhms         ! beginning (or only) time (HHMMSS) to process
  integer :: nymd_end     ! ending date (YYYYMMDD) to process
  integer :: nhms_end     ! ending time (HHMMSS) to process

! MPI Initialization
! ------------------
  call MP_init(ier)
        if(ier/=0) call MP_die(myname,'MP_init()',ier)
  call MP_comm_rank(MP_comm_world,myID,ier)
        if(ier/=0) call MP_die(myname,'MP_comm_rank()',ier)

! Parse command line
! ------------------
  call init_ ( mfiles, dynfiles, pick, &
               nymd, nhms, nymd_end, nhms_end, expid, notag, rcfile )

! Set up program options from rcfile
! ----------------------------------
  call pertenergy_Set ( MP_comm_world, root, rcfile, ier )

! Process perturbation file(s)
! ----------------------------
  call pertenergy_Process ( mfiles, dynfiles, expid, notag, nymd, nhms, nymd_end, nhms_end )

! Clean up
! --------
  call pertenergy_clean( )

  print *, 'Program done'
  call MP_finalize(ier)
        if(ier/=0) call MP_die(myname,'MP_finalized()',ier)
 contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!BOP
! !IROUTINE: Init_ --- Initialize pertenergy
!
! !DESCRIPTION: parses command line.
!
! !INTERFACE:
!
      subroutine Init_ ( mfiles, dynfiles, pick, &
                         nymd, nhms, nymd_end, nhms_end, expid, notag, rcfile )

      implicit NONE

      integer,       intent(in)  :: mfiles  ! max. number of eta files
                                            ! dynamics file names (eta)
      character*255, intent(out) :: dynfiles(mfiles)
      character*255, intent(out) :: expid
      logical,       intent(out) :: pick
      integer,       intent(out) :: nymd
      integer,       intent(out) :: nhms
      integer,       intent(out) :: nymd_end
      integer,       intent(out) :: nhms_end
      logical,       intent(out) :: notag

!
! !REVISION HISTORY:
!       15oct2004  Winslow/RT  Initial code, based on similar codes
!       08Aug2005  RG/EN       Adopted initadj.f90 for this stand-alone application
!       10Aug2008  Todling     Add various handles and safe-guards
!
!EOP
!BOC

      character*4, parameter :: myname = 'init'

      integer iret, i, iarg, argc, iargc, vectype, nfiles
      character(len=255) :: argv, res
      character(len=255) :: pertfile, ref_state_file, energy_file, rcfile
      character*10 str
      real eps_eer
      integer vnorm

!     Defaults
!     --------
      ref_state_file = 'ref_state.nc4'
      energy_file ='pertenergy.eta'    ! default filename for output pert
      pertfile    = 'pert.nc4'
      pick        = .false.
      nymd        = 0
      nhms        = 0
      nymd_end    = -1
      nhms_end    = -1
      rcfile      = 'null'
      expid       = 'null'
      notag       = .false.
      eps_eer     = -999.0
      vnorm       = -999

!     Parse command line
!     ------------------
      argc =  iargc()
      if ( argc .lt. 1 ) call usage()

      iarg = 0
      nfiles = 0

      do i = 1, 32767
         iarg = iarg + 1
         if ( iarg .gt. argc ) exit
         call GetArg ( iarg, argv )

         if (index(argv,'-g5' ) .gt. 0 ) then
            vectype = 5
         else if (index(argv,'-o' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, energy_file )
            nfiles = nfiles + 1
         else if (index(argv,'-pert' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, pertfile )
            nfiles = nfiles + 1
         else if (index(argv,'-ref' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, ref_state_file )
            nfiles = nfiles + 1
         else if (index(argv,'-pick' ) .gt. 0 ) then
            pick = .true.
            if ( iarg+2 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) nymd
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) nhms
         else if (index(argv,'-eps_eer' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) eps_eer
         else if (index(argv,'-epick' ) .gt. 0 ) then
            pick = .true.
            if ( iarg+2 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) nymd_end
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) nhms_end
            rcfile = 'pertenergy.rc'   ! reset default name for rc file
         else if (index(argv,'-expid' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, expid )
         else if (index(argv,'-notag' ) .gt. 0 ) then
            notag = .true.
         else if (index(argv,'-vnorm' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, argv )
            read(argv,*) vnorm
         else if (index(argv,'-rc' ) .gt. 0 ) then
            if ( iarg+1 .gt. argc ) call usage()
            iarg = iarg + 1
            call GetArg ( iArg, rcfile )
         end if

             if ( nfiles .gt. mfiles ) call die(myname,'too many eta files')
      end do

      if ( nfiles .lt. 2 .and. rcfile=='null' ) then
           print *
           print *, myname,': BOTH SENS and REF files are required'
           print *
           call usage()
      endif

      if ( expid=='null' .and. rcfile/='null' ) call usage()

      call pertutil_setparam ( 'vectype', vectype )
      if( vectype==5 ) then
          call pertutil_setparam ( 'tname', 'vT_' )
          call pertutil_setparam ( 'wgrid', 'a'   )
      endif

      if ( nymd_end==nhms_end .and. nymd_end==-1 ) then
           nymd_end= nymd
           nhms_end= nhms
      endif
      if ( rcfile=='null' ) then
          dynfiles(1) = pertfile
          dynfiles(2) = ref_state_file
          if (eps_eer < 0.0) then
             print*
             print*, myname, ': when no RC-file command line eps_eer must be used'
             print*
             call usage()
          else
             call pertutil_setparam ( 'eps_eer', real(eps_eer,r8) )
          endif
          if (vnorm < 0.0) then
             print*
             print*, myname, ': when no RC-file command line vnorm must be used'
             print*
             call usage()
          else
             call pertutil_setparam ( 'vnorm', vnorm )
          endif
      endif
      dynfiles(3) = energy_file

      end subroutine Init_

!.................................................................

      subroutine usage()
      print *
      print *, '  -------------------------------------------------------------------'
      print *, '  pertenergy - generates a perturbation vector from adj/sensitivity outputs '
      print *, '  -------------------------------------------------------------------'
      print *
      print *
      print *,'Usage: '
      print *,'  To process single file'
      print *,'    pertenergy.x [-o energy ] -pert pert -ref refstate -pick nymd nhms '
      print *
      print *,'  To process multiple files'
      print *,'    pertenergy.x [-o energy ] -expid expid -rc RCfile '
      print *, '                             -pick  nymdb nhmsb '
      print *, '                            -epick  nymde nhmse '
      print *
      print *, 'where'
      print *
      print *, '-o      energy       prefix for output nc4 filename'
      print *
      print *, '-ref    refstate    reference state dynamic vector filename'
      print *, '                    (default: ref_state.nc4)'
      print *
      print *, '-pick   nymd nhms   date and time at the reference state vector'
      print *, '                    (default: read single-time files in arg list)'
      print *
      print *, '-epick  nymde nhmse end date and time of files to process '
      print *, '                    (default: ignore it)'
      print *
      print *, '-pert   pert        filename of perturbation file'
      print *, '                    (default: pert.nc4)'
      print *
      print *, '-rc     rcfile      resource file to control processing, overrides -ref and -pert'
      print *, '                    REQUIRED to process multiple perturbation files'
      print *, '                    OPTIONAL to process a single perturbation file'
      print *, '                    (default name: pertenegy.rc)'
      print *
      print *, ' The following must be specified when RC not used:'
      print *
      print *, '-eps_eer EPS_EER    q-coeff for TE norm (default: none)'
      print *, '-vnorm   0,1,or 2   vertical weight in TE norm (default: none)'
      print *
      print *, ' Last updated: 29 Aug 2020; Todling '
      print *
      call MP_finalize(ier)
        if(ier/=0) call MP_die(myname,'MP_finalized()',ier)
      call exit(1)
      end subroutine usage

!.................................................................

 end program pertenergy


text:
MODULE EXTRA_CONSTANT
  USE CONSTANT, ONLY: C
  REAL, PARAMETER :: B = 1.0
END MODULE EXTRA_CONSTANT

MODULE REFERENCE_ADDITION
  USE ADDITION
END MODULE REFERENCE_ADDITION

max_stars_repo_name: thlux/fmodpy

max_issues_repo_name: thlux/fmodpy

max_forks_repo_name: thlux/fmodpy


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

!> Generate initial inflation files from namelist values.
!> This way an experiment can always start from a restart file
!> without having to alter the namelist between cycles 1 and 2.
!>
!> an alternative to running this program is to use the nco utilities thus:
!>
!> Here is an example using version 4.4.2 or later of the NCO tools:
!>   ncap2 -s "T=1.0;U=1.0;V=1.0" wrfinput_d01 prior_inflation_mean.nc
!>   ncap2 -s "T=0.6;U=0.6;V=0.6" wrfinput_d01 prior_inflation_sd.nc'

program fill_inflation_restart

use             types_mod, only : r8, i8, missing_r8

use         utilities_mod, only : register_module, error_handler, E_MSG, E_ERR, &
                                  find_namelist_in_file, check_namelist_read,   &
                                  nmlfileunit, do_nml_file, do_nml_term

use       assim_model_mod, only : static_init_assim_model

use      time_manager_mod, only : time_type

use  ensemble_manager_mod, only : init_ensemble_manager, ensemble_type, &
                                  end_ensemble_manager

use   state_vector_io_mod, only : state_vector_io_init, read_state, write_state

use   state_structure_mod, only : get_num_domains

use      io_filenames_mod, only : io_filenames_init, file_info_type,       &
                                  stage_metadata_type, get_stage_metadata, &
                                  get_restart_filename, set_file_metadata, &
                                  set_io_copy_flag, READ_COPY, WRITE_COPY

use             model_mod, only : get_model_size

use     mpi_utilities_mod, only : initialize_mpi_utilities, finalize_mpi_utilities

implicit none

! version controlled file description for error handling, do not edit
character(len=*), parameter :: source   = &
   "$URL$"
character(len=*), parameter :: revision = "$Revision$"
character(len=*), parameter :: revdate  = "$Date$"

! MAX_FILES is max number of domains
integer, parameter :: MAX_FILES = 10
integer, parameter :: ss_inflate_index    = 1
integer, parameter :: ss_inflate_sd_index = 2

!------------------------------------------------------------------
! The namelist variables
!------------------------------------------------------------------

logical            :: write_prior_inf = .FALSE.
real(r8)           :: prior_inf_mean = MISSING_R8
real(r8)           :: prior_inf_sd   = MISSING_R8
logical            :: write_post_inf = .FALSE.
real(r8)           :: post_inf_mean  = MISSING_R8
real(r8)           :: post_inf_sd    = MISSING_R8
logical            :: single_file    = .FALSE.
character(len=256) :: input_state_files(MAX_FILES)  = ''
logical            :: verbose        = .FALSE.

namelist /fill_inflation_restart_nml/              &
             prior_inf_mean, prior_inf_sd,         &
             post_inf_mean, post_inf_sd, verbose,  &
             write_prior_inf, write_post_inf,      &
             input_state_files, single_file

! io variables
integer                   :: iunit, io
type(file_info_type)      :: file_info_input
type(file_info_type)      :: file_info_output
type(stage_metadata_type) :: input_restart_files
type(stage_metadata_type) :: output_restart_files

! model state variables
type(ensemble_type)   :: ens_handle
type(time_type)       :: model_time
integer(i8)           :: model_size
integer               :: ncopies = 2    ! {prior,posterior}_inf_{mean,sd}
logical               :: read_time_from_file = .true.

! counter variables
integer :: idom, imem, ndomains

! message strings
character(len=512) :: my_base, my_desc, my_stage
character(len=512) :: string1, string2, string3

character(len=256), allocatable :: file_array_input(:,:)

!======================================================================
! start of executable code
!======================================================================

call initialize_modules_used()

call find_namelist_in_file("input.nml", "fill_inflation_restart_nml", iunit)
read(iunit, nml = fill_inflation_restart_nml, iostat = io)
call check_namelist_read(iunit, io, "fill_inflation_restart_nml")

! Record the namelist values used for the run ...
if (do_nml_file()) write(nmlfileunit, nml=fill_inflation_restart_nml)
if (do_nml_term()) write(     *     , nml=fill_inflation_restart_nml)

!----------------------------------------------------------------------
! Calling static_init_assim_model() is required, which also calls
! static_init_model(), so there is no need to explicitly call it.
!----------------------------------------------------------------------

call static_init_assim_model()

!----------------------------------------------------------------------
! initialization code, model size
!----------------------------------------------------------------------

model_size = get_model_size()

!----------------------------------------------------------------------
! read/write restart files
!----------------------------------------------------------------------

! Set up the ensemble storage and read in the restart file
call init_ensemble_manager(ens_handle, ncopies, model_size)

! Allocate space for file arrays.
! Contains a matrix of files (ncopies x ndomains)
! If perturbing from a single instance the number of
! input files does not have to be ens_size but rather
! a single file (or multiple files if more than one domain)

ndomains = get_num_domains()

allocate(file_array_input(ncopies, ndomains))

file_array_input  = RESHAPE(input_state_files,  (/1,  ndomains/))

if(single_file) then
  call error_handler(E_ERR, 'fill_inflation_restart: ', &
                     'single_file not yet supported, please contact DART', &
                     source, revision, revdate)
endif

call io_filenames_init(file_info_input,             &
                       ncopies      = ncopies,      &
                       cycling      = single_file,  &
                       single_file  = single_file,  &
                       restart_files = file_array_input(:,:))

! Read the template file to get the shape of netCDF file
! and its variables. It is possible to have multiple domains
! but only require one member.
write(my_base,'(A)') 'template'
write(my_desc,'(A)') 'template file'
call set_file_metadata(file_info_input,                  &
                       cnum     = 1,                     &
                       fnames   = file_array_input(1,:), &
                       basename = my_base,               &
                       desc     = my_desc)

call set_io_copy_flag(file_info_input, &
                      cnum    = 1,     &
                      io_flag = READ_COPY)

input_restart_files = get_stage_metadata(file_info_input)

imem = 1
do idom = 1, ndomains
   write(string1, *) '- Reading File : ', &
                     trim(get_restart_filename(input_restart_files, &
                                               copy   = imem,       &
                                               domain = idom))
   call error_handler(E_MSG, 'fill_inflation_restart: ', string1,   &
                      source, revision, revdate)
enddo

call read_state(ens_handle, file_info_input, read_time_from_file, model_time)

! Initialize file output to default inflation file names
call io_filenames_init(file_info_output,           &
                       ncopies      = ncopies,     &
                       cycling      = single_file, &
                       single_file  = single_file)

if(write_prior_inf) then
   call fill_inflation_files(prior_inf_mean, prior_inf_sd, 'prior')
endif

if(write_post_inf) then
   call fill_inflation_files(post_inf_mean,  post_inf_sd,   'post')
endif

deallocate(file_array_input)

call finalize_modules_used()

!======================================================================
contains
!======================================================================

!> fill inflation values in a separate file

subroutine fill_inflation_files(inf_mean, inf_sd, stage)
real(r8),         intent(in) :: inf_mean
real(r8),         intent(in) :: inf_sd
character(len=*), intent(in) :: stage

if (inf_mean == MISSING_R8 .or. inf_sd == MISSING_R8) then
   write(string1,*) 'you must specify both inflation mean and inflation standard deviation values'
   write(string2,*) 'you have "',trim(stage),'_inf_mean = ', inf_mean,'" and '
   write(string3,*) '         "',trim(stage),'_inf_sd   = ', inf_sd,  '"     '
   call error_handler(E_MSG, 'fill_inflation_restart: ', string1,  &
                      source, revision, revdate, text2=string2, text3=string3)
   return
endif
ens_handle%copies(ss_inflate_index   , :) = inf_mean
ens_handle%copies(ss_inflate_sd_index, :) = inf_sd

write(my_stage,'(3A)') 'input_', stage, 'inf'
write(my_base, '(A)')  'mean'
write(my_desc, '(2A)') stage, ' inflation mean'
call set_file_metadata(file_info_output,    &
                       cnum     = ss_inflate_index, &
                       stage    = my_stage, &
                       basename = my_base,  &
                       desc     = my_desc)

call set_io_copy_flag(file_info_output,    &
                      cnum    = 1,         &
                      io_flag = WRITE_COPY)

write(my_base, '(A)') 'sd'
write(my_desc, '(2A)') stage, ' inflation sd'
call set_file_metadata(file_info_output,    &
                       cnum     = ss_inflate_sd_index, &
                       stage    = my_stage, &
                       basename = my_base,  &
                       desc     = my_desc)

call set_io_copy_flag(file_info_output,    &
                      cnum    = 2,         &
                      io_flag = WRITE_COPY)

output_restart_files = get_stage_metadata(file_info_output)

do idom = 1, ndomains
   ! write inflation mean
   write(string1, *) '- Writing ',trim(stage), ' Mean Inflation File : "',   &
                     trim(get_restart_filename(output_restart_files,         &
                                               copy   = ss_inflate_index,    &
                                               domain = idom)), '"'
   write(string2, *) ' with value = ', inf_mean
   call error_handler(E_MSG, 'fill_inflation_restart: ', string1,            &
                      source, revision, revdate, text2=string2)

   ! write inflation stadard deviation
   write(string1, *) '- Writing ',trim(stage), ' SD   Inflation File : "',   &
                     trim(get_restart_filename(output_restart_files,         &
                                               copy   = ss_inflate_sd_index, &
                                               domain = idom)), '"'
   write(string2, *) ' with value = ', inf_sd
   call error_handler(E_MSG, 'fill_inflation_restart: ', string1,            &
                      source, revision, revdate, text2=string2)
enddo

call write_state(ens_handle, file_info_output)

end subroutine fill_inflation_files

!----------------------------------------------------------------------
!> initialize modules that need it

subroutine initialize_modules_used()

call initialize_mpi_utilities('fill_inflation_restart')

call register_module(source,revision,revdate)

call state_vector_io_init()

end subroutine initialize_modules_used

!----------------------------------------------------------------------
!> clean up before exiting

subroutine finalize_modules_used()

call end_ensemble_manager(ens_handle)

! this must be last, and you can't print/write anything
! after this is called.
call finalize_mpi_utilities()

end subroutine finalize_modules_used

end program fill_inflation_restart

max_stars_repo_name: nancycollins/DART

max_issues_repo_name: nancycollins/DART

max_forks_repo_name: nancycollins/DART


text:
!-----------------------------------------------------------------------
!     Function:      fblin_nag
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          02/21/2012
!     Description:   This subroutine calculates the RHS of the ODE for
!                    field line following (NAG)
!-----------------------------------------------------------------------
      SUBROUTINE fblin_nag(phi,q,qdot)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE fieldlines_grid
!      USE EZspline_obj
!      USE EZspline
!-----------------------------------------------------------------------
!     Input Variables
!          phi        phi angle
!          q          (q(1),q(2)) = (R,Z)
!          qdot       dq/dt
!-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION :: phi, q(2), qdot(2)
!-----------------------------------------------------------------------
!     Local Variables
!          ier        Error flag
!          r_temp     R
!          z_temp     Z
!          phi_temp   phi (radians)
!          br_temp    br/bphi evaluated from splines
!          bz_temp    bz/bphi evaluated from splines
!          br_dot     Random force for diffusion
!-----------------------------------------------------------------------
      INTEGER :: ier
      DOUBLE PRECISION :: r_temp, phi_temp, z_temp, br_temp, bz_temp
      ! For splines
      INTEGER :: i,j,k
      REAL*8 :: xparam, yparam, zparam, hx, hy, hz, hxi, hyi, hzi
      REAL*8 :: fval(1)
      INTEGER, PARAMETER :: ict(8)=(/1,0,0,0,0,0,0,0/)
      !INTEGER, parameter :: ict(10)=(/1,0,0,0,0,0,0,0,0,0/)
      REAL*8, PARAMETER :: one = 1

!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      ier    = 0
      r_temp = q(1)
      z_temp = q(2)
      phi_temp = MOD(phi,delta_phi)
      IF (phi_temp < 0) phi_temp = delta_phi + phi_temp
      br_temp = 0.0; bz_temp = 0.0
      !CALL EZspline_isInDomain(BR_spl,r_temp,phi_temp,z_temp,ier)
      !IF (ier == 0) THEN
      IF ((r_temp >= rmin-eps1) .and. (r_temp <= rmax+eps1) .and. &
          (phi_temp >= phimin-eps2) .and. (phi_temp <= phimax+eps2) .and. &
          (z_temp >= zmin-eps3) .and. (z_temp <= zmax+eps3)) THEN
         ! Get the gridpoint info
         i = MIN(MAX(COUNT(raxis < r_temp),1),nr-1)
         j = MIN(MAX(COUNT(phiaxis < phi_temp),1),nphi-1)
         k = MIN(MAX(COUNT(zaxis < z_temp),1),nz-1)
         hx     = raxis(i+1) - raxis(i)
         hy     = phiaxis(j+1) - phiaxis(j)
         hz     = zaxis(k+1) - zaxis(k)
         hxi    = one / hx
         hyi    = one / hy
         hzi    = one / hz
         xparam = (r_temp - raxis(i)) * hxi
         yparam = (phi_temp - phiaxis(j)) * hyi
         zparam = (z_temp - zaxis(k)) * hzi
         CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                         hx,hxi,hy,hyi,hz,hzi,&
                         BR4D(1,1,1,1),nr,nphi,nz)
         !CALL R8FVTRICUB(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
         !                hx,hxi,hy,hyi,hz,hzi,&
         !                BR4D(1,1,1,1),nr,nphi,nz)
         br_temp = fval(1)
         CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                         hx,hxi,hy,hyi,hz,hzi,&
                         BZ4D(1,1,1,1),nr,nphi,nz)
         !CALL R8FVTRICUB(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
         !                hx,hxi,hy,hyi,hz,hzi,&
         !                BZ4D(1,1,1,1),nr,nphi,nz)
         bz_temp = fval(1)
      END IF
      qdot(1) = br_temp
      qdot(2) = bz_temp
      RETURN
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE fblin_nag

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      USE PRIMEBAG	!Gain access to NEXTPRIME and ISPRIME.
Calculates the largest "left-truncatable" digit sequence that is a prime number, in various bases.
      INTEGER LBASE,MANY,ENUFF	!Some sizes.
      PARAMETER (LBASE = 13, MANY = 66666, ENUFF = 66)
      INTEGER NS,START(LBASE)	!A list of single-digit prime numbers for starters.
      INTEGER NH,LH		!Counters for the horde.
      INTEGER N,HORDEN(MANY)		!Numerical value of a digit sequence.
      INTEGER*1 HORDED(ENUFF,MANY)	!Single-digit values only.
      INTEGER B,D,DB	!The base, a digit, some power of the base.
      INTEGER L		!The length of the digit sequence: DB = B**L.
      INTEGER P		!Possibly a prime number.
      INTEGER I		!A stepper.

      MSG = 6	!I/O unit number for "standard output".
      IF (.NOT.GRASPPRIMEBAG(66)) STOP "Gan't grab my file!"	!Attempt in hope.
      NS = 0	!No starters.
      P = 1	!Start looking for some primes.
    1 P = NEXTPRIME(P)	!Thus skipping non-primes.
      IF (P.LE.LBASE) THEN	!Beyond the limit?
        NS = NS + 1		!No. Count another starter.
        START(NS) = P		!Save its value.
        GO TO 1			!And seek further.
      END IF		!One could insted prepare some values, the primes being well-known.
      WRITE (MSG,2) LBASE,NS,START(1:NS)	!But, parameterisation is easy enough.
    2 FORMAT ("Working in bases 3 to ",I0," there are ",I0,	!Announce the known.
     * " single-digit primes: ",666(I0:", "))	!The : sez stop if the list is exhausted.
      WRITE (MSG,3)		!Produce a heading for the tabular output.
    3 FORMAT (/"Base Digits Count Max. Value = (in base)")

   10 DO B = 3,LBASE	!Work through the bases.
        NH = 0			!The horde is empty.
        DO I = 1,NS		!Prepare the starters for base B.
          IF (START(I).GE.B) EXIT	!Like, they're single-digits in base B.
          NH = NH + 1			!So, count another in.
          HORDEN(NH) = START(I)		!Its numerical value.
          HORDED(1,NH) = START(I)	!Its digits. Just one.
        END DO			!On to the next single-digit prime number.
        L = 0	!Syncopation. The length of the digit sequences.
        DB = 1	!The power for the incoming digit.

   20   L = L + 1	!We're about to add another digit.
        IF (L.GE.ENUFF) STOP "Too many digits!"	!Hopefully, there's room.
        DB = DB*B	!The new power of B.
        IF (DB.LE.0) GO TO 29	!Integer overflow?
        LH = NH		!The live ones, awaiting extension.
        DO I = 1,LH	!Step through each starter.
          N = HORDEN(I)	!Grab its numerical value.
          DO D = 1,B - 1	!Consider all possible lead digits.
            P = D*DB + N		!Place it at the start of the number.
            IF (P.LE.0) GO TO 29	!Oh for IF OVERFLOW ...
            IF (ISPRIME(P)) THEN	!And if it is a prime,
              IF (NH.GE.MANY) STOP "Too many sequences!"	!Add a sequence.
              NH = NH + 1			!Count in a survivor.
              HORDEN(NH) = P			!The numerical value.
              HORDED(1:L,NH) = HORDED(1:L,I)	!The digits.
              HORDED(L + 1,NH) = D		!Plus the added high-order digit.
            END IF			!So much for persistent primality.
          END DO		!On to the next lead digit.
        END DO	!On to the next starter.

        N = NH - LH		!The number of entries added to the horde.
        IF (N.GT.0) THEN	!Were there any?
          DO I = 1,MIN(LH,N)		!Yes. Overwrite the starters.
            HORDEN(I) = HORDEN(NH)		!From the tail end of the horde.
            HORDED(1:L + 1,I) = HORDED(1:L + 1,NH)	!Digit sequences as well.
            NH = NH - 1				!One snipped off.
          END DO			!Thus fill the gap at the start.
          NH = N			!The new horde count.
          LH = NH			!All are starters for the next level.
          GO TO 20			!See how it goes.
        END IF			!So much for further progress.
        GO TO 30		!But if none, done.
   29   WRITE (MSG,28) B,L,NH,DB,P	!Curses. Offer some details.
   28   FORMAT (I4,I7,I6,28X,"Integer overflow!",2I12)
        CYCLE			!Or, GO TO the end of the loop.

   30   I = MAXLOC(HORDEN(1:NH),DIM = 1)	!Finger the mostest number.
        WRITE (MSG,31) B,L,NH,HORDEN(I),HORDED(L:1:-1,I)	!Results!
   31   FORMAT (I4,I7,I6,I11," = "666(I0:"."))	!See Format 3.

      END DO		!On to the next base.
      END	!Simple enough.

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData