text:
	SUBROUTINE GCURVE  ( itype, np, xarray, yarray, npout, xeval,
     +			     yeval, iret )
C************************************************************************
C* GCURVE								*
C* 									*
C* This subroutine fits a curve to a set of input points and then	*
C* evaluates the curve for values of x.  The input points must be	*
C* strictly monotonic in x:						*
C* 									*
C*     xarray (1) < xarray (2) < ... < xarray (np)			*
C*			or						*
C*     xarray (1) > xarray (2) > ... > xarray (np)			*
C* 									*
C* This subroutine simply fits the points to the curve.  It differs	*
C* from GCYEVL in that it does not map the points into the plotting	*
C* coordinate system before the fitting is done.			*
C* 									*
C* GCURVE  ( ITYPE, NP, XARRAY, YARRAY, NPOUT, XEVAL, YEVAL, IRET )	*
C* 									*
C* Input parameters:							*
C*	ITYPE		INTEGER		Type of curve			*
C*				  	  1 = piecewise linear		*
C*				  	  2 = cubic spline		*
C*	NP		INTEGER		Number of input points		*
C* 	XARRAY (NP)	REAL		X input coordinates 		*
C* 	YARRAY (NP)	REAL		Y input coordinates 		*
C*	NPOUT		INTEGER		Number of evaluations		*
C*	XEVAL (NPOUT)	REAL		X evaluation coordinates 	*
C*									*
C* Output parameters:							*
C* 	YEVAL (NPOUT)	REAL		Y evaluated coordinates 	*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 9/90	GEMPAK 5.0				*
C* T. Lee/GSC		 9/97	Fixed typo, NIMONO -> NOMONO		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'ERROR.PRM'
	INCLUDE		'CURVE.CMN'
	INCLUDE		'XYDEF.CMN'
C*
	REAL		xarray (*), yarray (*), xeval (*), yeval (*)
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
C
C*	Check for enough points.
C
	IF  ( ( np .lt. 2 ) .or. ( npout .lt. 1 ) ) iret = NOPNTS
	IF  ((np .gt. MAXPTS) .or. (npout .gt. MAXPTS)) iret = NIPNTS
	IF  ( iret .ne. NORMAL ) RETURN
C
C*	Check for type of curve to fit.
C
	CALL GCTYPE ( itype, jtype, jend, ier )
C
C*	Move data into internal arrays sorted in increasing order.
C*	Return with error if data is not strictly monotonic.
C
	delx = xarray (2) - xarray (1)
	IF  ( delx .lt. 0. ) THEN
	    ibegin = np
	    iend   = 1
	    inc    = -1
	  ELSE
	    ibegin = 1
	    iend   = np
	    inc    = 1
	END IF
C*
	x (1) = xarray ( ibegin )
	y (1) = yarray ( ibegin )
	nknt  = 1
	xx    = x (1)
	DO  i = ibegin+inc, iend, inc
	    tx = xarray (i)
	    ty = yarray (i)
	    IF  ( ( .not. ERMISS (tx) ) .and.
     +            ( .not. ERMISS (ty) ) ) THEN
		IF  ( tx .le. xx ) THEN
		    iret = NOMONO
		    RETURN
		  ELSE
		    nknt     = nknt + 1
		    x (nknt) = tx
		    y (nknt) = ty
		    xx       = tx
		END IF
	    END IF
	END DO
C
C*	Check to see if there are still at least two points.
C
	IF  ( nknt .lt. 2 ) THEN
	    iret = NOPNTS
	    RETURN
	END IF
C
C*	Generate coefficients for evaluating curve.
C
	IF  ( jtype .eq. 1 ) THEN
	    CALL GCFLIN ( x, y, nknt, coeffs, iret )
	  ELSE IF ( ( jtype .eq. 2 ) .and. ( nknt .ge. 3 ) ) THEN
	    CALL GCFSPL ( x, y, nknt, jend, coeffs, iret )
	  ELSE
	    DO  i = 1, npout
		yeval (i) = RMISSD
	    END DO
	    RETURN
	END IF
	IF  ( iret .ne. NORMAL ) THEN
	    iret = NORMAL
	    DO  i = 1, npout
		yeval (i) = RMISSD
	    END DO
	    RETURN
	END IF
C
C*	Transform evaluating points to normalized coordinates.
C
	DO  i = 1, npout
	    xout (i) = xeval (i)
	END DO
C
C*	Evaluate curve.
C
	IF  ( jtype .eq. 1 ) THEN
	    CALL GCELIN ( x, nknt, coeffs, xout, npout, yout, iret )
	  ELSE IF  ( jtype .eq. 2 ) THEN
	    CALL GCESPL ( x, nknt, coeffs, xout, npout, yout, iret )
	END IF
C
C*	Move points into output array.
C
	DO  i = 1, npout
	    yeval (i) = yout (i)
	END DO
C*
	RETURN
	END


text:
module Mod_nsc_EndElmope_im
   use typre

   use Mod_nsc_BaseElmope
   use Mod_nsc_ComputeVariables
   use Mod_nsc_ComputeConvectionCoefficients
   use Mod_nsc_ComputeConvectionJacobianGradientCoefficients
   use Mod_nsc_ComputeDiffusionCoefficients
   use Mod_nsc_ComputeNonLinearDiffusionCoefficients
   use Mod_nsc_ComputeReactionCoefficients
   use Mod_nsc_ComputeTaus
   use Mod_nsc_ComputeTestf
   use Mod_nsc_NonLinearDerivatives
   use Mod_nsc_ExternalForces

   use Mod_nsc_ComputeGpResidual
   use Mod_nsc_ComputeResidualProjection
   use Mod_nsc_InterpolateResidualProjection
   use Mod_nsc_ComputeGradientProjection
   use Mod_nsc_ArbitraryLagrangianEulerian
   use Mod_nsc_SubgridSpaceResidual
   use Mod_nsc_ComputeSubscales
   use Mod_nsc_ComputeNonLinearSubscales
   implicit none

contains

   !-------------------------------------------------------------------
   !SetPointers
   subroutine SetPointers
      implicit none

      !External Procedures
      procedure() :: NULLSUB

      integer(ip) :: kfl_nonlinear, nelty

      !Reset Procedure Composition
      call ResetProcedureComposition

      !Set All Pointers To NULLSUB (in Mod_nsc_BaseElmope)
      call SetPointersAndHooksToNULLSUB

      !Set the Specific Pointers for this subroutine
      ProcPointer_nsc_PreAssembly     => NULLSUB

      !------------------------------------------------------------------
      !Initialize the ProcedureFlags
      call SetPointersVariables(0)
      call SetPointersConvectionCoefficients(0)
      call SetPointersConvectionJacobianGradientCoefficients(0)
      call SetPointersDiffusionCoefficients(0)
      call SetPointersNonLinearDiffusionCoefficients(0)
      call SetPointersReactionCoefficients(0)
      call SetPointersExternalForces(0)
     !call SetPointersInterpolateGradients(0)
      call SetPointersComputeTaus(0)
      call SetPointersComputeTestf(0)

      call SetPointersInterpolateResidualProjection(0)
      call SetPointersComputeSubgridSpaceResidual(0)
      call SetPointersNonLinearDerivatives(0)
      call SetPointersArbitraryLagrangianEulerian(0)
      call SetPointersGetSubscales(0)

      !Also for the specific EndElmope Procedures
      call SetPointersComputeGpResidual(0)
      call SetPointersComputeResidualProjection(0)
      call SetPointersComputeGradientProjection(0)
      call SetPointersComputeNonLinearSubscales(0)
      !call SetPointersComputeTauSmoothing(0)
      call SetPointersComputeSubscales(0)

      !-----------------------------------------------------------------------
      !Now we set the required pointers

      !Non-linear elements
      call a%Mesh%IsNonLinearElement(ielem,kfl_nonlinear)
      if (kfl_nonlinear == 1) then
         call SetPointersNonLinearDerivatives(1)
      endif
      call SetPointersVariables(1)

      !Convection
      call SetPointersConvectionCoefficients(1)

       !Convection Jacobian Gradient Matrix Coefficients
      if (a%kfl_jacgr == 1) then
         call SetPointersConvectionJacobianGradientCoefficients(1)
      endif

      !Diffusion Exists
      if (a%kfl_visco == 1) then
          call SetPointersDiffusionCoefficients(1)
          !Non-linear elements
          if (kfl_nonlinear == 1) then
             call SetPointersNonLinearDiffusionCoefficients(1)
         endif
      endif

      !Reaction (transformed sources) Exists
      if (a%kfl_react == 1) then
         call SetPointersReactionCoefficients(1)
      endif

      !External forces
      call SetPointersExternalForces(1)

      !-----------------------------------------------------------
      !Shock Capturing
      if (a%kfl_shock /= 0 ) then
        if (a%kfl_shock == 1) then
             call SetPointersComputeGpResidual(1)
         else if (a%kfl_shock == 2) then
             call SetPointersComputeGradientProjection(1)
         else if (a%kfl_shock > 2) then
             call runend('Nsc_elmope: Other shock capturing methods not ready')
        endif
      endif
      !-----------------------------------------------------------
      !Residual Projection
      if (itask .eq. 'Endite' .and. a%kfl_repro >= 1) call SetPointersComputeResidualProjection(1)

      if ((itask .eq. 'Endste') .and. (a%npp_stepi(13)>0)) call SetPointersComputeGpResidual(1)

      !---------------------------------------------------------------
      !Dynamic, non-linear subscales
      !For Endste, only if tracking or transient, do not do it if non-linear
      if (itask .eq. 'Endste' .and. a%kfl_nolsg == 0 .and. (a%kfl_trasg == 1 .or. a%kfl_tacsg == 1)) then
         call SetPointersComputeSubscales(1)
      endif

      !For Endite, if non-linear subscales
      if (itask .eq. 'Endite' .and. a%kfl_nolsg == 1) then
         call SetPointersComputeNonLinearSubscales(1)
      endif

      !-----------------------------------------------------------
      !For adaptivity error indicators and subscales on the element boundaries
      if (itask .eq. 'Endste' .and. a%ErrorEstimatorTypeOfSubscales == 1) then
         call SetPointersComputeGradientProjection(1)
      endif

      !If more than one element type, then pointers should be reset when ielty changes!
      call a%Mesh%GetNelty(nelty)
      if (nelty > 1) ProcHook_nsc_OnIeltyChange => OnIeltyChange

      !Deallocate the procedure flags, so that they can be used in the next run
      !Initialize the ProcedureFlags
      call SetPointersVariables(100)
      call SetPointersConvectionCoefficients(100)
      call SetPointersConvectionJacobianGradientCoefficients(100)
      call SetPointersDiffusionCoefficients(100)
      call SetPointersNonLinearDiffusionCoefficients(100)
      call SetPointersReactionCoefficients(100)
      call SetPointersExternalForces(100)
      !call SetPointersInterpolateGradients(100)
      call SetPointersComputeTaus(100)
      call SetPointersComputeTestf(100)
      call SetPointersInterpolateResidualProjection(100)
      call SetPointersComputeSubgridSpaceResidual(100)
      call SetPointersNonLinearDerivatives(100)
      !call SetPointersTurbulenceModel(100)
      call SetPointersGetSubscales(100)
      call SetPointersArbitraryLagrangianEulerian(100)

      !Also for the specific EndElmope Procedures
      call SetPointersComputeGpResidual(100)
      call SetPointersComputeResidualProjection(100)
      call SetPointersComputeGradientProjection(100)
      call SetPointersComputeNonLinearSubscales(100)
      call SetPointersComputeSubscales(100)

   end subroutine

   !--------------------------------------------------------------------
   !Multiple type of elements
   subroutine OnIeltyChange
      implicit none

      if (e%ielty /= ielty0) then
         call SetPointers
         ielty0 = e%ielty
      endif
   end subroutine

end module

subroutine nsc_EndElmope_im(NSCompImplicitProblem,task)
   use Mod_nsc_BaseElmope
   use Mod_nsc_EndElmope_im
   use Mod_NSCompressibleImplicit
   implicit none
   class(NSCompressibleImplicitProblem), target :: NSCompImplicitProblem
   character(6) :: task
   integer(ip) :: aux_logic

   a=>NSCompImplicitProblem
   itask = task

   !Things to be done if endite
   if (itask .eq. 'Endite') then
      !Return if there is nothing to be done
      aux_logic = 0

      if (a%kfl_repro >= 1) aux_logic = 1
      if (a%kfl_nolsg == 1) aux_logic = 1
      if (a%kfl_shock /= 0 ) aux_logic = 1

      if (aux_logic == 0) return

   !Things to be done if endste
   elseif (itask .eq. 'Endste') then
      !Return if there is nothing to be don
      aux_logic = 0
      if (a%kfl_trasg /= 0 .and. a%kfl_nolsg == 0) aux_logic = aux_logic + 1
      if (a%kfl_tacsg /= 0 .and. a%kfl_nolsg == 0) aux_logic = aux_logic + 1

      if (aux_logic == 0) return

      !If only dissipation for postprocessing and not a postprocessing step return
!      if (a%kfl_dispa == 1 .and. aux_logic == 1) then
!         if (a%istep < a%npp_inits) return
!         if (a%npp_stepi(15)>0 ) then
!            if(mod(a%istep,a%npp_stepi(15))/=0) then
!               return
!            endif
!         endif
!      endif
   endif

   !We force closed rule for smoothing
   call a%Mesh%ElementAlloc(e,a%Memor,a%EndLoopQuadrature,'nsc_EnditeElmope')

   !Set the Pointers for execution
   ielem = 1
   call SetPointers

   !Hook
   call ProcHook_nsc_PreLoop

   call a%Memor%alloc(a%ndofn,e%mnode,a%ndofn,e%mnode,elmat,'elmat','nsc_EndElmope')

   !AllocateArrays in BaseElmope
   call AllocateBaseElmopeArrays

   !Hook
   call ProcHook_nsc_Initializations

   !do itest = 1,100
   call a%Mesh%GetNelem(nelem)
   elements : do ielem = 1,nelem
      !Load Element
      call a%Mesh%ElementLoad(ielem,e)

      !Hook
      call ProcHook_nsc_OnIeltyChange

      !Elmats to Zero
      call ProcHook_nsc_ElmatsToZero
      elmat=0.0_rp

      !Gathers
      call GatherBase
      !Hook
      call ProcHook_nsc_Gathers

      !Cartesian derivatives and Jacobian at center of gravity
      call e%elmdcg

      !Element length at center of gravity
      call e%elmlen

      call CalculateElmchl

      !Hook
      call ProcHook_nsc_PreGauss

      !Gauss Point Loop
      gauss_points : do igaus=1,e%pgaus
         e%igaus = igaus

         !Hook
         call ProcHook_nsc_InGauss

         dvol = e%weigp(e%igaus)*e%detjm

         !Interpolate
         call InterpolateBase
         !Hook
         call ProcHook_nsc_Interpolates

         !Physical Properties
         !Hook
         call ProcHook_nsc_PhysicalProp

         !Linear transport variables
         call ComputeLinearVariables

         !Nonlinear transport variables
         call ProcPointer_nsc_ComputeVariables

         !Sound speed
         call nsc_ComputeSoundSpeed(accph,accvh,gptem,gpspd)

         !Compute gradients
         call ComputeBaseGradients

         !Advection matrix
         call ProcPointer_nsc_ComputeConvectionCoefficients

         !Diffusion matrix
         call ProcPointer_nsc_ComputeDiffusionCoefficients

         !Other transport matrices
         call ProcPointer_nsc_ComputeTransportCoefficients

         !Compute the stability parameters, only if required
         call ProcHook_nsc_ComputeTaus

         !Compute the stabilization test function, only if required
         call ProcHook_nsc_ComputeTestf

         !Compute Elext, Temporal Derivatives
         elexd = 0.0_rp
         elexm = 0.0_rp
         elexe = 0.0_rp
         !Compute vector of external forces
         call ProcPointer_nsc_ExternalForces

         !Time integration
         eltemd = 0.0_rp
         eltemm = 0.0_rp
         elteme = 0.0_rp
         call nsc_TimeIntegrationToElTemp(e,1_ip,Integrator,a%dtinv,gpden,eltemd)
         call nsc_TimeIntegrationToElTemp(e,e%ndime,Integrator,a%dtinv,gpmom,eltemm)
         call nsc_TimeIntegrationToElTemp(e,1_ip,Integrator,a%dtinv,gpene,elteme)

         !InGaussElmats
         !Hook
         call ProcHook_nsc_InGaussElmats

         !InGaussElmats Assembly
         call ProcHook_nsc_InGaussElmatsAssembly

      enddo gauss_points

      !PreDirichlet
      call ProcHook_nsc_PreDirichlet

      !Assembly Endite
      !Hook
      call ProcHook_nsc_AssemblyEndite

   enddo elements
   !enddo

   !Hook
   call ProcHook_nsc_Finalizations

   !Deallocate arrays in BaseElmope
   call DeallocateBaseElmopeArrays

   call a%Memor%dealloc(a%ndofn,e%mnode,a%ndofn,e%mnode,elmat,'elmat','nsc_EndElmope')

   !Operations to be done after the Elemental Loop
   !Hook
   call ProcHook_nsc_PostLoop

   !Element
   call a%Mesh%ElementDeAlloc(e,a%Memor,a%EndLoopQuadrature,'nsc_EnditeElmope')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
C*26-NOV-90 11:35:29 MDS$SOURCE:[KKLARE]ZBLAS.FOR;5
! Subprograms based on "LINPACK User's Guide" by J.J. Dongarra et al.
! The Basic Linear Algebra Subprograms portion.
!
! root    prefix-suffix of name     Function
! -AMAX   IS- ID-     IC-     IZ-   index of max |x|  or  max  |x|+|y|
! -ASUM   S-  D-      SC-     DZ-   w := sum |x|  or  sum |x|+|y|
! -AXPY   S-  D-      C-      Z-    w := a*x+y
! -COPY   S-  D-      C-      Z-    y := x
! -DOT-   S-  D-  C-U C-C Z-U Z-C   w:=sum x*y, for -C  w:=conjg(x)*y
! -NRM2   S-  D-      SC-     DZ-   w := sqrt(sum |x|**2)
! -ROT    S-  D-      CS-     ZD-   apply rotation
! -ROTG   S-  D-      C-      Z-    setup Givens rotation
! -SCAL   S-  D-  CS- C-  ZD- Z-    x := ax
! -SWAP   S-  D-      C-      Z-    y :=: x
! -ROTG used by -CHEX and -CHUD.
! S-ROTG and CS-ROT used by C-SVDC etc.
!
! S- are single-precision.
! D- are double-precision.
! C- are single-precision complex.
! Z- are double-precision complex.
!
! N is the active number of elements.
! X and Y are vectors stepped in increments INCX and INCY.
! Dimension 1 is used for these arrays.
! LINPACK uses INCX=INCY=1 and special code for this case.
! Negative INCX or INCY can mean reversed order.
! A, B, C, and S are scalars.
!
!+ index of absolute-value maximum
      Integer Function ICAMAX(n,x,incx)
      Integer i,incx,j,n
      REAL xmax
      COMPLEX x(n*incx)
      If (n.le.0) Then
         ICAMAX=0
      ElseIf (n.eq.1) Then
         ICAMAX=1
      Else
         ICAMAX=1
         xmax=abs(REAL(x(1)))+abs(AIMAG(x(1)))
         If (incx.ne.1) Then
            j=1
            Do i=2,n
               j=j+incx
               If (abs(REAL(x(j)))+abs(AIMAG(x(j))).gt.xmax) Then
                xmax=abs(REAL(x(j)))+abs(AIMAG(x(j)))
                ICAMAX=i
               EndIf
            EndDo
         Else
            Do i=2,n
               If (abs(REAL(x(i)))+abs(AIMAG(x(i))).gt.xmax) Then
                  xmax=abs(REAL(x(i)))+abs(AIMAG(x(i)))
                  ICAMAX=i
               EndIf
            EndDo
         EndIf
      EndIf
      End
!+ Euclidean norm of vector
      REAL Function SCNRM2(n,x,incx)
      REAL SCASUM
      Integer i,incx,n
      REAL cutlo,cuthi,test
      COMPLEX x(1)
! machine precision      smallest-exp      largest-exp
! Prime            S 23      -128         127
! Prime            D 47      -32896(-32768-128) 32639(32767-128)
! Vax            S 24      -128         127
! Vax            D 56      -128         127
! DEC KL10      S 27      -128         127
! DEC KL10      D 62      -128         127 //worst case//
! DEC KI10      D 54      -128         127
! CDC-6600      S 48      -975(-1023+48)      1060(+1022+48) (+1023=infinity)
! CDC-6600      D 96      -975         1060
! The smallest number is 2**(smallest)*.5.
! The largest number is 2**(largest)*(1-2/2**precision).
! The expression eps+1 .gt. 1 is true for eps=2**(1-precision)
! unrounded or eps=2**(-precision) when rounded.
! This does not work on Prime where S-accumulator is 31 bits.
! For un-squared number halve the exponent.
! So we scale for maximum above largest/N (/2N for complex)
! And for maximum below CUTLO**2=smallest/eps.
! It is assumed that underflows are efficiently handled as zero.
! Round up CUTLO, down CUTHI.
!primecs          data cutlo,cuthi/1.570e-16,1.304e+19/
!primecd          data cutlo,cuthi/5.406e-4945,4.560e+4912/
!kl-10cs          data cutlo,cuthi/6.281e-16,1.304e+19/
!kl-10cd          data cutlo,cuthi/1.164e-10,1.304e+19/
!vaxcs          data cutlo,cuthi/8.881e-16,1.304e+19/
!vaxcd          data cutlo,cuthi/2.058e-11,1.304e+19/
      data cutlo,cuthi/8.232e-11,1.304e+19/
      SCNRM2=0
      If (n.le.0) Return
      test=SCASUM(n,x,incx)
      If (test.le.0) Return
!CIF
      If (test.le.cuthi.and.test*(n+n).ge.cutlo) Then
!CEND      If (test.le.cuthi.and.test*n.ge.cutlo) Then
         Do j=1,n*incx,incx
            SCNRM2=SCNRM2+REAL(x(j))**2
     1                  +AIMAG(x(j))**2
         EndDo
         SCNRM2=sqrt(SCNRM2)
      Else
         Do i=1,n*incx,incx
            SCNRM2=SCNRM2+(REAL(x(i))/test)**2
     1                  +(AIMAG(x(i))/test)**2
         EndDo
         SCNRM2=test*sqrt(SCNRM2)
      Endif
      End
!+ absolute values summed
      REAL Function SCASUM(n,x,incx)
      Integer i,incx,m,n
      COMPLEX x(1)
      SCASUM=0
      If (n.le.0) Then
      ElseIf (incx.ne.1) Then
         Do i=1,n*incx,incx
            SCASUM=SCASUM+abs(REAL(x(i)))+abs(AIMAG(x(i)))
         EndDo
      Else
         m=mod(n,6)
         Do i=1,m
            SCASUM=SCASUM+abs(REAL(x(i)))+abs(AIMAG(x(i)))
         EndDo
         Do i=m+1,n,6
            SCASUM=SCASUM+abs(REAL(x(i)))+abs(AIMAG(x(i)))
     1+abs(REAL(x(i+1)))+abs(AIMAG(x(i+1)))
     2+abs(REAL(x(i+2)))+abs(AIMAG(x(i+2)))
     3+abs(REAL(x(i+3)))+abs(AIMAG(x(i+3)))
     4+abs(REAL(x(i+4)))+abs(AIMAG(x(i+4)))
     5+abs(REAL(x(i+5)))+abs(AIMAG(x(i+5)))
         EndDo
      EndIf
      End
!+ y=a*x+y  scale and add, backwards if negative increment
      SUBROUTINE CAXPY(n,a,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      COMPLEX x(1),y(1),a
      If (n.le.0.or.a.eq.0) Then
      ElseIf (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            y(k)=y(k)+x(j)*a
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,4)
         Do i=1,m
            y(i)=y(i)+x(i)*a
         EndDo
         Do i=m+1,n,4
            y(i)=y(i)+x(i)*a
            y(i+1)=y(i+1)+x(i+1)*a
            y(i+2)=y(i+2)+x(i+2)*a
            y(i+3)=y(i+3)+x(i+3)*a
         EndDo
      EndIf
      End
!+ copy array, backwards if negative increment
      SUBROUTINE CCOPY(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      COMPLEX x(1),y(1)
      If (n.le.0) Then
      ElseIf (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            y(k)=x(j)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,7)
         Do i=1,m
            y(i)=x(i)
         EndDo
         Do i=m+1,n,7
            y(i)=x(i)
            y(i+1)=x(i+1)
            y(i+2)=x(i+2)
            y(i+3)=x(i+3)
            y(i+4)=x(i+4)
            y(i+5)=x(i+5)
            y(i+6)=x(i+6)
         EndDo
      EndIf
      End
!+ dot product
      COMPLEX Function CDOTU(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      COMPLEX x(1),y(1)
      CDOTU=0
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            CDOTU=CDOTU+x(j)*y(k)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            CDOTU=CDOTU+x(i)*y(i)
         Enddo
         Do i=m+1,n,5
            CDOTU=CDOTU+x(i)*y(i)+y(i+1)*x(i+1)+
     1     x(i+2)*y(i+2)+y(i+3)*x(i+3)+x(i+4)*y(i+4)
         EndDo
      EndIf
      End
!CIF
!+ dot product
      COMPLEX Function CDOTC(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      COMPLEX x(1),y(1)
      CDOTC=0
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            CDOTC=CDOTC+conjg(x(j))*y(k)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            CDOTC=CDOTC+conjg(x(i))*y(i)
         Enddo
         Do i=m+1,n,5
            CDOTC=CDOTC+conjg(x(i))*y(i)
     1+conjg(x(i+1))*y(i+1)
     2+conjg(x(i+2))*y(i+2)
     3+conjg(x(i+3))*y(i+3)
     4+conjg(x(i+4))*y(i+4)
         EndDo
      Endif
      End
!CEND
!+ apply a plane rotation, backwards if negative increment
      Subroutine CSROT(n,x,incx,y,incy,c,s)
      Integer i,incx,incy,j,k,n
      COMPLEX x(1),y(1),temp
      REAL c,s
      If (n.le.0.or.s.eq.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            temp=
     1CMPLX(REAL(x(j))*c+REAL(y(k))*s,AIMAG(x(j))*c+AIMAG(y(k))*s)
            y(k)=
     1CMPLX(-REAL(x(j))*s+REAL(y(k))*c,-AIMAG(x(j))*s+AIMAG(y(k))*c)
            x(j)=temp
            j=j+incx
            k=k+incy
         EndDo
      Else
         Do i=1,n
            temp=
     1 CMPLX(REAL(x(i))*c+REAL(y(i))*s,AIMAG(x(i))*c+AIMAG(y(i))*s)
            y(i)=
     1 CMPLX(-REAL(x(i))*s+REAL(y(i))*c,-AIMAG(x(i))*s+AIMAG(y(i))*c)
            x(i)=temp
         EndDo
      Endif
      End
!+ construct Givens plane rotation
      SUBROUTINE CROTG(a,b,c,s)
      REAL c,aa,bb
      COMPLEX a,b,s,r
      bb=abs(REAL(b))+abs(AIMAG(b))
      If (bb.gt.0) Then
         aa=abs(REAL(a))+abs(AIMAG(a))
         If (aa.gt.bb) Then
            r=a*sqrt(1+(bb/aa)**2)
            c=real(a/r)
            s=b/r
            b=s
            a=r
         Else
            r=b*sqrt(1+(aa/bb)**2)
            c=real(a/r)
            s=b/r
            a=r
            b=1
            If (c.ne.0) b=1/c
         Endif
      Else
            c=1
            s=0
      Endif
      End
!CIF
      Subroutine CSSCAL(n,a,x,incx)
      Integer i,incx,m,n
      REAL a
      COMPLEX x(1)
      If (n.le.0) Return
      If (incx.ne.1) Then
         Do i=1,n*incx,incx
            x(i)=CMPLX(REAL(x(i))*a,AIMAG(x(i))*a)
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            x(i)=CMPLX(REAL(x(i))*a,AIMAG(x(i))*a)
         EndDo
         Do i=m+1,n,5
            x(i)=CMPLX(REAL(x(i))*a,AIMAG(x(i))*a)
            x(i+1)=CMPLX(REAL(x(i+1))*a,AIMAG(x(i+1))*a)
            x(i+2)=CMPLX(REAL(x(i+2))*a,AIMAG(x(i+2))*a)
            x(i+3)=CMPLX(REAL(x(i+3))*a,AIMAG(x(i+3))*a)
            x(i+4)=CMPLX(REAL(x(i+4))*a,AIMAG(x(i+4))*a)
         EndDo
      Endif
      End
!CEND
!+ x=a*x  scale
      SUBROUTINE CSCAL(n,a,x,incx)
      Integer i,incx,m,n
      COMPLEX x(1),a
      If (n.le.0) Return
      If (incx.ne.1) Then
         Do i=1,n*incx,incx
            x(i)=x(i)*a
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            x(i)=x(i)*a
         EndDo
         Do i=m+1,n,5
            x(i)=x(i)*a
            x(i+1)=x(i+1)*a
            x(i+2)=x(i+2)*a
            x(i+3)=x(i+3)*a
            x(i+4)=x(i+4)*a
         EndDo
      EndIf
      End
!+ swap arrays, backwards if negative increment
      SUBROUTINE CSWAP(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      COMPLEX x(1),y(1),temp
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            temp=x(j)
            x(j)=y(k)
            y(k)=temp
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,3)
         Do i=1,m
            temp=x(i)
            x(i)=y(i)
            y(i)=temp
         EndDo
         Do i=m+1,n,3
            temp=x(i)
            x(i)=y(i)
            y(i)=temp
            temp=x(i+1)
            x(i+1)=y(i+1)
            y(i+1)=temp
            temp=x(i+2)
            x(i+2)=y(i+2)
            y(i+2)=temp
         EndDo
      EndIf
      End

ext: f

max_stars_repo_name: fluffynukeit/mdsplus

max_issues_repo_name: fluffynukeit/mdsplus

max_forks_repo_name: fluffynukeit/mdsplus


module parameters
use kind_parameters
 !! Constant Parameters

 real(kind=dp), parameter:: pi = 3.1415926535897932384626433832795028841971693993751058209d0

 !! Iteration variables
 integer:: itermax, counter, maxcase, qval
 integer*4:: icv, jcv, ipmax, jpmax, N, isol, iteration
 integer, allocatable:: iter_num(:)

 !! Other parameters
 real(kind=dp):: Lx, Ly, reynolds_number, inlet_velocity, lid_velocity, hydraulic_dia, totalmb, mass_ratio, v_ref, rho_ref, time_ref, mu_ref, bond_number
 real(kind=dp):: beta, tambient, Rayleigh, alpha, gravity

 !! For specifying boundary conditions
 real(kind=dp):: u_west, u_east, u_south, u_north, v_west, v_east, v_south, v_north
 real(kind=dp):: t_west, t_east, t_south, t_north

!! Property Flag
 integer*4:: prop_flag

!! Under-relaxation parameters, residual variables

 real(kind=dp):: alphau, alphav, alphap, alphat, alphapc, alpha_mg
 real(kind=dp):: epsilonu, epsilonv, epsilonp, epsilont

! Residual Variables
 real(kind=dp):: residualu, residualv, residualp, residualt
 integer*4 :: is_converged

! Monitor convergence and output steps variables

integer*4 :: output_step, monitor_step
logical :: write_VTK

! Debug variables

integer*4 :: debug, parallel_debug

!! CPU Time computing variables

 real(kind=dp):: tstart, tend

  !! Print files variables
 character*50:: filename, string, input_file, prog_name,  grid_type, mg_handle
 logical:: is_input_data, variable_property

!! Unsteady variables

 real(kind=dp):: time_step, max_time, time, tlocal1, tlocal2
 integer:: tcounter

! Constants for Runge-Kutta 2nd order implicit method
real(kind=dp):: A21 = 0.0d0
real(kind=dp):: A22 = 1.0d0

! Parameters for level set method

real(kind=dp):: CFL, epsilon, delta_phi, surface_tension
real(kind=dp) :: rho_g, rho_l, mu_g, mu_l, kcond_g, kcond_l, cp_g, cp_l
real(kind=dp):: radius, center_x, center_y, L_perimeter, global_error, volume_error, error_sum, H_phi_exact, H_phi_numerical, vol_exact, vol_sum, vol_error
real(kind=dp):: alphax1, alphax2, alphax3, alphay1, alphay2, alphay3, phi_xx, phi_xy, phi_yy
real(kind=dp):: tau_redistance, phi_heaviside, phi_delta, gamma_ls, beta_ls
integer:: iter_redistance

end module parameters


text:
      SUBROUTINE PHO_DQMASS(I,J,K,L,Am82,Am102)
C**********************************************************************
C
C     determine minimal masses corresponding to the input flavours
C     (diquark a-diquark string system)
C
C     input: I,J,K,L   quark flavours (PDG convention)
C
C     output: AM82     mass of two octett baryons
C             AM102    mass of two decuplett baryons
C
C**********************************************************************

      IMPLICIT NONE

      SAVE

      INTEGER I , J , K , L
      DOUBLE PRECISION Am82 , Am102

C  input/output channels
      INCLUDE 'inc/poinou'
C  event debugging information
      INCLUDE 'inc/podebg'
C  general particle data
      INCLUDE 'inc/popar2'

C  local variables
      INTEGER ii , jj , kk , ll

      ii = ABS(I)
      kk = ABS(K)
      jj = ABS(J)
      ll = ABS(L)

      Am82 = XM_bb82_list(ii,jj,kk,ll)
      Am102 = XM_bb102_list(ii,jj,kk,ll)

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
************************************************************************
* (Static) Grid adaption, multpliple-level version
*
* Direct version without using old SETLEV-method for switching
* levels. Works on the parameters given by the initialization
* subroutine of the multiple-grid version of the grid adaption
* (INGAST/INGAMG/INGASM).
*
* This subroutine performs grid-adaption according to information
* in the triangulation structures. The grid adaption takes place
* on the levels configured by INGASM, so the caller has not to
* take any influence concerning levels of adaption etc.
*
* By specifying another NLMAX than the standard one, the caller
* can specify the level of grid adaption, on which the parameters
* configured by INGASM (i.e. in the standard implementation, the
* parameters of the .DAT file) act.
*
* In:
*   TRIAS  - array [1..SZTRIA,1..NLMAX] of integer
*            Triangulation structures on all levels
*   NLMIN  - Minimum level in TRIAS
*   NLMAX  - Maximum level in TRIAS
*   IGEOM  - array [1..*] of integer
*   DGEOM  - array [1..*] of double
*            Integer- and double-precision parameter blocks with
*            geometry information. Passed to boundary
*            routines. Not used in this routine.
*   IDATA  - array [1..*] of integer
*   DDATA  - array [1..*] of double
*            User defined integer- and double-precision parameter
*            blocks. Passed to callback routines like monitor
*            function. Not used in this routine.
************************************************************************

      SUBROUTINE XSMMGW (TRIAS,NLMIN,NLMAX,IGEOM,DGEOM,IDATA,DDATA)

      IMPLICIT NONE

      INCLUDE 'cmem.inc'
      INCLUDE 'cout.inc'
      INCLUDE 'cfileout.inc'

      INCLUDE 'cbasicmg.inc'

      INCLUDE 'stria.inc'

C     parameters

      INTEGER TRIAS(SZTRIA,NNLEV),NLMIN,NLMAX,IGEOM(*),IDATA(*)
      DOUBLE PRECISION DGEOM(*),DDATA(*)

C     Variables for the static method:

      INCLUDE 'sgridadaptgeneralparams.inc'
      INCLUDE 'sgridadaptsolverparams.inc'
      INCLUDE 'sgridadaptstaticparams.inc'
      INCLUDE 'sgridadaptmgstaticparams.inc'

      INTEGER GASPRI(SZGACI+SZGLSI+SZGASI),GASPMI(SZGAMI)
      DOUBLE PRECISION GASPRD(SZGACD+SZGLSD+SZGASD)

C     local variables

      INTEGER II,ILEV,ILVSLV,ICONF
      DOUBLE PRECISION TGES
      DOUBLE PRECISION TGTOT,TGLSPR,TGLS,TGODE,TGGPR,TGGSM,TGMSM
      DOUBLE PRECISION TGGAPR,TGGRC
      INTEGER NLSYS ,NLSITE,NCODE ,NMFEVL
      CHARACTER CSTR*255

C     Initialize the parameter arrays with the variables
C     in the COMMON blocks from the DAT file.

      CALL INGAST (GASPRI, GASPRD)

C     Initialize the parameter block of the multiple-grid adaption
C     method

      CALL INGAMG (GASPMI)

C     Clear timing information for output

      TGTOT  = 0.0
      TGLSPR = 0.0
      TGLS   = 0.0
      TGODE  = 0.0
      TGGPR  = 0.0
      TGGSM  = 0.0
      TGMSM  = 0.0
      TGGAPR = 0.0
      TGGRC  = 0.0

      NLSYS  = 0
      NLSITE = 0
      NCODE  = 0
      NMFEVL = 0

C     Loop through all adaption configurations

      DO ICONF = 1,GASPMI(OIMLGAS)

C       Switch to the configuration.
C       Obtain ILEV and ILVSLV, i.e. the information on which
C       levels to adapt and to solve the LGS that occur:

        CALL INGASM (ICONF,GASPRI,GASPMI,ILEV,ILVSLV)

C       Switch to the level where we should adapt:

        II=1
        IF (ILEV.LE.0) THEN
          ILEV = NLMAX+ILEV
        END IF
        IF (ILEV.GT.NLMAX) ILEV = NLMAX

        IF (ILVSLV.LE.0) THEN
          ILVSLV = NLMAX+ILVSLV
        END IF
        IF (ILVSLV.GT.NLMAX) ILVSLV = NLMAX
        IF (ILVSLV.LT.NLMIN) ILVSLV = NLMIN

        WRITE(CSTR,'(A)')
     *              '=============================================='
        CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
        WRITE(CSTR,'(A,I5)') 'Grid deformation configuration: ',ICONF
        CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
        WRITE(CSTR,'(A,I5)') 'Adaptation on level:            ',ILEV
        CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
        WRITE(CSTR,*)
        CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)

C       No grid adaption if we drop below the minimum level.

        IF (ILEV.GE.NLMIN) THEN

C         Adapt the grid and recreate missing information on all levels

          CALL XSMMGS (TRIAS,NLMIN,NLMAX,ILEV,ILVSLV,
     *                 GASPRI,GASPRD,IGEOM,DGEOM,IDATA,DDATA)

          TGTOT  = TGTOT  + GASPRD(OTGTOT )
          TGLSPR = TGLSPR + GASPRD(OTGLSPR)
          TGLS   = TGLS   + GASPRD(OTGLS  )
          TGODE  = TGODE  + GASPRD(OTGODE )
          TGGPR  = TGGPR  + GASPRD(OTGGPR )
          TGGSM  = TGGSM  + GASPRD(OTGGSM )
          TGMSM  = TGMSM  + GASPRD(OTGMSM )
          TGGAPR = TGGAPR + GASPRD(OTGGAPR)
          TGGRC  = TGGRC  + GASPRD(OTGGRC )

          NLSYS  = NLSYS  + GASPRI(ONLSYS )
          NLSITE = NLSITE + GASPRI(ONLSITE)
          NCODE  = NCODE  + GASPRI(ONCODE )
          NMFEVL = NMFEVL + GASPRI(ONMFEVL)

        ELSE

          WRITE(CSTR,'(A)') 'Grid adaption skipped on this level.'
          CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
          WRITE(CSTR,'(I3,A,I3)') ILEV,' = ILEV < NLMIN = ',NLMIN
          CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)

        END IF

      END DO

C     Finally some output, finish.

      TGES=TGTOT/100D0
      WRITE(CSTR,'(A)')
     *            '=============================================='
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A)') 'Grid deformation statistics.'
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,*)
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Absolute / Rel. time    : ',
     *             TGTOT,100D0
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Grid prolongation/restr.: ',
     *             TGGPR,TGGPR/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Preparation of grid ad. : ',
     *             TGGAPR,TGGAPR/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Setting up linear sys.  : ',
     *             TGLSPR,TGLSPR/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Solving Linear Systems  : ',
     *             TGLS,TGLS/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Moving grid             : ',
     *             TGODE,TGODE/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Grid smoothing:         : ',
     *             TGGSM,TGGSM/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Monitor fct. smoothing  : ',
     *             TGMSM,TGMSM/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,F16.6,F14.6)') 'Gradient calculation    : ',
     *             TGGRC,TGGRC/TGES
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,*)
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,I12)') '#linear systems         : ',
     *             NLSYS
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,I12)') '#iterations in lin. slv.: ',
     *             NLSITE
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,I12)') '#ODE''s solved           : ',
     *             NCODE
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,'(A,I12)') '#monitor function eval. : ',
     *             NMFEVL
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE(CSTR,*)
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)

      WRITE (CSTR,'(A)') 'Grid deformation completed.'
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE (CSTR,'(A)')
     *             '=============================================='
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)
      WRITE (CSTR,*)
      CALL CNOUTS (MT,MTERM,MFILE,.TRUE.,CSTR)

      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      subroutine initav(array,nmax,value)
c
c     This subroutine initializes the real*8 array "array" to the value
c     of "value" over the first nmax positions. Normally, nmax would
c     be the dimension of the 1D "array". However, nmax could be less
c     than the true dimension. Also, array could actually have more
c     than one dimension. For example, if the array is really "a",
c     having dimensions of i1 and i2, this subroutine could be used by
c     calling it in the following manner: call initar(a,(i1*i2),value).
c     Use this subroutine with caution if the two arrays do not have
c     identical dimensioning or nmax is not the product of the true
c     (declared) dimensions.
c
c     To initialize a real*8 array to a zero value, use EQLIBU/initaz.f
c     instead.
c
c     This subroutine is called by:
c
c       Any
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       array  = the array
c       nmax   = dimension or pseudo-dimension of array
c       value  = the initialization value
c
c     Output:
c
c       array  = array, with the first nmax positions set to the
c                initialization value
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nmax
c
      real*8 array(nmax),value
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,ileft
c
      real*8 vx
c
c-----------------------------------------------------------------------
c
c     Note the assignment of the value to a local variable.
c     Note that the loop is unrolled.
c
      vx = value
      ileft = (nmax/8)*8
c
      do i = 1,ileft,8
        array(i) = vx
        array(i + 1) = vx
        array(i + 2) = vx
        array(i + 3) = vx
        array(i + 4) = vx
        array(i + 5) = vx
        array(i + 6) = vx
        array(i + 7) = vx
      enddo
c
      do i = ileft + 1,nmax
        array(i) = vx
      enddo
c
      end


text:
module parameters

    ! Sizes for arrays

    integer :: MaxBlades, MaxSegPerBlade, MaxSegEndPerBlade, MaxSegEnds, MaxSeg     ! Maximums for blade geometry arrays
    integer :: MaxStruts                                                            ! Maximums for strut arrays
    integer :: MaxAirfoilSect, MaxReVals, MaxAOAVals                ! Maximums for airfoil coeff data arrays
    integer :: MaxRevs                              ! Max revolutions
    integer :: MaxTimeStepPerRev                            ! Max time increments in a revolution
    integer :: MaxTimeSteps                                                         ! Max total time increments (MaxRevs * MaxTimeStepPerRev)
    integer :: MaxWakeNodes                             ! Max number of points in wake for an element (MaxRevs * MaxTimeStepPerRev)
    integer :: MaxNLIters                               ! Max iterations to perform in converging non linear component of blade element system

end module parameters

max_stars_repo_name: ebranlard/CACTUS

max_issues_repo_name: ebranlard/CACTUS

max_forks_repo_name: ebranlard/CACTUS


text:
!==========================================================================
elemental function gsw_enthalpy_sso_0 (p)
!==========================================================================
!  This function calculates enthalpy at the Standard Ocean Salinity, SSO,
!  and at a Conservative Temperature of zero degrees C, as a function of
!  pressure, p, in dbar, using a streamlined version of the
!  computationally-efficient expression for specific volume, that is, a
!  streamlined version of the code "gsw_enthalpy(SA,CT,p)".
!==========================================================================

use gsw_mod_teos10_constants, only : db2pa

use gsw_mod_specvol_coefficients

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: p

real (r8) :: gsw_enthalpy_sso_0

real (r8) :: dynamic_enthalpy_sso_0_p, z

z = p*1e-4_r8

dynamic_enthalpy_sso_0_p = &
      z*( 9.726613854843870e-4_r8 + z*(-2.252956605630465e-5_r8 &
    + z*( 2.376909655387404e-6_r8 + z*(-1.664294869986011e-7_r8 &
    + z*(-5.988108894465758e-9_r8 + z*(h006 + h007*z))))))

gsw_enthalpy_sso_0 = dynamic_enthalpy_sso_0_p*db2pa*1e4_r8

return
end function

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

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
 SUBROUTINE POLATES6(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
                     IGDTNUMO,IGDTMPLO,IGDTLENO, &
                     MI,MO,KM,IBI,LI,GI, &
                     NO,RLAT,RLON,IBO,LO,GO,IRET)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! $Revision: 74685 $
!
! SUBPROGRAM:  POLATES6   INTERPOLATE SCALAR FIELDS (BUDGET)
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM PERFORMS BUDGET INTERPOLATION
!           FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
!           THE ALGORITHM SIMPLY COMPUTES (WEIGHTED) AVERAGES
!           OF NEIGHBOR POINTS ARRANGED IN A SQUARE BOX
!           CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING
!           NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS.
!           OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS
!           FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2
!           (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED;
!           FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS
!           POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1))
!           WHICH DEFAULTS TO ALL 1 (IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!           ANOTHER OPTION IS THE MINIMUM PERCENTAGE FOR MASK,
!           I.E. PERCENT VALID INPUT DATA REQUIRED TO MAKE OUTPUT DATA,
!           (IPOPT(3+IPOPT(1)) WHICH DEFAULTS TO 50 (IF -1).
!           ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
!           THE CODE RECOGNIZES THE FOLLOWING PROJECTIONS, WHERE
!           "IGDTNUMI/O" IS THE GRIB 2 GRID DEFINTION TEMPLATE NUMBER
!           FOR THE INPUT AND OUTPUT GRIDS, RESPECTIVELY:
!             (IGDTNUMI/O=00) EQUIDISTANT CYLINDRICAL
!             (IGDTNUMI/O=01) ROTATED EQUIDISTANT CYLINDRICAL. "E" AND
!                             NON-"E" STAGGERED
!             (IGDTNUMI/O=10) MERCATOR CYLINDRICAL
!             (IGDTNUMI/O=20) POLAR STEREOGRAPHIC AZIMUTHAL
!             (IGDTNUMI/O=30) LAMBERT CONFORMAL CONICAL
!             (IGDTNUMI/O=40) GAUSSIAN CYLINDRICAL
!           AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
!           AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
!           INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
!           OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
!           EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
!           THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
!
! PROGRAM HISTORY LOG:
!   96-04-10  IREDELL
!   96-10-04  IREDELL  NEIGHBOR POINTS NOT BILINEAR INTERPOLATION
! 1999-04-08  IREDELL  SPLIT IJKGDS INTO TWO PIECES
! 2001-06-18  IREDELL  INCLUDE MINIMUM MASK PERCENTAGE OPTION
! 2015-01-27  GAYNO    REPLACE CALLS TO GDSWIZ WITH NEW MERGED
!                      VERSION OF GDSWZD.
! 2015-07-13  GAYNO    CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAYS
!                      WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAYS.
!
! USAGE:    CALL POLATES6(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
!                         IGDTNUMO,IGDTMPLO,IGDTLENO, &
!                         MI,MO,KM,IBI,LI,GI, &
!                         NO,RLAT,RLON,IBO,LO,GO,IRET)
!
!   INPUT ARGUMENT LIST:
!     IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
!                IPOPT(1) IS NUMBER OF RADIUS POINTS
!                (DEFAULTS TO 2 IF IPOPT(1)=-1);
!                IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
!                (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!                IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
!                (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
!     IGDTNUMI - INTEGER GRID DEFINITION TEMPLATE NUMBER - INPUT GRID.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE:
!                  00 - EQUIDISTANT CYLINDRICAL
!                  01 - ROTATED EQUIDISTANT CYLINDRICAL.  "E"
!                       AND NON-"E" STAGGERED
!                  10 - MERCATOR CYCLINDRICAL
!                  20 - POLAR STEREOGRAPHIC AZIMUTHAL
!                  30 - LAMBERT CONFORMAL CONICAL
!                  40 - GAUSSIAN EQUIDISTANT CYCLINDRICAL
!     IGDTMPLI - INTEGER (IGDTLENI) GRID DEFINITION TEMPLATE ARRAY -
!                INPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
!                OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                (SECTION 3 INFO).  SEE COMMENTS IN ROUTINE
!                IPOLATES FOR COMPLETE DEFINITION.
!     IGDTLENI - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
!                TEMPLATE ARRAY - INPUT GRID.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     IGDTNUMO - INTEGER GRID DEFINITION TEMPLATE NUMBER - OUTPUT GRID.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                SAME DEFINITION AS "IGDTNUMI".
!     IGDTMPLO - INTEGER (IGDTLENO) GRID DEFINITION TEMPLATE ARRAY -
!                OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
!                OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                (SECTION 3 INFO).  SEE COMMENTS IN ROUTINE
!                IPOLATES FOR COMPLETE DEFINITION.
!     IGDTLENO - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
!                TEMPLATE ARRAY - OUTPUT GRID.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF INPUT GRID FIELDS IF KM=1
!     MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
!     KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
!     IBI      - INTEGER (KM) INPUT BITMAP FLAGS
!     LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
!     GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
!
!   OUTPUT ARGUMENT LIST:
!     NO       - INTEGER NUMBER OF OUTPUT POINTS
!     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
!     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
!     IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
!     LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
!     GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
!     IRET     - INTEGER RETURN CODE
!                0    SUCCESSFUL INTERPOLATION
!                2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
!                3    UNRECOGNIZED OUTPUT GRID
!                31   INVALID UNDEFINED OUTPUT GRID
!                32   INVALID BUDGET METHOD PARAMETERS
!
! SUBPROGRAMS CALLED:
!   GDSWZD       GRID DESCRIPTION SECTION WIZARD
!   IJKGDS0      SET UP PARAMETERS FOR IJKGDS1
!   IJKGDS1      RETURN FIELD POSITION FOR A GIVEN GRID POINT
!   POLFIXS      MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 USE GDSWZD_MOD
!
 IMPLICIT NONE
!
 INTEGER,        INTENT(IN   ) :: IGDTNUMI, IGDTLENI
 INTEGER,        INTENT(IN   ) :: IGDTMPLI(IGDTLENI)
 INTEGER,        INTENT(IN   ) :: IGDTNUMO, IGDTLENO
 INTEGER,        INTENT(IN   ) :: IGDTMPLO(IGDTLENO)
 INTEGER,        INTENT(IN   ) :: IBI(KM), IPOPT(20), KM, MI, MO
 INTEGER,        INTENT(  OUT) :: IBO(KM), IRET, NO
!
 LOGICAL*1,      INTENT(IN   ) :: LI(MI,KM)
 LOGICAL*1,      INTENT(  OUT) :: LO(MO,KM)
!
 REAL,           INTENT(IN   ) :: GI(MI,KM)
 REAL,           INTENT(  OUT) :: GO(MO,KM), RLAT(MO), RLON(MO)
!
 REAL,       PARAMETER         :: FILL=-9999.
!
 INTEGER                       :: IB, I1, IJKGDS1, IJKGDSA(20)
 INTEGER                       :: JB, J1, K, LB, LSW, MP, N
 INTEGER                       :: N11(MO), NB, NB1, NB2, NB3, NB4, NV
!
 REAL                          :: PMP,RLOB(MO),RLAB(MO)
 REAL                          :: WB, WO(MO,KM), XI, YI
 REAL                          :: XPTB(MO),YPTB(MO),XPTS(MO),YPTS(MO)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
 IRET=0
 IF(IGDTNUMO.GE.0) THEN
   CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
   IF(NO.EQ.0) IRET=3
 ELSE
   IRET=31
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SET PARAMETERS
 NB1=IPOPT(1)
 IF(NB1.EQ.-1) NB1=2
 IF(IRET.EQ.0.AND.NB1.LT.0) IRET=32
 LSW=1
 IF(IPOPT(1).EQ.-1.OR.IPOPT(2).EQ.-1) LSW=0
 IF(IRET.EQ.0.AND.LSW.EQ.1.AND.NB1.GT.15) IRET=32
 MP=IPOPT(3+IPOPT(1))
 IF(MP.EQ.-1.OR.MP.EQ.0) MP=50
 IF(MP.LT.0.OR.MP.GT.100) IRET=32
 PMP=MP*0.01
 IF(IRET.EQ.0) THEN
   NB2=2*NB1+1
   NB3=NB2*NB2
   NB4=NB3
   IF(LSW.EQ.1) THEN
     NB4=IPOPT(2)
     DO IB=1,NB1
       NB4=NB4+8*IB*IPOPT(2+IB)
     ENDDO
   ENDIF
 ELSE
   NB2=0
   NB3=0
   NB4=0
 ENDIF
 DO K=1,KM
   DO N=1,NO
     GO(N,K)=0.
     WO(N,K)=0.
   ENDDO
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
 CALL IJKGDS0(IGDTNUMI,IGDTMPLI,IGDTLENI,IJKGDSA)
 DO NB=1,NB3
!  LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
   JB=(NB-1)/NB2-NB1
   IB=NB-(JB+NB1)*NB2-NB1-1
   LB=MAX(ABS(IB),ABS(JB))
   WB=1
   IF(LSW.EQ.1) WB=IPOPT(2+LB)
   IF(WB.NE.0) THEN
     DO N=1,NO
       XPTB(N)=XPTS(N)+IB/REAL(NB2)
       YPTB(N)=YPTS(N)+JB/REAL(NB2)
     ENDDO
     CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     CALL GDSWZD(IGDTNUMI,IGDTMPLI,IGDTLENI,-1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     IF(IRET.EQ.0.AND.NV.EQ.0.AND.LB.EQ.0) IRET=2
     DO N=1,NO
       XI=XPTB(N)
       YI=YPTB(N)
       IF(XI.NE.FILL.AND.YI.NE.FILL) THEN
         I1=NINT(XI)
         J1=NINT(YI)
         N11(N)=IJKGDS1(I1,J1,IJKGDSA)
       ELSE
         N11(N)=0
       ENDIF
     ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  INTERPOLATE WITH OR WITHOUT BITMAPS
     DO K=1,KM
       DO N=1,NO
         IF(N11(N).GT.0) THEN
           IF(IBI(K).EQ.0.OR.LI(N11(N),K)) THEN
             GO(N,K)=GO(N,K)+WB*GI(N11(N),K)
             WO(N,K)=WO(N,K)+WB
           ENDIF
         ENDIF
       ENDDO
     ENDDO
   ENDIF
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE OUTPUT BITMAPS AND FIELDS
 DO K=1,KM
   IBO(K)=IBI(K)
   DO N=1,NO
     LO(N,K)=WO(N,K).GE.PMP*NB4
     IF(LO(N,K)) THEN
       GO(N,K)=GO(N,K)/WO(N,K)
     ELSE
       IBO(K)=1
       GO(N,K)=0.
     ENDIF
   ENDDO
 ENDDO
 IF(IGDTNUMO.EQ.0) CALL POLFIXS(NO,MO,KM,RLAT,IBO,LO,GO)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE POLATES6

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM