!
!     -------------------------------------------------------------
!      N O N R E L A T 3
!     -------------------------------------------------------------
!                                                                  *
!     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
!     OF TWO PARTICLE OPERATOR IN CASE :       N'1 = N1 - 1        *
!                                              N'2 = N2 + 1        *
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville            February 1994   *
!
      SUBROUTINE NONRELAT3(IA, IB, IC, ID, IIRE)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE OPERAT_C
      USE MEDEFN_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  11:17:07  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE nonrelat31_I
      USE ittk_I
      USE nonrelat32_I
      USE nonrelat33_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: IA
      INTEGER  :: IB
      INTEGER  :: IC
      INTEGER  :: ID
      INTEGER  :: IIRE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
      IF (IHSH <= 1) RETURN
      IF (IB == ID) THEN
         IF (IA==IB .OR. IC==IB) THEN
            IF (IA == IC) GO TO 10
            IF (IC == IB) THEN
               CALL NONRELAT31 (IC, IA, IA, IB, IC, ID, IIRE)
            ELSE
               IF (IIRE == 0) RETURN
               IF (ISOTOP == 1) THEN
                  IF (ITTK(LJ(IA),LJ(IA),1) == 0) RETURN
                  IF (ITTK(LJ(IC),LJ(IA),1) == 0) RETURN
               ENDIF
               CALL NONRELAT32 (IC, IA, IA, IB, IC, ID)
            ENDIF
         ELSE
            IF (IIRE == 0) RETURN
            CALL NONRELAT33 (IC, IA, IB, 1, IA, IB, IC, ID)
         ENDIF
      ELSE IF (IA == IC) THEN
         IF (IB==IA .OR. ID==IA) THEN
            IF (IB == ID) GO TO 10
            IF (ID == IA) THEN
               CALL NONRELAT31 (ID, IB, IA, IB, IC, ID, IIRE)
            ELSE
               IF (IIRE == 0) RETURN
               IF (ISOTOP == 1) THEN
                  IF (ITTK(LJ(IB),LJ(IB),1) == 0) RETURN
                  IF (ITTK(LJ(ID),LJ(IB),1) == 0) RETURN
               ENDIF
               CALL NONRELAT32 (ID, IB, IA, IB, IC, ID)
            ENDIF
         ELSE
            IF (IIRE == 0) RETURN
            CALL NONRELAT33 (ID, IB, IA, 1, IA, IB, IC, ID)
         ENDIF
      ELSE IF (IA == ID) THEN
         IF (IB==IA .OR. IC==IA) THEN
            IF (IB == IC) GO TO 10
            IF (IC == ID) THEN
               CALL NONRELAT31 (IC, IB, IA, IB, IC, ID, IIRE)
            ELSE
               IF (IIRE == 0) RETURN
               IF (ISOTOP == 1) THEN
                  IF (ITTK(LJ(IB),LJ(IB),1) == 0) RETURN
                  IF (ITTK(LJ(IC),LJ(IB),1) == 0) RETURN
               ENDIF
               CALL NONRELAT32 (IC, IB, IA, IB, IC, ID)
            ENDIF
         ELSE
            IF (IIRE == 0) RETURN
            CALL NONRELAT33 (IC, IB, IA, 2, IA, IB, ID, IC)
         ENDIF
      ELSE IF (IB == IC) THEN
         IF (IA==IB .OR. ID==IB) THEN
            IF (IA == ID) GO TO 10
            IF (ID == IB) THEN
               CALL NONRELAT31 (ID, IA, IA, IB, IC, ID, IIRE)
            ELSE
               IF (IIRE == 0) RETURN
               IF (ISOTOP == 1) THEN
                  IF (ITTK(LJ(IA),LJ(IA),1) == 0) RETURN
                  IF (ITTK(LJ(ID),LJ(IA),1) == 0) RETURN
               ENDIF
               CALL NONRELAT32 (ID, IA, IA, IB, IC, ID)
            ENDIF
         ELSE
            IF (IIRE == 0) RETURN
            CALL NONRELAT33 (ID, IA, IB, 2, IA, IB, ID, IC)
         ENDIF
      ENDIF
      RETURN
   10 CONTINUE
      WRITE (6, '(A)') ' ERRO IN NONRELAT3'
      STOP
      END SUBROUTINE NONRELAT3


text:
      SUBROUTINE DCHUR1(UPLO,TRANS,N,M,R,LDR,X,Z,LDZ,Y,RHO,C,S,INFO)
*     .. Scalar Arguments ..
      CHARACTER UPLO,TRANS
      INTEGER N,M,LDR,LDZ,INFO
*     .. Array Arguments ..
      DOUBLE PRECISION R(LDR,*)
      DOUBLE PRECISION X(*)
      DOUBLE PRECISION Z(LDZ,*)
      DOUBLE PRECISION Y(*)
      DOUBLE PRECISION RHO(*)
      DOUBLE PRECISION C(*)
      DOUBLE PRECISION S(*)
*
*  Purpose
*  =======
*
*  DCHUR1 performs a rank 1 update of a Cholesky factor R and
*  optionnally, of associated vectors Z
*
*  If UPLO = 'U', this solves the problems
*      R2'*R2 = R1'*R1 - x*x'
*      R2'*op(Z2) = R1'*op(Z1) - x*y'
*  Or if UPLO='L', the transposed versions
*      R2*R2'=R1*R1' - x*x'
*      R2*op(Z2)=R1*op(Z1) - x*y'
*  where R1 is R on entry, R2 is R on exit, ' means transpose,
*  Z1 is Z on entry, Z2 is Z on exit, op is transpose or identity.
*
*  The problem can also be written if UPLO='U' as finding U such that
*
*                              (R  op(Z))     (R2  op(Z2))
*                         U  * (        )  =  (          ) ,
*                              (X      Y)     ( 0     Y2 )
*
*  Or its transposed form if UPLO='L'
*
*                         (R      X)          (R2      0 )
*                         (        ) * U'  =  (          ) ,
*                         (op(Z)' Y)          (op(Z2)' Y2)
*
*  Where U is a product of Givens plane rotations U=U(N)*...*U(1)
*  And U(I) is a rotation in the plane (I,N+1)
*
*                        (     c(i)      s(i) )
*                  U(I)= (                    ) .
*                        (    -s(i)      c(i) )
*
*  When op(Z) are the factors of a least square problems,
*  and RHO is the norm of residuals of least square solution,
*  this is equivalent to adding a new observation (X,Y),
*  and update the factored solution R,Z,RHO of the problem.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  R is a upper triangular factor of R'*R
*          = 'L':  R is a lower triangular factor of R*R'.
*
*  TRANS    (input) CHARACTER*1
*          = 'N':  Z is a N-BY-M Matrix and op=identity.
*          = 'T':  Z is a M-BY-N Matrix and op=transpose
*
*  N       (input) INTEGER
*          The order of the matrix R.  N >= 0.
*
*  M       (input) INTEGER
*          The number of columns the matrix Z.  M >= 0.
*          In the case M = 0, the array Z,Y,RHO are not referenced
*
*  R       (input/output) DOUBLE PRECISION array, dimension (LDR,N)
*          On entry, the triangular matrix R.  If UPLO = 'U', the leading
*          N-by-N upper triangular part of R contains the upper
*          triangular part of the matrix R, and the strictly lower
*          triangular part of R is not referenced.  If UPLO = 'L', the
*          leading N-by-N lower triangular part of R contains the lower
*          triangular part of the matrix R, and the strictly upper
*          triangular part of R is not referenced.
*
*          On exit, if INFO = 0, R contains the update Cholesky factor
*
*  LDR     (input) INTEGER
*          The leading dimension of the array R, LDR >= max(1,N).
*
*  X       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, X contains the vector X used to update R.
*          On exit, X is destructed.
*
*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,k)
*          On entry, Z contains a list of M vectors of dimension N
*          to be updated with R.
*          If TRANS = 'N', Z is considered to be N-by-M (column vectors)
*          If TRANS = 'T', Z is considered to be M-by-N (row vectors)
*          On exit, if INFO = 0, Z contains the updated matrix Z
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.
*          If TRANS = 'N' LDZ >= max(1,N).
*          If TRANS = 'T' LDZ >= max(1,M).
*
*  Y       (input/output) DOUBLE PRECISION array, dimension (M)
*          On entry, Y contains the scalar quantities used to update Z.
*          on exit, Y contains the updated values Y2.
*
*  RHO     (input/output) DOUBLE PRECISION array, dimension (M)
*          on entry, RHO contains the norm of residual vectors.
*          on exit, RHO contains the updated norm of residual vectors.
*
*  C       (output) DOUBLE PRECISION array, dimension (N)
*          C contains the cosines of performed Givens rotations.
*
*  S       (output) DOUBLE PRECISION array, dimension (N)
*          S contains the sines of performed Givens rotations.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, the ith residual RHO was negative.
*
*
*  Note:
*  This is a variant of LINPACKs' dchud [1] & [2]
*  but with the following differences taking [3] into account
*  - uses BLAS drot and drotg for Given's rotations
*  - forces positive elements on the diagonal of R
*  - accepts upper or lower Cholesky factor R
*
*  Author: Nicolas Cellier
*
*  References:
*  [1] Linpack dchud http://www.netlib.org/linpack/dchud.f
*  [2] LINPACK User Guide http://www.uploading.com/files/F89HJTGH/linguide.zip.html
*  [3] Low Rank Updates for the Cholesky Decomposition - Matthias Seeger
*      http://infoscience.epfl.ch/record/161468/files/cholupdate.pdf
*
*     .. Local Scalars ..
      LOGICAL UPPER,DOTRAN
      INTEGER I
      DOUBLE PRECISION SC,W
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
*     .. External Subroutines ..
      EXTERNAL DROT,DROTG
*     .. Intrinsic Functions ..
      INTRINSIC DSQRT
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO=0
      UPPER=LSAME(UPLO,'U')
      DOTRAN=LSAME(TRANS,'T')
      IF( .NOT. UPPER .AND. .NOT. LSAME(UPLO,'L') ) THEN
        INFO=-1
      ELSE IF( .NOT.DOTRAN .AND. .NOT. LSAME(TRANS,'N') ) THEN
        INFO=-2
      ELSE IF( N.LT.0 ) THEN
        INFO=-3
      ELSE IF( M.LT.0 ) THEN
        INFO=-4
      ELSE IF( LDR.LT.N ) THEN
        INFO=-6
      ELSE IF (M.GT.0) THEN
        IF(DOTRAN) THEN
          IF( LDZ.LT.M ) INFO=-9
        ELSE
          IF( LDZ.LT.N ) INFO=-9
        ENDIF
      ENDIF
      IF(INFO.NE.0) RETURN
*
*     Quick return if possible
*
      IF(N.EQ.0) RETURN
*
* ------ FIRST PART: UPDATE R ------
*
      IF(UPPER) THEN
*       Compute givens rotations to annihilate x
        DO 10 I=1,N
          CALL DROTG(R(I,I), X(I), C(I), S(I))
*         Force positive diagonal values
          IF(R(I,I).LT.0.0D0) THEN
            R(I,I)=-R(I,I)
            C(I)=-C(I)
            S(I)=-S(I)
          ENDIF
*         Apply rotation on each row
          IF(I.LT.N) THEN
            CALL DROT(N-I,R(I,I+1),LDR,X(I+1),1,C(I),S(I))
          ENDIF
 10     CONTINUE
      ELSE
*       Compute givens rotations to annihilate x
        DO 20 I=1,N
          CALL DROTG(R(I,I), X(I), C(I), S(I))
*         Force positive diagonal values
          IF(R(I,I).LT.0.0D0) THEN
            R(I,I)=-R(I,I)
            C(I)=-C(I)
            S(I)=-S(I)
          ENDIF
*         Apply rotation on each column
          IF(I.LT.N) THEN
            CALL DROT(N-I,R(I+1,I),1,X(I+1),1,C(I),S(I))
          ENDIF
 20     CONTINUE
      ENDIF
*
* ------  SECOND PART: UPDATE Z AND RHO IF REQUESTED ------
*
      IF(M.GT.0) THEN
        IF(DOTRAN) THEN
*         Update colums of z
          DO 30 I=1,N
            CALL DROT(M,Z(1,I),1,Y(1),1,C(I),S(I));
 30       CONTINUE
        ELSE
*         Update rows of z
          DO 40 I=1,N
            CALL DROT(M,Z(I,1),LDZ,Y(1),1,C(I),S(I));
 40       CONTINUE
        ENDIF
*       Update rho
        DO 50 I=1,M
          W=DABS(Y(I))
          IF(RHO(I).GT.0.0D0) THEN
            SC=RHO(I)+W
            RHO(I)=SC*DSQRT((RHO(I)/SC)**2+(W/SC)**2)
          ELSE
            INFO=I
            GOTO 100
          ENDIF
 50     CONTINUE
      ENDIF
 100  CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: nick-terry/Splitting-GP

max_issues_repo_name: nick-terry/Splitting-GP

max_forks_repo_name: nick-terry/Splitting-GP


text:
! -*- Mode: Fortran; -*-
!
!  (C) 2001 by Argonne National Laboratory.
!      See COPYRIGHT in top-level directory.
!
      program main
      implicit none

      include 'mpif.h'

!     Fortran equivalent of perf.c

      integer SIZE
      parameter (SIZE=1048576*4)
!     read/write size per node in bytes

      integer buf(SIZE/4), j, mynod, nprocs, ntimes, flag, foo
      double precision stim, read_tim, write_tim, new_read_tim
      double precision new_write_tim, min_read_tim, min_write_tim
      double precision read_bw, write_bw
      integer fh, status(MPI_STATUS_SIZE), ierr, argc, iargc, i
      character*1024 str    ! used to store the filename
      integer*8 offset

      ntimes = 5
      min_read_tim = 10000000.0D0
      min_write_tim = 10000000.0D0

      call MPI_INIT(ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, foo, ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, mynod, ierr)

!     process 0 takes the file name as a command-line argument and
!     broadcasts it to other processes

      if (mynod .eq. 0) then
         argc = iargc()
         i = 0
         call getarg(i,str)
         do while ((i .lt. argc) .and. (str .ne. '-fname'))
            i = i + 1
            call getarg(i,str)
         end do
         if (i .ge. argc) then
            print *
            print *, '*#  Usage: fperf -fname filename'
            print *
            call MPI_ABORT(MPI_COMM_WORLD, 1, ierr)
         end if

         i = i + 1
         call getarg(i,str)
         call MPI_BCAST(str, 1024, MPI_CHARACTER, 0,                    &
     &        MPI_COMM_WORLD, ierr)
         print *, 'Access size per process = ', SIZE, ' bytes',         &
     &        ', ntimes = ', ntimes
      else
         call MPI_BCAST(str, 1024, MPI_CHARACTER, 0,                    &
     &        MPI_COMM_WORLD, ierr)
      end if

      offset = mynod*SIZE
      do j=1, ntimes
         call MPI_FILE_OPEN(MPI_COMM_WORLD, str,                        &
     &        MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

         call MPI_FILE_SEEK(fh, offset, MPI_SEEK_SET, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         stim = MPI_WTIME()
         call MPI_FILE_WRITE(fh, buf, SIZE, MPI_BYTE, status, ierr)
         write_tim = MPI_WTIME() - stim

         call MPI_FILE_CLOSE(fh, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)

         call MPI_FILE_OPEN(MPI_COMM_WORLD, str,                        &
     &        MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

         call MPI_FILE_SEEK(fh, offset, MPI_SEEK_SET, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         stim = MPI_WTIME()
         call MPI_FILE_READ(fh, buf, SIZE, MPI_BYTE, status, ierr)
         read_tim = MPI_WTIME() - stim

         call MPI_FILE_CLOSE(fh, ierr)

         call MPI_ALLREDUCE(write_tim, new_write_tim, 1,                &
     &        MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE(read_tim, new_read_tim, 1,                  &
     &        MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)

         if (new_read_tim .lt. min_read_tim) then
            min_read_tim = new_read_tim
         end if
         if (new_write_tim .lt. min_write_tim) then
            min_write_tim = new_write_tim
         end if
      end do

      if (mynod .eq. 0) then
         read_bw = (SIZE*nprocs*1.0D0)/(min_read_tim*1000000.0D0)
         write_bw = (SIZE*nprocs*1.0D0)/(min_write_tim*1000000.0D0)
         print *, 'Write bandwidth without file sync = ',               &
     &        write_bw, ' Mbytes/sec'
         print *, 'Read bandwidth without prior file sync = ',          &
     &        read_bw, ' Mbytes/sec'
      end if

      min_read_tim = 10000000.0D0
      min_write_tim = 10000000.0D0

      flag = 0
      do j=1, ntimes
         call MPI_FILE_OPEN(MPI_COMM_WORLD, str,                        &
     &        MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

         call MPI_FILE_SEEK(fh, offset, MPI_SEEK_SET, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         stim = MPI_WTIME()
         call MPI_FILE_WRITE(fh, buf, SIZE, MPI_BYTE, status, ierr)
         call MPI_FILE_SYNC(fh, ierr)
         write_tim = MPI_WTIME() - stim
         if (ierr .eq. MPI_ERR_UNKNOWN) then
            flag = 1
         end if

         call MPI_FILE_CLOSE(fh, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)

         call MPI_FILE_OPEN(MPI_COMM_WORLD, str,                        &
     &        MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

         call MPI_FILE_SEEK(fh, offset, MPI_SEEK_SET, ierr)

         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         stim = MPI_WTIME()
         call MPI_FILE_READ(fh, buf, SIZE, MPI_BYTE, status, ierr)
         read_tim = MPI_WTIME() - stim

         call MPI_FILE_CLOSE(fh, ierr)

         call MPI_ALLREDUCE(write_tim, new_write_tim, 1,                &
     &        MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE(read_tim, new_read_tim, 1,                  &
     &        MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)

         if (new_read_tim .lt. min_read_tim) then
            min_read_tim = new_read_tim
         end if
         if (new_write_tim .lt. min_write_tim) then
            min_write_tim = new_write_tim
         end if

      end do

      if (mynod .eq. 0) then
         if (flag .eq. 1) then
            print *, 'MPI_FILE_SYNC returns error.'
         else
            read_bw = (SIZE*nprocs*1.0D0)/(min_read_tim*1000000.0D0)
            write_bw = (SIZE*nprocs*1.0D0)/(min_write_tim*1000000.0D0)
            print *, 'Write bandwidth including file sync = ',          &
     &           write_bw, ' Mbytes/sec'
            print *, 'Read bandwidth after file sync = ',               &
     &           read_bw, ' Mbytes/sec'
         end if
      end if

      call MPI_FINALIZE(ierr)

      stop
      end


text:
program NUDFT
implicit none

    integer, parameter :: general_kind=4

    real(kind=general_kind), parameter :: pi=atan(1.)*4.0

	real(kind=general_kind), parameter :: T_low=0.1,T_hi=10.0,dT=0.0001
    real(kind=general_kind), allocatable :: dat(:,:),x(:),y(:),T(:)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    call import_OGLE_data(5,dat)

    allocate(x(size(dat(:,1))),y(size(dat(:,1))))
	x=dat(:,1)
	y=dat(:,2)
    !normalize101
    y=(y-mean(y))/std(y)

    print*,NFT(5.612)

contains

    real(kind=general_kind) function mean(arr)
        implicit none
        real(kind=general_kind), intent(in), dimension(:) :: arr

        mean=sum(arr)/size(arr)
    end function mean

    real(kind=general_kind) function std(arr)
        implicit none
        real(kind=general_kind), intent(in), dimension(:) :: arr

        std=sqrt(sum((arr-mean(arr))**2)/(size(arr)-1))
    end function std

    subroutine print_array(arr,arr_size)
        implicit none
        real(kind=general_kind), intent(in) :: arr(*)
        integer, intent(in) :: arr_size
        integer :: i
        do i = 1,arr_size
            print*,arr(i)
        end do
    end subroutine print_array

    integer function count_file_lines(filename)
        implicit none
        character (len=21), intent(in) :: filename

        count_file_lines=0
        open(50,file=filename)
        do
            read(50,*,end=10)
            count_file_lines=count_file_lines+1
        end do
        10 close(50)
    end function count_file_lines

    character(len=4) function str4(k)
        implicit none
        ! left-pad an integer with 4 zeros
        integer, intent(in) :: k
        write (str4, '(I4.4)') k
    end function str4

    character(len=21) function OGLE_file_string(k)
        implicit none
        ! ogle filename generator
        integer, intent(in) :: k
        OGLE_file_string="OGLE-LMC-CEP-"//str4(k)//".dat"
    end function OGLE_file_string

    subroutine import_OGLE_data(k,M)
        implicit none
        integer, intent(in) :: k
		character (len=21) :: filename
		real(kind=general_kind), allocatable :: M(:,:)
		integer i,row,col

        filename=OGLE_file_string(k)

		row=count_file_lines(filename)
		col=3 ! specific for OGLE files
		allocate(M(row,col))

		open(50,file=filename)
		do i=1,row
			read(50,*) M(i,:)
		end do
		close(50)
	end subroutine import_OGLE_data

    complex(kind=general_kind) function NFT(T_)
        implicit none
        ! complex nonuniform fourier transform
        real(kind=general_kind), intent(in) :: T_
        NFT=sum(y*exp(cmplx(0.0,-2*pi)*x/T_))
    end function NFT

end program NUDFT

max_stars_repo_name: santiagohenao/NUDFT-III-Fortran

max_issues_repo_name: santiagohenao/NUDFT-III-Fortran

max_forks_repo_name: santiagohenao/NUDFT-III-Fortran


text:
subroutine u_wmap(lundat    ,header    ,runid     ,itime     ,idate     , &
                & timnow    ,mmax      ,nmax      ,kmax      ,nmaxus    , &
                & nsrc      ,ltur      ,lmaxd     ,lstsci    ,lsal      , &
                & ltem      ,ilay      ,zmodel    ,kfu       ,kfv       , &
                & kcs       ,kfs       ,kfumin    ,kfvmin    ,kfumax    , &
                & kfvmax    ,kfsmin    ,kfsmax    ,mnksrc    ,xz        , &
                & yz        ,dps       ,rho       ,s1        ,u1        , &
                & v1        ,w1        ,wphy      ,r1        ,p1        , &
                & rtur1     ,taubpu    ,taubpv    ,taubsu    ,taubsv    , &
                & vicww     ,dicww     ,vicuv     ,rich      ,umnldf    , &
                & vmnldf    )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program 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 version 3.
!
!  This program 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/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: u_wmap.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/plugin_user/src/u_wmap.f90 $
!!--description-----------------------------------------------------------------
!
! Example routine to write vertical eddy viscosity at
! layer=ILAY, water elevation + height at all waterlevel
! points (map) to a user defined output file
!
! Eddy viscosity are averaged to centre of the cell
! Values of KCx: 0= inactive;      1= active; 2= open boundary
! Values of KFx: 0= temporary dry; 1= wet
! WARNING: RBUFF0 and RBUFF1 uses the same memory location
! see call in postpr
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
   !
    implicit none
!
! Global variables
!
    integer                                         , intent(in) :: idate
    integer                                         , intent(in) :: ilay
    integer                                         , intent(in) :: itime
    integer                                         , intent(in) :: kmax !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: lmaxd !  Description and declaration in dimens.igs
    integer                                         , intent(in) :: lsal !  Description and declaration in dimens.igs
    integer                                         , intent(in) :: lstsci !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: ltem !  Description and declaration in dimens.igs
    integer                                         , intent(in) :: ltur !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: lundat
    integer                                         , intent(in) :: mmax !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: nmax !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: nmaxus !  Description and declaration in esm_alloc_int.f90
    integer                                         , intent(in) :: nsrc !  Description and declaration in esm_alloc_int.f90
    integer, dimension(7, nsrc)                     , intent(in) :: mnksrc !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kcs !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfs !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfsmax !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfsmin !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfu !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfumax !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfumin !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfv !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfvmax !  Description and declaration in esm_alloc_int.f90
    integer, dimension(nmax, -1:mmax + 2)           , intent(in) :: kfvmin !  Description and declaration in esm_alloc_int.f90
    logical                                         , intent(in) :: zmodel !  Description and declaration in procs.igs
    real(fp)                                        , intent(in) :: timnow
    real(prec), dimension(nmax, -1:mmax + 2)              , intent(in) :: dps !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: s1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: taubpu !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: taubpv !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: taubsu !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: taubsv !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: umnldf !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: vmnldf !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: xz !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2)              , intent(in) :: yz !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, 0:kmax)      , intent(in) :: dicww !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, 0:kmax)      , intent(in) :: rich !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, 0:kmax)      , intent(in) :: vicww !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, 0:kmax)      , intent(in) :: w1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, 0:kmax, ltur), intent(in) :: rtur1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax + 2)    , intent(in) :: vicuv !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax)        , intent(in) :: p1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax)        , intent(in) :: rho !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax)        , intent(in) :: u1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax)        , intent(in) :: v1 !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax)        , intent(in) :: wphy !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmax, -1:mmax + 2, kmax, lstsci), intent(in) :: r1 !  Description and declaration in esm_alloc_real.f90
    character(*)                                    , intent(in) :: runid
    character(131), dimension(10)                                :: header !  Description and declaration in postpr.igs
!
! Local variables
!
    integer :: i
    integer :: m
    integer :: n
    logical :: lvar
    real(fp):: vw
    real(fp):: wh
    real(fp):: xydef
    real(fp):: zdef
!
!! executable statements -------------------------------------------------------
!
    !
    ! Initialisation
    !
    xydef = 999999.
    zdef = 999.999
    !
    ! Header;
    !
    header(1) = '* MAP output of run:' // runid
    header(2) = '* Simulation results at simulation time: '
    header(3) = '* column    1 : x-coordinates'
    header(4) = '* column    2 : y-coordinates'
    header(5) = '* column    3 : water level (m)'
    header(6) = '* column    4 : water height (m)'
    header(7) = '* column    5 : vert. eddy visc. (m2/s) at layer: '
    write (header(2)(42:), '(i8,1x,i6)') idate, itime
    write (header(7)(51:), '(i5)') ilay
    !
    do i = 1, 7
       write (lundat, '(  a)') header(i)
    enddo
    write (lundat, '(4i6)') mmax*nmaxus, 5, mmax, nmaxus
    !
    ! determine computed values
    ! if kcs = 1 and kfs = 1 then active an wet point; else point
    ! inactive
    !
    do n = 1, nmaxus
       do m = 1, mmax
          wh = s1(n, m) + real(dps(n, m),fp)
          lvar = ((kcs(n, m)==1) .and. (kfs(n, m)==1))
          vw = zdef
          if (ilay>=1 .and. ilay<=kmax .and. lvar) then
             vw = (vicww(n, m, ilay) + vicww(n, m, ilay - 1))/2.
          endif
          !
          ! Write values to output file
          !
          if (lvar) then
             write (lundat, '(2(e13.6,1x),5(1x,g13.6) )') &
                 & xz(n, m), yz(n, m), wh, s1(n, m), vw
          else
             write (lundat, '(2(e13.6,1x),5(1x,g13.6) )') &
                 & (xydef, i = 1, 2), (zdef, i = 1, 3)
          endif
       enddo
    enddo
end subroutine u_wmap

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


FUNCTION real_to_char ( value )
!
!  Purpose:
!    To convert a real value into a 12-character string, with the
!    number printed in as readable a format as possible considering
!    its range.  This routine prints out the number according to the
!    following rules:
!       1.  value > 9999999.                ES12.5
!       2.  value < -999999.                ES12.5
!       3.  0.    <  ABS(value) < 0.01      ES12.5
!       4.  value = 0.0                     F12.4
!       5.  Otherwise                       F12.4
!
!  Record of revisions:
!      Date       Programmer          Description of change
!      ====       ==========          =====================
!    11/25/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Data dictionary: declare calling parameter types & definitions
REAL, INTENT(IN) :: value           ! value to convert to char form
CHARACTER (len=12) :: real_to_char  ! Output character string

! Data dictionary: declare local variable types & definitions
CHARACTER(len=9) :: fmt             ! Format descriptor
CHARACTER(len=12) :: string         ! Output string

! Clear string before use
string = ' '

! Select proper format
IF ( value > 9999999. ) THEN
   fmt = '(ES12.5)'
ELSE IF ( value < -999999. ) THEN
   fmt = '(ES12.5)'
ELSE IF ( value == 0. ) THEN
   fmt = '(F12.4)'
ELSE IF ( ABS(value) < 0.01 ) THEN
   fmt = '(ES12.5)'
ELSE
   fmt = '(F12.4)'
END IF

! Convert to character form.
WRITE (string,fmt) value
real_to_char = string

END FUNCTION real_to_char


text:
      function temp (x, T)
      use stel_kinds
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      real(rprec) :: x, temp, T(0:10)
C-----------------------------------------------
      temp = (T(0) + x*(T(1) + x*(T(2) + x*(T(3) + x*(T(4) +
     1     x*(T(5) + x*(T(6) + x*(T(7) + x*(T(8) + x*(T(9) + x*T(10)
     2   ))))))))))
      end function temp

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      SUBROUTINE zgtdts (IFLTAB, CA, NA, CB, NB, CC, NC, CF, NF,
     * JULS, ISTIME, JULE, IETIME, INTL, IHEAD, NHEAD, IOFSET,
     * VALUES, NVALS, CUNITS, CTYPE, ISTAT)
C
C     Replaced by zrrts6
C
      INTEGER IFLTAB(*)
      INTEGER IHEAD(*)
      REAL VALUES(*)
      CHARACTER CA*(*), CB*(*), CC*(*), CF*(*), CUNITS*(*), CTYPE*(*)
      CHARACTER CPATH*392, CE*12, CDATE*9, CTIME*4
      LOGICAL LQREAD
C
C
      CE = ' '
      I = 2
      CALL zgintl6 (INTL, CE, N, I)
      IF  (I.NE.0) GO TO 900
C
      ND = 0
      NE = 12
      CPATH = ' '
      CALL zfpn (CA, NA, CB, NB, CC, NC, ' ', ND, CE, NE,
     * CF, NF, CPATH, NPATH)
C
      CALL DATCLN (JULS, ISTIME, JUL, IT)
      CALL JULDAT (JUL, 104, CDATE, N)
      N = M2IHM (IT, CTIME)
C
      N = NOPERS (INTL,0,JULS,ISTIME,JULE,IETIME) + 1
      IF (N.GT.NVALS) N = NVALS
      NVALS = N
      IF (NVALS.LE.0) GO TO 910
C
C
      CALL zrrtsx ( IFLTAB, CPATH(1:NPATH), CDATE, CTIME, NVALS,
     * VALUES, IQUAL, .FALSE., LQREAD, CUNITS, CTYPE, IUHEAD, 0,
     * JHEAD, IOFSET, JCOMP, ISTAT)
C
C     Fix up header
      !IF ((ISTAT.LE.4).AND.(NHEAD.GE.20)) THEN
      !IHEAD(20) = IOFSET
      !ENDIF
C
C
 800  CONTINUE
      RETURN
C
 900  CONTINUE
      !WRITE (MUNIT, 901) INTL
 901  FORMAT (/,' -----DSS---zgtdts6;  ERROR:  Unrecognizable Time',
     * ' Interval',/' Interval Specified:',I6,/)
      ISTAT = 12
      GO TO 800
C
 910  CONTINUE
      !WRITE (MUNIT, 911) NVALS
 911  FORMAT (/,' -----DSS---zgtdts6;  ERROR:  Illegal Time Window',
     * ' Given',/' Number of Values:',I6,/)
      !WRITE (MUNIT, 912) CDATE, CTIME, JULS
 912  FORMAT (' Starting Date: ',A,',  Time: ',A,I8)
      !CALL DATCLN (JULE, IETIME, JUL, IT)
      !CALL JULDAT (JUL, 104, CDATE, N)
      N = M2IHM (IT, CTIME)
      !WRITE (MUNIT, 913) CDATE, CTIME, JULE
 913  FORMAT (' Ending Date: ',A,',  Time: ',A,I8,/)
      ISTAT = 15
      GO TO 800
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
./verilog/tb_dltdc_model.sv
./verilog/dltdc_model.sv
./verilog/functions.v

max_stars_repo_name: idea-fasoc/fasoc

max_issues_repo_name: idea-fasoc/fasoc

max_forks_repo_name: idea-fasoc/fasoc


text:
module implicit_mod
!||zn, pjk 5/13/99
!following vb's changes to old implicit_mod.
!implicit_init works with global fields, everything else on domain

use              fms_mod, only: error_mesg, FATAL, write_version_number

use        constants_mod, only: rdgas, radius, cp_air

use press_and_geopot_mod, only: pressure_variables

use    matrix_invert_mod, only: invert

use       transforms_mod, only: get_spec_domain

implicit none

private

public :: implicit_init, implicit_correction, implicit_end

interface linear_geopotential
   module procedure   linear_geopotential_1d, linear_geopotential_3d
end interface

interface linear_tp_tendency
   module procedure   linear_tp_tendency_1d, linear_tp_tendency_3d
end interface

logical :: module_is_initialized = .false.

character(len=128), parameter :: version = '$Id: implicit.f90,v 13.0 2006/03/28 21:17:54 fms Exp $'
character(len=128), parameter :: tagname = '$Name: latest $'

real,    allocatable, dimension(:)   :: ref_temperature_implicit
real,    allocatable, dimension(:)   :: ref_ln_p_half, ref_ln_p_full, del_ln_p_half, del_ln_p_full
real,    allocatable, dimension(:,:) :: eigen
integer, allocatable, dimension(:,:) :: wavenumber

real :: alpha, ref_surf_p_implicit

real, allocatable, dimension(:,:)   :: div_mat
real, allocatable, dimension(:)     :: h
real, allocatable, dimension(:,:,:) :: wave_matrix
real, allocatable, dimension(:)     :: pk, bk, dpk, dbk

real :: dt = 0.0
real :: xi

integer :: num_levels, num_total_wavenumbers
character(len=64) :: vert_difference_option
integer :: ms, me, ns, ne

contains

!------------------------------------------------------------------------

subroutine implicit_init(pk_in, bk_in, ref_temperature_implicit_in, ref_surf_p_implicit_in, &
                         num_total_wavenumbers_in, eigen_in, wavenumber_in, alpha_in,       &
                         vert_difference_option_in)

real,    intent(in), dimension(:)     :: pk_in, bk_in, ref_temperature_implicit_in
real,    intent(in), dimension(0:,0:) :: eigen_in
integer, intent(in), dimension(0:,0:) :: wavenumber_in

real,    intent(in) :: ref_surf_p_implicit_in, alpha_in
integer, intent(in) :: num_total_wavenumbers_in
character(len=*), intent(in) :: vert_difference_option_in

real, dimension(size(ref_temperature_implicit_in,1)  ) :: p_full_work_1, p_full_work_2, ln_p_full_work_1, ln_p_full_work_2
real, dimension(size(ref_temperature_implicit_in,1)+1) :: p_half_work_1, p_half_work_2, ln_p_half_work_1, ln_p_half_work_2

real :: surface_p_1, eps=1.e-5
integer :: k

if(module_is_initialized) return

call write_version_number(version, tagname)

call get_spec_domain(ms, me, ns, ne)

num_levels = size(ref_temperature_implicit_in,1)
num_total_wavenumbers = num_total_wavenumbers_in

allocate(ref_temperature_implicit(num_levels))
allocate(ref_ln_p_half(num_levels+1))
allocate(ref_ln_p_full(num_levels))
allocate(del_ln_p_half(num_levels+1))
allocate(del_ln_p_full(num_levels))
allocate(div_mat(num_levels,num_levels))
allocate(wave_matrix(num_levels,num_levels,0:num_total_wavenumbers))
allocate(h(num_levels))
allocate(eigen(0:size(eigen_in,1)-1,0:size(eigen_in,2)-1))
allocate(wavenumber(0:size(eigen_in,1)-1,0:size(eigen_in,2)-1))

ref_temperature_implicit = ref_temperature_implicit_in
ref_surf_p_implicit = ref_surf_p_implicit_in
eigen         = eigen_in
wavenumber    = wavenumber_in
alpha         = alpha_in
vert_difference_option = vert_difference_option_in

allocate( pk(size(pk_in(:))),  bk(size(bk_in,1)))
allocate(dpk(num_levels  ), dbk(num_levels  ))

pk = pk_in
bk = bk_in

do k = 1,num_levels
  dpk(k) = pk(k+1) - pk(k)
  dbk(k) = bk(k+1) - bk(k)
end do

call pressure_variables(p_half_work_1, ref_ln_p_half, p_full_work_1, ref_ln_p_full, ref_surf_p_implicit)

!del_ln_p_half = dln_p_half_dps(ref_surf_p_implicit)
!del_ln_p_full = dln_p_full_dps(ref_surf_p_implicit)

!  functions dln_p_half_dps and del_ln_p_full are used nowhere else and
!  they perform very simple calculations, so it only makes the code more
!  confusing to do these calculations in functions.

do k=2,num_levels+1
  del_ln_p_half(k) = bk(k)/(pk(k)+bk(k)*ref_surf_p_implicit)
enddo

if(pk(1).eq.0.0) then
  del_ln_p_half(1) = 1.0/ref_surf_p_implicit
else
  del_ln_p_half(1) = bk(1)/(pk(1)+bk(1)*ref_surf_p_implicit)
endif

surface_p_1 = ref_surf_p_implicit*(1.0 - 0.5*eps)
call pressure_variables(p_half_work_1, ln_p_half_work_1, p_full_work_1, ln_p_full_work_1, surface_p_1)

surface_p_1 = ref_surf_p_implicit*(1.0 + 0.5*eps)
call pressure_variables(p_half_work_2, ln_p_half_work_2, p_full_work_2, ln_p_full_work_2, surface_p_1)

del_ln_p_full = (ln_p_full_work_2 - ln_p_full_work_1)/(eps*ref_surf_p_implicit)

call build_matrix

module_is_initialized = .true.

return
end subroutine implicit_init

!------------------------------------------------------------------------

subroutine build_matrix

real, dimension(num_levels) :: input, zero, dt_t
real, dimension(num_levels+1) :: zero_1
real :: dt_p

real, dimension(num_levels, num_levels) :: gamma, tau
real, dimension(num_levels) :: nu, h1, h2

integer :: k, kk, kkk

zero  = 0.0
zero_1 = 0.0

tau = 0.0
nu = 0.0
do k=1,num_levels
  input    = 0.0
  input(k) = 1.0

  call linear_tp_tendency(input, ref_temperature_implicit, ref_surf_p_implicit, &
                          ref_ln_p_half, ref_ln_p_full, dt_p, dt_t)
  nu(k)    = - dt_p
  tau(:,k) = - dt_t

  gamma(:,k) = linear_geopotential(input,zero_1,zero, &
                    ref_temperature_implicit, ref_ln_p_half, ref_ln_p_full)
end do

h1 = pres_grad_funct(ref_temperature_implicit, ref_ln_p_half, ref_ln_p_full, ref_surf_p_implicit)

h2 = linear_geopotential(zero,del_ln_p_half,del_ln_p_full, &
                    ref_temperature_implicit, ref_ln_p_half, ref_ln_p_full)

h = h1 + h2

do k = 1,num_levels
  do kk = 1,num_levels
    div_mat(k,kk)= h(k)*nu(kk)
    do kkk = 1,num_levels
      div_mat(k,kk)=div_mat(k,kk)+gamma(k,kkk)*tau(kkk,kk)
    enddo
  enddo
enddo

return
end subroutine build_matrix

!-------------------------------------------------------------------------

subroutine build_wave_matrices

real :: factor, det
integer :: k, L

wave_matrix = 0.0
do L = 0,num_total_wavenumbers
  factor = xi*xi*L*(L+1)/radius**2
  do k = 1,num_levels
    wave_matrix(k,k,L) = 1.0
  end do
  wave_matrix(:,:,L) = wave_matrix(:,:,L) + factor*div_mat(:,:)
  call invert(wave_matrix(:,:,L), det)
end do

return
end subroutine build_wave_matrices

!---------------------------------------------------------------------

subroutine implicit_correction(dt_divs, dt_ts, dt_ln_ps, divs, ts, ln_ps, dt_in, previous, current)
!in parallel, wavenumber (0,0) actually means (ms,ns)
!it is necessary to use the (0:,0:) dimensioning here because of the m-n loop below
complex, intent(inout), dimension(0:,0:,:) :: dt_divs, dt_ts
complex, intent(inout), dimension(0:,0:) :: dt_ln_ps
complex, intent(in), dimension(0:,0:,:,:) :: divs, ts
complex, intent(in), dimension(0:,0:,:) :: ln_ps
real, intent(in) :: dt_in
integer,intent(in) :: previous, current

complex, dimension(0:size(divs,1)-1,0:size(divs,2)-1,size(divs,3)) :: dt_ts_temp
complex, dimension(num_levels) ::  work
complex, dimension(0:size(divs,1)-1,0:size(divs,2)-1) :: dt_ps_temp
integer :: m,n,L

if(.not.module_is_initialized) then
  call error_mesg('implicit_correction','failed to initialize implicit', FATAL)
end if

if(dt_in.ne.dt) then
   dt = dt_in
   xi = dt*alpha
   call build_wave_matrices
end if

call adjust_dt_divs(dt_divs, dt_ts, dt_ln_ps, divs, ts, ln_ps, previous, current)

do n = 0,size(divs,2)-1
  do m = 0,size(divs,1)-1
!||zn, pjk 5/13/99, following vb's dmsm/implicit.f90
    L = wavenumber(m+ms,n+ns)
    if(L.le.num_total_wavenumbers) then
      work = matmul(wave_matrix(:,:,L),dt_divs(m,n,:))
      dt_divs(m,n,:) = work
    end if
  end do
end do

call linear_tp_tendency(dt_divs, ref_temperature_implicit, ref_surf_p_implicit, &
                        ref_ln_p_half, ref_ln_p_full, dt_ps_temp, dt_ts_temp)

dt_ts    = dt_ts    + xi*dt_ts_temp
dt_ln_ps = dt_ln_ps + xi*dt_ps_temp/ref_surf_p_implicit

return
end subroutine implicit_correction

!---------------------------------------------------------------------
subroutine adjust_dt_divs(dt_divs, dt_ts, dt_ln_ps, divs, ts, ln_ps, previous, current)
  !||zn, pjk 5/13/99, following vb's dmsm/implicit.f90. (0:,0:) -> (:,:)
complex, intent(inout), dimension(:,:,:) :: dt_divs
complex, intent(inout), dimension(:,:,:) :: dt_ts
complex, intent(inout), dimension(:,:) :: dt_ln_ps
complex, intent(in), dimension(:,:,:,:) :: divs, ts
complex, intent(in), dimension(:,:,:) :: ln_ps
integer, intent(in) :: previous, current

complex, dimension(size(ts,1),size(ts,2),size(ts,3)) :: divs_temp, ts_temp, geopot, dt_ts_temp, ln_p_full_temp
complex, dimension(size(ts,1),size(ts,2),size(ts,3)+1) :: ln_p_half_temp
complex, dimension(size(ts,1),size(ts,2)) :: ps_temp, dt_ps_temp

integer :: k

divs_temp = divs(:,:,:,previous) - divs(:,:,:,current)

call linear_tp_tendency(divs_temp, ref_temperature_implicit, ref_surf_p_implicit, &
                        ref_ln_p_half, ref_ln_p_full, dt_ps_temp, dt_ts_temp)

dt_ts    = dt_ts    + dt_ts_temp
dt_ln_ps = dt_ln_ps + dt_ps_temp/ref_surf_p_implicit

ts_temp =    ts(:,:,:,previous) -    ts(:,:,:,current) + xi*dt_ts
ps_temp = ln_ps(:,:,  previous) - ln_ps(:,:,  current) + xi*dt_ln_ps

ln_p_full_temp = (0.,0.)
ln_p_half_temp = (0.,0.)

geopot = linear_geopotential(ts_temp,ln_p_half_temp,ln_p_full_temp, ref_temperature_implicit, ref_ln_p_half, ref_ln_p_full)

do k=1,num_levels
  dt_divs(:,:,k)=dt_divs(:,:,k)+eigen(ms:me,ns:ne)*(geopot(:,:,k)+h(k)*ps_temp(:,:)*ref_surf_p_implicit)
enddo

return
end subroutine adjust_dt_divs

!-----------------------------------------------------------------------

function linear_geopotential_3d(del_t, del_ln_p_half, del_ln_p_full, t, ln_p_half, ln_p_full) result (geopot)
!||zn, pjk 5/13/99, following vb's dmsm/implicit.f90. (0:,0:) -> (:,:)

complex, intent (in), dimension (:,:,:) :: del_t
complex, intent (in), dimension (:,:,:) :: del_ln_p_half
complex, intent (in), dimension (:,:,:) :: del_ln_p_full
real,    intent (in), dimension (:)       :: t
real,    intent (in), dimension (:)       :: ln_p_half
real,    intent (in), dimension (:)       :: ln_p_full

complex, dimension (size(del_t,1),size(del_t,2),size(del_t,3))   :: geopot
complex, dimension (size(del_t,1),size(del_t,2),size(del_t,3)+1) :: geopot_half

integer :: k

geopot_half(:,:,num_levels+1) = (0.,0.)

do k=num_levels,2,-1
  geopot_half(:,:,k) = geopot_half(:,:,k+1)                                &
    + rdgas*(del_t(:,:,k)*(ln_p_half(k+1)         - ln_p_half(k))    &
    +t(k)        *(del_ln_p_half(:,:,k+1) - del_ln_p_half(:,:,k)))
enddo

do k=1,num_levels
  geopot(:,:,k) = geopot_half(:,:,k+1)                                     &
    + rdgas*(del_t(:,:,k)*(ln_p_half(k+1)     - ln_p_full(k))        &
    +t(k)        *(del_ln_p_half(:,:,k+1) - del_ln_p_full(:,:,k)))
enddo

return
end function linear_geopotential_3d
!-----------------------------------------------------------------------

function linear_geopotential_1d(del_t,del_ln_p_half,del_ln_p_full, t, ln_p_half, ln_p_full) result (geopot)

real, intent (in), dimension (:) :: del_t
real, intent (in), dimension (:) :: del_ln_p_half
real, intent (in), dimension (:) :: del_ln_p_full
real, intent (in), dimension (:) :: t
real, intent (in), dimension (:) :: ln_p_half
real, intent (in), dimension (:) :: ln_p_full

real, dimension (size(del_t,1)) :: geopot

complex, dimension(0:0, 0:0, size(del_t,1)  ) :: del_t_3d, geopot_3d
complex, dimension(0:0, 0:0, size(del_t,1)+1) :: del_ln_p_half_3d
complex, dimension(0:0, 0:0, size(del_t,1)  ) :: del_ln_p_full_3d

del_t_3d(0,0,:)         = cmplx(del_t,         0.0)
del_ln_p_half_3d(0,0,:) = cmplx(del_ln_p_half, 0.0)
del_ln_p_full_3d(0,0,:) = cmplx(del_ln_p_full, 0.0)

geopot_3d = linear_geopotential(del_t_3d, del_ln_p_half_3d, del_ln_p_full_3d, t, ln_p_half, ln_p_full)

geopot = real(geopot_3d(0,0,:))

return
end function linear_geopotential_1d
!-----------------------------------------------------------------------

function pres_grad_funct(tg, ln_p_half, ln_p_full, p_surf) result(x)

real, intent (in), dimension (:) :: tg, ln_p_half, ln_p_full
real, intent (in)                :: p_surf
real, dimension (size(tg,1)) :: x

real :: dlog_1, dlog_2
integer :: k

if(trim(vert_difference_option) == 'simmons_and_burridge') then
  do k=1,size(ln_p_full,1)
    dlog_1 = ln_p_half(k+1) - ln_p_full(k)
    dlog_2 = ln_p_full(k)   - ln_p_half(k)
    x(k) = rdgas*tg(k) * (bk(k+1)*dlog_1 + bk(k)*dlog_2)/(dpk(k) + dbk(k)*p_surf)
  end do
else if(trim(vert_difference_option) == 'mcm') then
  do k=1,size(ln_p_full,1)
    x(k) = rdgas*tg(k)/p_surf
  end do
endif

return
end function pres_grad_funct
!-----------------------------------------------------------------------

subroutine linear_tp_tendency_3d (div, t_ref, p_surf_ref, ln_p_half_ref, ln_p_full_ref, dt_p_surf, dt_t)

complex, intent (in), dimension (0:,0:,:) :: div
real,    intent (in), dimension (:) :: t_ref, ln_p_half_ref, ln_p_full_ref
real,    intent (in)  :: p_surf_ref

complex, intent (out), dimension (0:,0:)   :: dt_p_surf
complex, intent (out), dimension (0:,0:,:) :: dt_t

real ::  dp, dp_inv, dlog_1, dlog_3, p_full_ref
complex , dimension(0:size(div,1)-1,0:size(div,2)-1) :: dmean, dmean_tot
complex , dimension(0:size(div,1)-1,0:size(div,2)-1,size(div,3)+1) :: vert_vel, temp

real :: kappa
integer :: k

kappa = rdgas/cp_air

dmean_tot = 0.0

if(vert_difference_option == 'simmons_and_burridge') then

  do k=1,num_levels
    dp = dpk(k) + dbk(k)*p_surf_ref
    dp_inv = 1/dp
    dlog_1 = ln_p_half_ref(k+1) - ln_p_full_ref(k)
    dlog_3 = ln_p_half_ref(k+1) - ln_p_half_ref(k)
    dmean = div(:,:,k)*dp
    dt_t(:,:,k) =  - kappa*t_ref(k)*(dmean_tot*dlog_3 + dmean*dlog_1)*dp_inv
    dmean_tot = dmean_tot + dmean
    vert_vel(:,:,k+1) = - dmean_tot
  enddo

else if (vert_difference_option == 'mcm') then

  do k=1,num_levels
    dp = dpk(k) + dbk(k)*p_surf_ref
    dmean = div(:,:,k)*dp
    p_full_ref = 0.5*(pk(k+1) + pk(k)) + 0.5*(bk(k+1) + bk(k))*p_surf_ref
    dt_t(:,:,k) = - (kappa*t_ref(k)/p_full_ref)*(dmean_tot + 0.5*dmean)
    dmean_tot = dmean_tot + dmean
    vert_vel(:,:,k+1) = - dmean_tot
  enddo

endif

dt_p_surf = - dmean_tot

do k=2,num_levels
  vert_vel(:,:,k) = vert_vel(:,:,k) + dmean_tot*bk(k)
enddo

do k=2,num_levels
  temp(:,:,k) = - vert_vel(:,:,k)*(t_ref(k) - t_ref(k-1))
enddo

temp(:,:,1)            = 0.0
temp(:,:,num_levels+1) = 0.0

do k=1,num_levels
  dp = dpk(k) + dbk(k)*p_surf_ref
  dp_inv = 1/dp
  dt_t(:,:,k) = dt_t(:,:,k) + .5*dp_inv*(temp(:,:,k+1)+temp(:,:,k))
enddo

return
end subroutine linear_tp_tendency_3d
!-----------------------------------------------------------------------

subroutine linear_tp_tendency_1d(div, t_ref, p_surf_ref, ln_p_half_ref, ln_p_full_ref, dt_p_surf, dt_t)

real, intent(in), dimension (:) :: div
real, intent(in), dimension (:) :: t_ref, ln_p_half_ref, ln_p_full_ref
real, intent(in)  :: p_surf_ref

real, intent(out) :: dt_p_surf
real, intent(out), dimension (:) :: dt_t

complex, dimension(0:0,0:0,size(div,1)) :: div_3d, dt_t_3d
complex, dimension(0:0,0:0)           :: dt_p_surf_3d

div_3d(0,0,:) = cmplx(div,0.)

call linear_tp_tendency(div_3d, t_ref, p_surf_ref, ln_p_half_ref, ln_p_full_ref, dt_p_surf_3d, dt_t_3d)

dt_t      = real(dt_t_3d(0,0,:))
dt_p_surf = real(dt_p_surf_3d(0,0))

return
end subroutine linear_tp_tendency_1d
!-----------------------------------------------------------------------
subroutine implicit_end

if(.not.module_is_initialized) return

deallocate(ref_temperature_implicit)
deallocate(ref_ln_p_half, ref_ln_p_full, del_ln_p_half, del_ln_p_full)
deallocate(div_mat, wave_matrix, h, eigen, wavenumber, pk,  bk, dpk, dbk)
module_is_initialized = .false.

return
end subroutine implicit_end
!-----------------------------------------------------------------------

end module implicit_mod

max_stars_repo_name: JoyMonteiro/FMSDynamicalCore


text:
C    @(#)jacbfc.f	20.7 10/13/99
      subroutine jacbfc
C
C     Zero subscripts are used in arrays korder, kolum
C
      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/amtrx.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/beta.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/dcsln.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/ikk.inc'
      include 'ipfinc/optim1.inc'
      include 'ipfinc/tran.inc'
      include 'ipfinc/intbus.inc'
C
      kta=kt+ntota
      if (kt.gt.ntot) go to 290
      if (msw.eq.0) go to 100
      if (ikk(1,kta).eq.2) go to 110
      lp = 0
      kolum(0) = 0
      korder(lp) = 0
      max = 0
      go to 330

  100 continue
      jt=jckikk(kta,3)
c
c     Disable area interchange to improve initial convergence
c
c     if (jt.gt.0) go to 212

  110 max = 0
      kl=0
      lp=0
      lsw=1
      ek=e(kt)
      aik=0.0
      do 180 l = km(kt), km(kt)-1+kmlen(kt)
      mt = ikmu(l)
      mta=mt+ntota
      lp=lp+1
      em=e(mt)
      bkm = bkmu(l)
      rh=ek*em*bkm
      aik=aik+rh
      go to (130,140,150), lsw
  120 lp=lp-1
      go to 180

  130 if (mt.le.nbslck) go to 120
      lsw=2

  140 if (mt.lt.kt) go to 150
      max = kta
      kl=lp
      lsw=3
      korder(lp)=lp+1
      kolum(lp)=kta
      lp=lp+1

  150 if (msw.eq.0) go to 170
      if (ikk(1,mta).eq.2) go to 160
      lp = lp - 1
      go to 180

  160 rh = bkm
  170 max = mta
      kolum(lp) = mta
      rowh(lp)=-rh
      korder(lp)=lp+1
  180 continue
      if (kl.ne.0) go to 190
      kl=lp+1
      lp=kl
      kolum(lp)=kta
      korder(lp) = lp + 1
      max = kta
  190 if (msw.ne.0) aik = -bkku(kt)
      rowh(kl) = aik
C
C     Test for ideal bus tie branch
C
      ix = jckikk (kta,10)
      if (ix.eq.0) go to 210
C
C     Add a small residual value to diagonal to replace the branch
C     removed.
C
      rowh(kl) = rowh(kl) + 0.001
C
C     Add a column for tie entity
C
  200 lp = lp + 1
      korder(lp) = lp + 1
      iax = txtie(8,ix)
      max = iabs (iax)
      kolum(lp) = max
      rh = 1.0
      if (iax.lt.0) rh = -1.0
      if (msw.ne.0) rh = rh/ek
      rowh(lp) = rh
C
C     Check for additional bus ties
C
      ix = ix + 1
      if (ix.gt.ntxtie) go to 210
      if (txtie(1,ix).eq.kt) go to 200
  210 continue
      korder(lp)=0
      mend=lp+1
      korder(0)=1
      go to 330
C
C     Area Interchange Constraints
C
  212 max=0
      kolum(0)=0
      korder(0)=0
      lp=0
      mend=1
      korder(1)=1
      j1=karea(3,jt)
      j2=j1+karea(4,jt)-1
      do 280 j=j1,j2
      ix=kaloc(j)
      iax=iabs(ix)
      k1=tie(1,iax)
      k2=tie(7,iax)
      ka1=tie(2,iax)
      ka2=tie(8,iax)
      kdc=tie(9,iax)
      if (kdc.ne.0) go to 280
C
C     Process phase shifter tie lines. (The following logic is
C     vulnerable if K1 and K2 both have two or more bus ties!)
C
      i1 = jckikk(k1+ntota,10)
      i2 = jckikk(k2+ntota,10)
      if (i1 .eq. 0 .or. i2 .eq. 0) then
         go to 228
      else if (txtie(1,i1) .eq. txtie(2,i2) .and.
     1         txtie(2,i1) .eq. txtie(1,i2)) then
      else
         do 214 i1 = 1, ntxtie
         if (txtie(1,i1) .eq. k1 .and. txtie(2,i1) .eq. k2)
     1      go to 216
  214    continue
         go to 228
  216    continue
      endif

      if (ix .lt. 0) go to 280

      kx = dabs (txtie(8,i1))
      rh = 1.0
      if (txtie(8,i1) .lt. 0) rh = -rh
      if (ka1.ne.jt) rh = -rh

      lp = 0
      go to 222

  221 if (kolum(lp) - kx) 222,225,224

  222 ko=lp
      lp=korder(lp)
      if (lp.ne.0) go to 221
      max=kx

  224 korder(mend)=korder(ko)
      korder(ko)=mend
      kolum(mend)=kx
      ko=mend
      lp=mend
      mend=mend+1
      rowh(lp)=0.0

  225 continue
      rowh(lp) = rowh(lp) + rh
      ko = 0
      lp = korder(ko)
      go to 280

  228 mt=k1
      if (ix.lt.0) mt=k2
      if (mt.le.nbslck) go to 280
      mta=mt+ntota
  230 if (kolum(lp)-mta) 240,260,250
  240 ko=lp
      lp=korder(lp)
      if (lp.ne.0) go to 230
      max=mta
  250 korder(mend)=korder(ko)
      korder(ko)=mend
      kolum(mend)=mta
      ko=mend
      lp=mend
      mend=mend+1
      rowh(lp)=0
      go to 270

  260 ko=lp
  270 rh=e(k1)*e(k2)*tie(6,iax)
      if (ix.lt.0) rh=-rh
      if (ka1.ne.jt) rh = -rh
      rowh(lp)=rowh(lp)+rh
  280 continue
      go to 330
C
C     Add bus tie constraints
C
  290 ix = ikk(5,kta)
      if (ix.le.0) call erexit
      k1 = dmin1 (txtie(1,ix),txtie(2,ix))
      k2 = dmax1 (txtie(1,ix),txtie(2,ix))
      lp = 0
      if (k1.le.nbslck) go to 310
      if (msw.eq.0) go to 300
      if (ikk(1,k1+ntota).eq.1) go to 310
  300 lp=lp+1
      ko=lp
      kolum(lp) = k1 + ntota
      max = kolum(lp)
      korder(lp) = lp + 1
      if (k1 .eq. txtie(1,ix)) then
         rowh(lp) = 1.0
      else
         rowh(lp) = -1.0
         if (msw .ne .0) rowh(lp) = -txtie(3,ix)
      endif
      if (k1 .eq. k2) go to 320
  310 k1 = k2
      if (msw.eq.0) go to 300
      if (ikk(1,k1+ntota).eq.2) go to 300

  320 if (msw .eq. 0) then

         lp = lp + 1
         ko = lp
         kolum(lp) = kta
         max = kolum(lp)
         korder(lp) = lp + 1

         k1 = txtie(1,ix)
         k2 = txtie(2,ix)

         cx = e(k1) * e(k2) * txtie(10,ix)
         rowh(lp) = 1.0 / cx

      else

         if (lp .eq. 0) then
C
C           A bus tie spans two PV nodes. Change into pseudo-constraint.
C
            lp = 1
            ko = lp
            kolum(lp) = kta
            max = kolum(lp)
            rowh(lp) = 1.0
         endif

      endif
      mend = lp + 1
      korder(0) = 1
      korder(lp) = 0
  330 continue
      if (kolum(lp) .ne. max) then
         lp = 0
  340    ko = lp
         lp = korder(lp)
         if (lp .eq. 0) call erexit
         if (kolum(lp) .lt. max) go to 340
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
module sm_module

  use module_with_ptr_object
  use const_scalar_func_type
  use scalar_func_factories_EXPORTS_NOTHING
  use my_func_type

contains

  subroutine solid_mech_init
    call define_ptr_object_and_contents
    call get_ptr_object_func
  end subroutine

  subroutine define_ptr_object_and_contents

    integer :: stat
    class(scalar_func), allocatable :: unused1, unused2, dummy, rho
    character(80) :: errmsg
    character(:), allocatable :: unused3

    allocate(const_scalar_func :: dummy)

    call alloc_my_func (rho, stat, errmsg)
    if (stat /= 0) errmsg = 'x' // trim(errmsg) ! NOT EXECUTED
    call move_func_to_ptr_component (rho)

  end subroutine

  subroutine get_ptr_object_func
    class(scalar_func), allocatable :: f
    call get_copy_of_ptr_func_component (f)
  end subroutine

end module sm_module

program main
  use sm_module
  call solid_mech_init
end program

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
C    @(#)pflend.f	20.3 2/13/96
      subroutine pflend
C
C       Powerflow_end job processing
C
C       Added case ID to PRGMON call - AHS 4/22/82

      include 'ipfinc/blank.inc'
      include 'ipfinc/jobctl.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/usranl.inc'

      character job*30

      job = 'FINIS PWRFLW '//prgvsn

C     Add the case ID to the job label

      job(21:30) = chase1(1)
      call prtime('POWERFLOW')
C
C
      if ( kspare(16) .eq. -1 ) then
           close ( unit = mfich, status = 'DELETE' )
      else
           close ( unit = mfich, status = 'KEEP' )
      endif

      if ( kspare(17) .eq. 0 .and. batch) then
           close ( unit = dbug, status = 'DELETE' )
      else
           close ( unit = dbug, status = 'KEEP' )
      endif

      if (numusr .gt. 0 .and. lunusr .gt. 0) then
          close (unit = lunusr, status = 'DELETE')
      endif
C
      return
      end


text:
      SUBROUTINE ADVECT(DTA,FORCED)
C
C THIS IS THE MAIN COMPUTATIONAL PROCEDURE FOR THE ADVECTION
C EQUATION THAT ADVANCES ONE TIME-
C STEP. THIS TIMESTEPPING IS DONE FOR THE PROGNOSTIC VARIABLE
C GEOPOTENTIAL USING A LEAPFROG TIMESTEPPING ALGORITHM.
C DIVERGENCE AND VORTICITY ARE JUST COPIED FOR SIMPLICITY.
C THE OLD TIMELEVEL IS LNM1, THE DERIVATIVE IS EVALUATED AT
C TIMELEVEL LN, AND THE COMPUTED VALUES GO INTO TIMELEVEL LNP1.
C THE NONLINEAR PRODUCTS ARE EVALUATED AT THE GRIDPOINTS, AND
C A FFT TRANSFORM (FFT991) IS USED FOR EACH LATITUDE.
C SHTRNS INVERSE TRANSFORMS TO GRIDPOINT SPACE.
C
C CALLED BY: COMP1
C CALLS: FFT991, SHTRNS
C
C REVISIONS:
C 7-10-92 CHANGE TO CCM CODING CONVENTIONS (R. JAKOB)
C
C---- Model Parameters -------------------------------------------------
C
      INCLUDE 'params.i'
C
C---- Common Blocks ----------------------------------------------------
C
C TRANSFORM ARRAYS
      INCLUDE 'trnsfm.i'
C TIME DEPENDENT FIELDS
      INCLUDE 'tdvars.i'
C WORKSPACE
      INCLUDE 'wrkspc.i'
C
C------------ Arguments ------------------------------------------------
C
C     Input
C
C TIMESTEP IN SECONDS
      REAL DTA
C FLAG FOR FORCING TERMS
      LOGICAL FORCED
C
C------ Local Variables ------------------------------------------------
C
C     EQUIVALENCE TO STORAGE POOL IN COMMON BLOCK WRKSPC WHERE POSSIBLE
C
      REAL TRMC(NLON+2,NLAT), TRMD(NLON+2,NLAT),
     $     TRMG(NLON+2,NLAT)
      COMPLEX TRMCFC(NFC,NLAT), TRMDFC(NFC,NLAT),
     $        TRMGFC(NFC,NLAT)
      EQUIVALENCE (ETAFCG(1,1), WS1(1,1)),
     $            (DIVFCG(1,1), WS2(1,1)),
     $            (TRMC(1,1), TRMCFC(1,1), WS3(1,1)),
     $            (TRMD(1,1), TRMDFC(1,1), WS4(1,1)),
     $            (PHIFCG(1,1), WS5(1,1)),
     $            (TRMG(1,1), TRMGFC(1,1), WS6(1,1))
C
C     FORCING TERMS (PHYSICAL & SPECTRAL)
C
      REAL ETAFCG(NLON+2,NLAT), DIVFCG(NLON+2,NLAT),
     $     PHIFCG(NLON+2,NLAT)
C
C     TEMPORARIES
C
      INTEGER I,J,L
C
C----- Executable Statements -------------------------------------------
C
C     ADVECTION ROUTINE
C
C     ZERO TERMS FOR LATER ACCUMULATIONS
C
      DO 815 L=1,NALP
         PHISC(L) = (0.0,0.0)
  815 CONTINUE
C
C     COMPUTE THE GEOPOTENTIAL FORCING TERM
C
      IF (FORCED) THEN
         CALL FORCE(ETAFCG,DIVFCG,PHIFCG)
      ENDIF
C
C     EVALUATE NON-LINEAR ADVECTION TERMS
C     AND OLD TIMELEVEL OF DIVERGENCE & GEOPOTENTIAL
C
      DO 20 J=1,NLAT
         DO 10 I=1,NLON
            TRMC(I,J) = UCOS(I,J)*PHI(I,J,LN)
            TRMD(I,J) = VCOS(I,J)*PHI(I,J,LN)
            TRMG(I,J) = PHI(I,J,LNM1) - DTA*PHIBAR*DIV(I,J,LN)
            IF (FORCED) THEN
               TRMG(I,J) = TRMG(I,J) + DTA*PHIFCG(I,J)
            ENDIF
   10    CONTINUE
   20 CONTINUE
C
C     FOURIER TRANSFORM NON-LINEAR TERMS IN PLACE
C
      CALL FFT991(TRMC,FFTWS2,TRIGS,IFAX,1,NLON+2,NLON,NLAT,-1)
      CALL FFT991(TRMD,FFTWS2,TRIGS,IFAX,1,NLON+2,NLON,NLAT,-1)
      CALL FFT991(TRMG,FFTWS2,TRIGS,IFAX,1,NLON+2,NLON,NLAT,-1)
C
C     FWD LEGENDRE TRANSFORM NONLINEAR TERMS IN CONTINUITY EQUATION
C     (BECOMES ADVECTION EQUATION FOR NONDIVERGENT FIELD)
C
      CALL FTRNPE (DTA, TRMCFC, TRMDFC, TRMGFC, PHISC)
C
C     TRANSFORM ADVECTED GEOPOTENTIAL BACK TO GRID SPACE
C
      CALL SHTRNS(LVLS,LNP1,+1,PHI,PHISC)
C
C     COPY VORTICITY & DIVERGENCE FIELD
C     (SO THERE AREN'T MORE CHANGES TO PROGRAM STRUCTURE)
C
      DO 710 J=1,NLAT
         DO 700 I=1,NLON
            DIV(I,J,LNP1) = DIV(I,J,LNM1)
            ZETA(I,J,LNP1) = ZETA(I,J,LNM1)
  700    CONTINUE
  710 CONTINUE
C
C     FINISHED FOR ADVECTION EQUATION
C
      RETURN
C
      END


text:
C    Copyright(C) 1988-2017 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
      SUBROUTINE HEADER (NDIM, TITLE, NUMEL, NUMNP, AXI, GENFIL)
C
      include 'nu_io.blk'
      CHARACTER*16 FORM(3)
      CHARACTER*80 TITLE
      CHARACTER*(*) GENFIL
      LOGICAL AXI, FIRST
      CHARACTER*6 DIMEN(3)

      DATA DIMEN/'One',   'Two',   'Three'/
      DATA FORM /', Planar',', Axisymmetric',' '/
      DATA FIRST /.TRUE./
C
      IO = IHARD
      IF (FIRST) THEN
          WRITE (IO, 10) GENFIL(:LENSTR(GENFIL))
   10     FORMAT (5X,'Genesis: ',A/)
          WRITE (IO, 20) TITLE(:LENSTR(TITLE))
   20     FORMAT (5X,'Title:   ',A/)
          ILAB = 2
          IF (NDIM .EQ. 2 .AND. .NOT. AXI) ILAB = 1
          IF (NDIM .EQ. 3) ILAB = 3
          WRITE (IO, 30) NUMNP,
     *        NUMEL, DIMEN(NDIM)(:LENSTR(DIMEN(NDIM))),
     *        FORM(ILAB)(:LENSTR(FORM(ILAB)))
   30     FORMAT (5X,'Number of Nodes:    ',I10,/
     *        5X,'Number of Elements: ',I10/
     *        5X,A,'-Dimensional Mesh',A)
          FIRST = .FALSE.
      END IF
      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
! ***********************************************************************
!
!   Copyright (C) 2010  Bill Paxton
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa 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 library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************

      module run_star_extras

      use star_lib
      use star_def
      use const_def
      use crlibm_lib

      ! use star_utils, only: tau_eff

      implicit none

      ! these routines are called by the standard run_star check_model
      contains

      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         ! this is the place to set any procedure pointers you want to change
         ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)

         ! Uncomment these lines if you wish to use the functions in this file,
         ! otherwise we use a null_ version which does nothing.
         s% extras_startup => extras_startup
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns

         ! Once you have set the function pointers you want,
         ! then uncomment this (or set it in your star_job inlist)
         ! to disable the printed warning message,
          s% job% warn_run_star_extras =.false.

      end subroutine extras_controls

      ! None of the following functions are called unless you set their
      ! function point in extras_control.

      integer function extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_startup = 0
         if (.not. restart) then
            call alloc_extra_info(s)
         else ! it is a restart
            call unpack_extra_info(s)
         end if
      end function extras_startup

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going
         if (.false. .and. s% star_mass_h1 < 0.35d0) then
            ! stop when star hydrogen mass drops to specified level
            extras_check_model = terminate
            write(*, *) 'have reached desired hydrogen mass'
            return
         end if

         ! if you want to check multiple conditions, it can be useful
         ! to set a different termination code depending on which
         ! condition was triggered.  MESA provides 9 customizeable
         ! termination codes, named t_xtra1 .. t_xtra9.  You can
         ! customize the messages that will be printed upon exit by
         ! setting the corresponding termination_code_str value.
         ! termination_code_str(t_xtra1) = 'my termination condition'

         ! Switch atmosphere option below a certain temperature.
         if(s% Teff < s% x_ctrl(10)) then
            s% which_atm_option = 'WD_tau_25_tables'
            s% tau_factor = 1
         end if

         ! by default, indicate where (in the code) MESA terminated
         if (extras_check_model == terminate) s% termination_code = t_extras_check_model
      end function extras_check_model

      integer function how_many_extra_history_columns(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 18
      end function how_many_extra_history_columns

      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
         use ionization_lib, only: eval_typical_charge
         use chem_lib, only: chem_get_iso_id
         use chem_def
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         integer :: k, k_cz_bot, itrace
         real(dp) :: cz_xm, edv_iso
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         ! Id in the NET for the iso we wish to trace.
         itrace = s% net_iso(ica40) ! e.g. s% net_iso(ica40) for Calcium

         !note: do NOT add the extras names to history_columns.list
         ! the history_columns.list is only for the built-in log column options.
         ! it must not include the new column names you are adding here.

         ! Mass exterior to tau_eff = 5
         ! Turn this one into a dummy by default, since compiling requires star_utils module.
         names(1) = 'photosphere5_xm_eff'
         vals(1) = -1
         !do k = 2,s% nz ! first k is not a tau_eff
         !   if(tau_eff(s,k) > 5d0) then
         !      vals(1) = s% xmstar*sum(s% dq(1:k-1))/Msun
         !      exit
         !   end if
         !end do

         ! Mass exterior to tau = 5
         names(2) = 'photosphere5_xm'
         do k = 2,s% nz ! first k is not a tau_eff
            if(s% tau(k) > 5d0) then
               vals(2) = s% xmstar*sum(s% dq(1:k-1))/Msun
               exit
            end if
         end do

         ! Mass exterior to tau = 2/3
         names(3) = 'photosphere23_xm'
         do k = 2,s% nz ! first k is not a tau_eff
            if(s% tau(k) > 2d0/3d0) then
               vals(3) = s% xmstar*sum(s% dq(1:k-1))/Msun
               exit
            end if
         end do

         ! Mass exterior to tau = 1
         names(4) = 'photosphere1_xm'
         do k = 2,s% nz ! first k is not a tau_eff
            if(s% tau(k) > 1d0) then
               vals(4) = s% xmstar*sum(s% dq(1:k-1))/Msun
               exit
            end if
         end do

         names(5) = 'evan_cz_xm'
         ! Find cell where switches from convective to not
         k_cz_bot = 1
         do k=2,s% nz
            if(s% mixing_type(k-1) == convective_mixing .and. &
                 s% mixing_type(k) /= convective_mixing) then
               k_cz_bot = k - 1
               exit
            end if
         end do
         ! print *, k_cz_bot
         ! print *, s% cz_bdy_dq(k_cz_bot)
         ! vals(5) = s% xmstar*(sum(s% dq(1:k_cz_bot-1)) - s% cz_bdy_dq(k_cz_bot))/Msun

         ! k_cz_bot is the last convective cell,
         ! the face of k_cz_bot + 1 is the face at the boundary of the convective zone
         ! (roughly, can also try to interpolate within cells to find exact boundary)
         cz_xm = s% xmstar*sum(s% dq(1:k_cz_bot))/Msun
         vals(5) = cz_xm

         ! Also use location of base of the convection zone to evaluate tyipcal charges there.
         names(6) = 'base_charge_Na'
         names(7) = 'base_charge_Mg'
         names(8) = 'base_charge_Si'
         names(9) = 'base_charge_Ca'
         names(10) = 'base_charge_Fe'
         vals(6) = s% typical_charge(s% net_iso(ina23),k_cz_bot)
         vals(7) = s% typical_charge(s% net_iso(img24),k_cz_bot)
         vals(8) = s% typical_charge(s% net_iso(isi28),k_cz_bot)
         vals(9) = s% typical_charge(s% net_iso(ica40),k_cz_bot)
         vals(10) = s% typical_charge(s% net_iso(ife56),k_cz_bot)

         ! Most quantities here are evaluated at the face of k_cz_bot+1,
         ! the outer edge of the first non-convective cell.
         ! See MESA1 fig 9.
         edv_iso = abs(s% edv(itrace,k_cz_bot+1))
         names(11) = 'inst_diffusion_timescale'
         vals(11) = cz_xm*Msun / ( pi4*pow2(s% r(k_cz_bot+1))*s% rho_face(k_cz_bot+1)*edv_iso)

         names(12) = 'steady_diffusion_timescale'
         ! Backs up a bit into convective zone to evaluate mass fraction.
         vals(12) = -1
         ! vals(12) = &
         !      s% xa(itrace,k_cz_bot-2) * &
         !      cz_xm*Msun / (s% star_mdot * Msun/secyer)

         names(13) = 'cz_bot_zone'
         vals(13) = k_cz_bot

         names(14) = 'edv_base'
         vals(14) = edv_iso

         ! I need MESA timescales for the 4 dominant elements to see how my ratios compare.
         itrace = s% net_iso(img24)
         edv_iso = abs(s% edv(itrace,k_cz_bot+1))
         names(15) = 'inst_Mg'
         vals(15) = cz_xm*Msun / ( pi4*pow2(s% r(k_cz_bot+1))*s% rho_face(k_cz_bot+1)*edv_iso)

         itrace = s% net_iso(isi28)
         edv_iso = abs(s% edv(itrace,k_cz_bot+1))
         names(16) = 'inst_Si'
         vals(16) = cz_xm*Msun / ( pi4*pow2(s% r(k_cz_bot+1))*s% rho_face(k_cz_bot+1)*edv_iso)

         itrace = s% net_iso(io16)
         edv_iso = abs(s% edv(itrace,k_cz_bot+1))
         names(17) = 'inst_O'
         vals(17) = cz_xm*Msun / ( pi4*pow2(s% r(k_cz_bot+1))*s% rho_face(k_cz_bot+1)*edv_iso)

         itrace = s% net_iso(ife56)
         edv_iso = abs(s% edv(itrace,k_cz_bot+1))
         names(18) = 'inst_Fe'
         vals(18) = cz_xm*Msun / ( pi4*pow2(s% r(k_cz_bot+1))*s% rho_face(k_cz_bot+1)*edv_iso)

      end subroutine data_for_extra_history_columns

      integer function how_many_extra_profile_columns(id, id_extra)
         use star_def, only: star_info
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 3
      end function how_many_extra_profile_columns

      subroutine data_for_extra_profile_columns(id, id_extra, n, nz, names, vals, ierr)
         use star_def, only: star_info, maxlen_profile_column_name
         use const_def, only: dp
         use chem_def
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: j, k, nc, m, hi, oi, fei
         ! These must have dimension matching the net (pollution3.net) + 1 for electrons.
         real(dp), dimension(18) :: A, charge, nd
         real(dp), dimension(18,18) :: Kdiff,zdiff,zdiff1,zdiff2
         real(dp) :: kappa, DHFe, DHO
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         !note: do NOT add the extra names to profile_columns.list
         ! the profile_columns.list is only for the built-in profile column options.
         ! it must not include the new column names you are adding here.

         ! here is an example for adding a profile column
         !if (n /= 1) stop 'data_for_extra_profile_columns'
         !names(1) = 'beta'
         !do k = 1, nz
         !   vals(k,1) = s% Pgas(k)/s% P(k)
         !end do

         ! print *, "saving extra profile_columns"
         ! Number of species in pollution3.net
         nc = 17
         m = nc+1
         if(nc /= s% species) then
            stop "net doesn't match assumed size in run_star_extras"
         end if

         names(1) = 'kappa_SM'
         names(2) = 'D12_H_Fe'
         names(3) = 'D12_H_O'

         A(1:nc) = chem_isos% Z_plus_N(s% chem_id(1:nc))
         A(m) = me/amu

         ! indices for various species
         hi = s% net_iso(ih1)
         oi = s% net_iso(io16)
         fei = s% net_iso(ife56)
         ! print *, "Iron index:", fei

         do k=1,s%nz

            do j=1,nc
               charge(j) = max( 1d0, s% typical_charge(j,k) )
               nd(j) = s% rho(k) * s% xa(j,k) / (A(j)*amu)
            end do
            ! print *, "Charge for 5th net iso: ", charge(5)
            ! Electron entries
            charge(m) = -1d0
            nd(m) = 0d0
            do j=1,nc
               nd(m) = nd(m) + charge(j)*nd(j)
            end do

            call get_SM_coeffs( &
                 nc,m, &
                 s% rho(k), s% T(k), &
                 A,charge,nd, &
                 Kdiff,zdiff,zdiff1,zdiff2,kappa)

            ! D12 = ( n1 n2/(n1 + n2) ) k T / K12
            DHFe = nd(hi)*nd(fei)*boltzm*s% T(k) / ( (nd(hi) + nd(fei)) * Kdiff(hi,fei))
            DHO = nd(hi)*nd(oi)*boltzm*s% T(k) / ( (nd(hi) + nd(oi)) * Kdiff(hi,oi))

            vals(k,1) = kappa
            vals(k,2) = DHFe
            vals(k,3) = DHO
         end do

      end subroutine data_for_extra_profile_columns

      ! returns either keep_going or terminate.
      ! note: cannot request retry or backup; extras_check_model can do that.
      integer function extras_finish_step(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_finish_step = keep_going
         call store_extra_info(s)

         ! to save a profile,
            ! s% need_to_save_profiles_now = .true.
         ! to update the star log,
            ! s% need_to_update_history_now = .true.

         ! see extras_check_model for information about custom termination codes
         ! by default, indicate where (in the code) MESA terminated
         if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step
      end function extras_finish_step

      subroutine extras_after_evolve(id, id_extra, ierr)
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine extras_after_evolve

      ! routines for saving and restoring extra data so can do restarts

         ! put these defs at the top and delete from the following routines
         !integer, parameter :: extra_info_alloc = 1
         !integer, parameter :: extra_info_get = 2
         !integer, parameter :: extra_info_put = 3

      subroutine alloc_extra_info(s)
         integer, parameter :: extra_info_alloc = 1
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_alloc)
      end subroutine alloc_extra_info

      subroutine unpack_extra_info(s)
         integer, parameter :: extra_info_get = 2
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_get)
      end subroutine unpack_extra_info

      subroutine store_extra_info(s)
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_put)
      end subroutine store_extra_info

      subroutine move_extra_info(s,op)
         integer, parameter :: extra_info_alloc = 1
         integer, parameter :: extra_info_get = 2
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         integer, intent(in) :: op

         integer :: i, j, num_ints, num_dbls, ierr

         i = 0
         ! call move_int or move_flg
         num_ints = i

         i = 0
         ! call move_dbl

         num_dbls = i

         if (op /= extra_info_alloc) return
         if (num_ints == 0 .and. num_dbls == 0) return

         ierr = 0
         call star_alloc_extras(s% id, num_ints, num_dbls, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in star_alloc_extras'
            write(*,*) 'alloc_extras num_ints', num_ints
            write(*,*) 'alloc_extras num_dbls', num_dbls
            stop 1
         end if

         contains

         subroutine move_dbl(dbl)
            real(dp) :: dbl
            i = i+1
            select case (op)
            case (extra_info_get)
               dbl = s% extra_work(i)
            case (extra_info_put)
               s% extra_work(i) = dbl
            end select
         end subroutine move_dbl

         subroutine move_int(int)
            integer :: int
            i = i+1
            select case (op)
            case (extra_info_get)
               int = s% extra_iwork(i)
            case (extra_info_put)
               s% extra_iwork(i) = int
            end select
         end subroutine move_int

         subroutine move_flg(flg)
            logical :: flg
            i = i+1
            select case (op)
            case (extra_info_get)
               flg = (s% extra_iwork(i) /= 0)
            case (extra_info_put)
               if (flg) then
                  s% extra_iwork(i) = 1
               else
                  s% extra_iwork(i) = 0
               end if
            end select
         end subroutine move_flg

       end subroutine move_extra_info

       ! Calculate coefficients given in Appendix C.3 of Stanton & Murillo, PR E 93, 043203 (2016)
      subroutine get_SM_coeffs(nc,m,rho,T,A,Z,nd,Kdiff,zdiff,zdiff1,zdiff2,kappa)
        integer, intent(in) :: nc, m
        real(dp), intent(in) :: rho, T
        real(dp), dimension(:), intent(in) :: A, Z, nd ! m
        real(dp), intent(out) :: kappa ! ion separation relative to electron screening length.
        real(dp), dimension(:,:), intent(inout) :: Kdiff,zdiff,zdiff1,zdiff2 ! (m,m) but electron entries shouldn't be used.
        ! Only the ion terms are modified by this routine right now.

        real(dp), dimension(nc,nc) :: g_plasma, mu
        real(dp), dimension(nc,nc,2,3) :: Keff, Omega ! Indexed by all species and different orders
        ! Fitting coefficients from Stanton & Murillo
        real(dp), dimension(2,3) :: a1,a2,a3,a4,a5,b0,b1,b2,b3,b4
        real(dp) :: lambda ! The screening length
        integer :: i,j,k,facmo,no,mo ! Last two are used for different orders of collisions.

        ! Initialize all to 0, since some entries never get set or used.
        a1(:,:) = 0d0
        a2(:,:) = 0d0
        a3(:,:) = 0d0
        a4(:,:) = 0d0
        a5(:,:) = 0d0
        b0(:,:) = 0d0
        b1(:,:) = 0d0
        b2(:,:) = 0d0
        b3(:,:) = 0d0
        b4(:,:) = 0d0

        ! Table IV, coefficients for fits (C23)-(C24)
        a1(1,1) = 1.4660d0
        a1(1,2) = 0.52094d0
        a1(1,3) = 0.30346d0
        a1(2,2) = 0.85401d0

        a2(1,1) = -1.7836d0
        a2(1,2) = 0.25153d0
        a2(1,3) = 0.23739d0
        a2(2,2) = -0.22898d0

        a3(1,1) = 1.4313d0
        a3(1,2) = -1.1337d0
        a3(1,3) = -0.62167d0
        a3(2,2) = -0.60059d0

        a4(1,1) = -0.55833d0
        a4(1,2) = 1.2155d0
        a4(1,3) = 0.56110d0
        a4(2,2) = 0.80591d0

        a5(1,1) = 0.061162d0
        a5(1,2) = -0.43784d0
        a5(1,3) = -0.18046d0
        a5(2,2) = -0.30555d0

        b0(1,1) = 0.081033d0
        b0(1,2) = 0.20572d0
        b0(1,3) = 0.68375d0
        b0(2,2) = 0.43475d0

        b1(1,1) = -0.091336d0
        b1(1,2) = -0.16536d0
        b1(1,3) = -0.38459d0
        b1(2,2) = -0.21147d0

        b2(1,1) = 0.051760d0
        b2(1,2) = 0.061572d0
        b2(1,3) = 0.10711d0
        b2(2,2) = 0.11116d0

        b3(1,1) = -0.50026d0
        b3(1,2) = -0.12770d0
        b3(1,3) = 0.10649d0
        b3(2,2) = 0.19665d0

        b4(1,1) = 0.17044d0
        b4(1,2) = 0.066993d0
        b4(1,3) = 0.028760d0
        b4(2,2) = 0.15195d0

        ! Get the screening length
        call lam_SM(nc,m,rho,T,Z,nd,lambda,kappa)

        ! Calculate g_plasma
        do i=1,nc
           do j=1,nc
              g_plasma(i,j) = Z(i)*Z(j)*qe*qe/(lambda*boltzm*T)
           end do
        end do

        ! Calculate Keff using Eqns (C22-C24)
        do i=1,nc
           do j=1,nc
              if( g_plasma(i,j) < 0d0) then ! Don't calculate for attractive potentials, set to 0
                 Keff(i,j,:,:) = 0d0
              else if( g_plasma(i,j) < 1d0) then ! Use eqn C23 for weakly coupled
                 do no=1,2
                    do mo=1,3 ! Implementing the (m-1)! term with a simple if statement.
                       if(mo .eq. 3) then
                          facmo = 2 ! (3-1)!
                       else
                          facmo = 1 ! (1-1)! and (2-1)!
                       end if
                       Keff(i,j,no,mo) = (-1d0*no/4d0)*facmo*safe_log_cr( &
                            a1(no,mo)*g_plasma(i,j) &
                            + a2(no,mo)*pow2(g_plasma(i,j)) &
                            + a3(no,mo)*pow3(g_plasma(i,j)) &
                            + a4(no,mo)*pow4(g_plasma(i,j)) &
                            + a5(no,mo)*pow5(g_plasma(i,j)) )
                    end do
                 end do
              else ! Use eqn C24 for strongly coupled
                 do no=1,2
                    do mo=1,3
                       Keff(i,j,no,mo) = &
                            (b0(no,mo) + b1(no,mo)*safe_log_cr(g_plasma(i,j)) &
                            + b2(no,mo)*pow2(safe_log_cr(g_plasma(i,j))) ) / &
                            (1d0 + b3(no,mo)*g_plasma(i,j) + b4(no,mo)*pow2(g_plasma(i,j)) )
                    end do
                 end do
              end if
           end do
        end do

        ! Calculate the collision integrals using (C19)
        do i=1,nc
           do j=1,nc
              mu(i,j) = amu*A(i)*A(j)/(A(i) + A(j)) ! Reduced mass for collision
              do no=1,2
                 do mo=1,3
                    Omega(i,j,no,mo) = &
                         sqrt(2d0*pi/(mu(i,j))) * &
                         (pow2(Z(i)*Z(j)*qe*qe)/pow_cr(boltzm*T,1.5d0)) * &
                         Keff(i,j,no,mo)
                 end do
              end do
           end do
        end do

        ! Resistance coefficient is (16/3) ni nj mu Omega11
        ! (compare Paquette eqn 22 with SM eq B8)
        ! The zdiffs are given in terms of collision integrals in Paquette eqns (14-16), (23-25)
        do i=1,nc
           do j=1,nc
              Kdiff(i,j) = (16d0/3d0)*nd(i)*nd(j)*mu(i,j)*Omega(i,j,1,1)
              zdiff(i,j) = 1d0 - 2d0*Omega(i,j,1,2)/(5d0*Omega(i,j,1,1))
              zdiff1(i,j) = 2.5d0 + (2d0*Omega(i,j,1,3) - 10d0*Omega(i,j,1,2))/(5d0*Omega(i,j,1,1))
              zdiff2(i,j) = Omega(i,j,2,2)/Omega(i,j,1,1)
           end do
        end do

        ! Note that SM don't give fits for attractive potentials, so we haven't
        ! touched the electron entries. They exit this routine unchanged, so they
        ! either need to be initialized before this routine is called or somehow
        ! calculated later.

      end subroutine get_SM_coeffs

      ! Screening Length according to Stanton & Murillo
      subroutine lam_SM(nc,m,rho,T,Z,nd,lam_eff,kappa)
        integer, intent(in) :: nc, m
        real(dp), intent(in) :: rho, T
        real(dp), dimension(:), intent(in) :: Z, nd ! m. charges, number densities of all species
        real(dp), intent(out) :: lam_eff
        real(dp), intent(out) :: kappa ! ion separation relative to electron screening length.

        real(dp) :: tiny_n, ne, EF, lam_e, lam_sum, rhotot, ni_tot
        real(dp), dimension(nc) :: ai, gam_is, lam_i ! ion sphere radius, coupling, screening for each type of ion
        integer :: i

        tiny_n = 1d-20 ! g/cc
        ne = nd(m)
        ! Electron Fermi energy
        EF = (hbar*hbar*pow_cr(3d0*pi*pi*ne,2d0/3d0))/(2d0*me)
        ! Electron screening length accounting for degeneracy correction
        lam_e = pow_cr(pi4*qe*qe*ne/sqrt(pow2(boltzm*T) + pow2((2d0/3d0)*EF)),-0.5d0)

        ! Compute kappa
        ni_tot = 0d0
        do i = 1,nc
           ni_tot = ni_tot + nd(i)
        end do
        kappa = pow_cr(3d0/(pi4*ni_tot),1d0/3d0)/lam_e

        rhotot = 0d0
        do i = 1,nc
           ! rhotot is a CHARGE density. Just trying to follow the notation of SM eq (34)
           rhotot = rhotot + Z(i)*qe*nd(i) ! = qe*ne?
        end do
        do i = 1,nc
           ai(i) = pow_cr(3d0*Z(i)*qe/(pi4*rhotot),1d0/3d0)
           gam_is(i) = Z(i)*Z(i)*qe*qe/(ai(i)*boltzm*T)
           ! Number densities that are 0 or tiny cause screening length to diverge.
           ! This is physical; nothing is there to screen anything.
           ! But numerically, I don't want to rely on fortran to handle division by zero and infinity,
           ! so just set these screening lengths to be huge and they won't
           ! contribute anything to overall screening length.
           if(nd(i) < tiny_n) then
              lam_i(i) = 1d99 ! cm. This won't contribute to any screening.
           else
              lam_i(i) = sqrt(boltzm*T/( pi4*pow2(Z(i)*qe)*nd(i) ))
           end if
        end do

        lam_sum = 1d0/pow2(lam_e) ! The electron part of the screening length.
        do i = 1,nc ! Sum over all the ions.
           lam_sum = lam_sum + pow_cr( pow2(lam_i(i))*(1d0+3d0*gam_is(i)), -1d0)
        end do
        lam_eff = pow_cr(lam_sum,-0.5d0)
      end subroutine lam_SM

      end module run_star_extras


text:
!RUN: %flang -g -S -emit-llvm %s -o - | FileCheck %s

!CHECK: !DIDerivedType(tag: DW_TAG_member, name: "proc1"
!CHECK-SAME: baseType: [[SPTRTYPE:![0-9]+]]
!CHECK: [[SPTRTYPE]] = !DIDerivedType(tag: DW_TAG_pointer_type, baseType: [[SUBTYPE:![0-9]+]]
!CHECK: [[SUBTYPE]] = !DISubroutineType(types: [[ARGLIST:![0-9]+]]
!CHECK: [[ARGLIST]] = !{null, [[ARG1:![0-9]+]]
!CHECK: [[ARG1]] = !DICompositeType(tag: DW_TAG_array_type
!CHECK-SAME: dataLocation:
!CHECK-SAME: associated:

program main
  interface
    subroutine sub1(arg11)
      integer, pointer :: arg11(:)
    end subroutine
  end interface
  type type1
    procedure(sub1), pointer, nopass :: proc1
  end type
  type type2
    type(type1) :: mem1
  end type type2
  type(type2) :: arg1
end program main

max_stars_repo_name: vigbalu/flang

max_issues_repo_name: vigbalu/flang

max_forks_repo_name: vigbalu/flang


text:
! Program :
! Author  : wansooha@gmail.com
! Date    :

        program test_cfg_write
        use cfgio_mod
        implicit none
        type(cfg_t):: cfg
        character(len=128):: str
        logical:: ls
        integer n
        cfg=parse_cfg("test_intpl.cfg")
        call cfg%get("section 2",'n',n)
        print*,'n=',n
        call cfg%get("section1",'intpl',str)
        print*,'intpl='//trim(str)
        call cfg%get("section1",'fs',str)
        print*,'fs='//trim(str)
        str=cfg%gets("section 2",'fi')
        print*,'fi='//trim(str)
        call cfg%get("section 2",'fv',str)
        print*,'fv='//trim(str)
        call cfg%get("section 3",'ls1',ls)
        print*,'ls1=',ls
        call cfg%get("section 3",'ls2',ls)
        print*,'ls2=',ls

        call cfg%write()
        call cfg%finalize()

        end

max_stars_repo_name: baklanovp/cfgio

max_issues_repo_name: baklanovp/cfgio

max_forks_repo_name: baklanovp/cfgio


text:
!> @mainpage KORC documentation
!! @author Leopoldo Carbajal @note KORC and its documentation are in
!! constant evolution. This makes it possible that there are issues
!! that need to be solved.  If you find an issue please report it
!! immediatly through the "issues" section of the github repository.
!!
!! @section s1 Introduction The Kinetic Orbit Runaway electrons Code
!! (KORC) is a full-orbit particle tracer that evolves relativistic
!! electrons in both analytical and grid-based pre-computed electric
!! and magnetic fields. KORC includes the effects of: the acceleration
!! of the electrons due to the electric field, synchrotron radiation
!! energy losses, and collisions with the background plasma containing
!! high-Z impurities.
!!
!! For details about the equations of motion of the relativistic
!! electrons followed by KOqRC we refer the user to Carbajal et
!! al. Phys. Plasmas 24, 042512 (2017) and Carbajal
!! and del-Castillo-Negrete, Nuclear Fusion, submitted
!! (2018). Also, for details about the KORC's synchrotron
!! radiation synthetic diagnostic we refer the user to Carbajal
!! and del-Castillo-Negrete, Plasma Phys. Controll. Fusion 59,
!! 124001 (2017).
!!
!! KORC is a modular Fortran 95 code that uses a hybrid MPI + open MP
!! parallelization paradigm to exploit multi-core nodes systems, such
!! as Cori and Edison NERSC systems (www.nersc.gov).
!!
!!
!! @section s2 Installation @warning KORC has been installed and
!! tested in systems with OS 10.6 and higher, Ubuntu 14.04, and SuSe
!! 12. Though we have paid special attention to portability when
!! developing KORC, there is no guarantee that KORC can be compiled,
!! ran, and is accurate in other systems not listed above. As a rule
!! of thumb, we recommend to perform any convenient benchmark tests
!! when compiling and running KORC in a new system, this to make sure
!! that the external libraries and compilers do not modify the
!! simulation results.
!!
!! @subsection s1s2 Getting started To compile and run KORC you will
!! need the following:

The HDF5
!! library for I/O.

The PSPLINES library for
!! interpolating the pre-computed electric and magnetic fields as well
!! aqs the plasma profiles.

The GNU compilers suite, or

!!
The INTEL compilers suite.

The corresponding open MP
!! library

The MPI or open MPI library

 Git, for obtaining the latest version
!! of KORC.

!!
!! The HDF5 and PSPLINES libraries need to be installed BEFORE
!! compiling KORC, or if these are already present in the system, they
!! need to be loaded to your environment and their installation paths
!! need to be added to the Makefile accordingly. We refer the user to
!! the documentation of HDF5 and PSPLINES for the specifics about
!! their installation.
!!
!! If you are installing the HDF5 and PSPLINES libraries in your
!! system, we recommend to perform a local installation. Then enter
!! the absolute path of the folder containing the "lib" and
!! "bin" folders of each library to the PSPLINE_INSTALL and
!! HDF5_INSTALL variables of the Makefile.
!!
!!
!! @note PSPLINES is known for having precision issues sometimes when
!! compiled using the INTEL compilers suite in Linux systems.  Please,
!! double-check that the interpolations are giving the correct
!! numbers.
!!
!! @subsection s2s2 Cloning and compiling KORC To obtain the latest
!! version of KORC you will need to clone the Github repository to the
!! system where you want to run KORC.  This can be done as follows:
!!

 Make sure you have all the needed external libraries and
!! compilers in place.
Using the terminal of your system,
!! cd to the directory where you want to save KORC.

!!
Type the following in the terminal "git
!! https://github.com/ORNL-Fusion/KORC.git". This will copy all
!! the KORC files from the Github repository to your local system.
!!
cd to the KORC-FO folder.
In the terminal type:
!! "./compile.sh GNU" if you are using GNU compilers, or
!! "./compile.sh INTEL" if you are using INTEL compilers.
!!

!!
!! @section s3 Running KORC If you are running KORC in a system that
!! uses a resource manager such as TORQUE or SLURM, you will need to
!! request the number of nodes and cores per node that you will need,
!! this along with any other resources such as wall time, memory per
!! node, etc. Decide the number of MPI processes (nmpi) and the number
!! of open MP threads per MPI process (nomp) that you will use in your
!! simulation.  Then, on the terminal or in your batch file type the
!! following:
!!
!!  mpirun -np nmpi -x OMP_NUM_THREADS=nomp ./bin/KORC
!! "full_path_to_input_file.korc" "full_path_to_the_outputs_folder"
!!
!!
!! The above example assumes that you are using the open MPI wrapper
!! mpirun to execute the binary "KORC" generated after
!! compiling KORC. This wrapper can be replaced by any other wrapper
!! that you might want to use. Notice that when using a different MPI
!! wrapper the way to pass the number of MPI processes (nmpi) and the
!! number of open MP threads per MPI (nomp) to KORC might vary.
!!
!! @section s4 Generating the KORC's input file.
!! KORC uses...

max_stars_repo_name: ORNL-Fusion/KORC

max_issues_repo_name: ORNL-Fusion/KORC

max_forks_repo_name: ORNL-Fusion/KORC