text:
C $Header: /u/gcmpack/MITgcm/verification/tutorial_plume_on_slope/code/CPP_OPTIONS.h,v 1.1 2006/04/20 16:45:12 jmc Exp $
C $Name: checkpoint65t $

C CPP flags controlling particular source code features

C o Shortwave heating as extra term in external_forcing.F
C Note: this should be a run-time option

C o Include/exclude phi_hyd calculation code

C o Include/exclude call to S/R CONVECT

C o Include/exclude call to S/R CALC_DIFFUSIVITY

C o Include/exclude Implicit vertical advection code

C o Include/exclude nonHydrostatic code

C o Include pressure loading code

C o Use "Exact Convervation" of fluid in Free-Surface formulation
C   so that d/dt(eta) is exactly equal to - Div.Transport

C o Allow the use of Non-Linear Free-Surface formulation
C   this implies that surface thickness (hFactors) vary with time

C o Use "OLD" UV discretisation near boundaries (*not* recommended)
C   Note - only works with  #undef NO_SLIP_LATERAL  in calc_mom_rhs.F
C          because the old code did not have no-slip BCs

C o Execution environment support options
CBOP
C     !ROUTINE: CPP_EEOPTIONS.h
C     !INTERFACE:
C     include "CPP_EEOPTIONS.h"
C
C     !DESCRIPTION:
C     *==========================================================*
C     | CPP\_EEOPTIONS.h                                         |
C     *==========================================================*
C     | C preprocessor "execution environment" supporting        |
C     | flags. Use this file to set flags controlling the        |
C     | execution environment in which a model runs - as opposed |
C     | to the dynamical problem the model solves.               |
C     | Note: Many options are implemented with both compile time|
C     |       and run-time switches. This allows options to be   |
C     |       removed altogether, made optional at run-time or   |
C     |       to be permanently enabled. This convention helps   |
C     |       with the data-dependence analysis performed by the |
C     |       adjoint model compiler. This data dependency       |
C     |       analysis can be upset by runtime switches that it  |
C     |       is unable to recoginise as being fixed for the     |
C     |       duration of an integration.                        |
C     |       A reasonable way to use these flags is to          |
C     |       set all options as selectable at runtime but then  |
C     |       once an experimental configuration has been        |
C     |       identified, rebuild the code with the appropriate  |
C     |       options set at compile time.                       |
C     *==========================================================*
CEOP

C     In general the following convention applies:
C     ALLOW  - indicates an feature will be included but it may
C     CAN      have a run-time flag to allow it to be switched
C              on and off.
C              If ALLOW or CAN directives are "undef'd" this generally
C              means that the feature will not be available i.e. it
C              will not be included in the compiled code and so no
C              run-time option to use the feature will be available.
C
C     ALWAYS - indicates the choice will be fixed at compile time
C              so no run-time option will be present

C=== Macro related options ===
C--   Control storage of floating point operands
C     On many systems it improves performance only to use
C     8-byte precision for time stepped variables.
C     Constant in time terms ( geometric factors etc.. )
C     can use 4-byte precision, reducing memory utilisation and
C     boosting performance because of a smaller working set size.
C     However, on vector CRAY systems this degrades performance.
C     Enable to switch REAL4_IS_SLOW from genmake2 (with LET_RS_BE_REAL4):

C--   Control use of "double" precision constants.
C     Use D0 where it means REAL*8 but not where it means REAL*16

C--   Enable some old macro conventions for backward compatibility

C=== IO related options ===
C--   Flag used to indicate whether Fortran formatted write
C     and read are threadsafe. On SGI the routines can be thread
C     safe, on Sun it is not possible - if you are unsure then
C     undef this option.

C--   Flag used to indicate whether Binary write to Local file (i.e.,
C     a different file for each tile) and read are thread-safe.

C--   Flag to turn off the writing of error message to ioUnit zero

C--   Alternative formulation of BYTESWAP, faster than
C     compiler flag -byteswapio on the Altix.

C--   Flag to turn on old default of opening scratch files with the
C     STATUS='SCRATCH' option. This method, while perfectly FORTRAN-standard,
C     caused filename conflicts on some multi-node/multi-processor platforms
C     in the past and has been replace by something (hopefully) more robust.

C--   Flag defined for eeboot_minimal.F, eeset_parms.F and open_copy_data_file.F
C     to write STDOUT, STDERR and scratch files from process 0 only.
C WARNING: to use only when absolutely confident that the setup is working
C     since any message (error/warning/print) from any proc <> 0 will be lost.

C=== MPI, EXCH and GLOBAL_SUM related options ===
C--   Flag turns off MPI_SEND ready_to_receive polling in the
C     gather_* subroutines to speed up integrations.

C--   Control MPI based parallel processing
CXXX We no longer select the use of MPI via this file (CPP_EEOPTIONS.h)
CXXX To use MPI, use an appropriate genmake2 options file or use
CXXX genmake2 -mpi .
CXXX #undef  1

C--   Control use of communication that might overlap computation.
C     Under MPI selects/deselects "non-blocking" sends and receives.
C--   Control use of communication that is atomic to computation.
C     Under MPI selects/deselects "blocking" sends and receives.

C--   Control XY periodicity in processor to grid mappings
C     Note: Model code does not need to know whether a domain is
C           periodic because it has overlap regions for every box.
C           Model assume that these values have been
C           filled in some way.

C--   disconnect tiles (no exchange between tiles, just fill-in edges
C     assuming locally periodic subdomain)

C--   Always cumulate tile local-sum in the same order by applying MPI allreduce
C     to array of tiles ; can get slower with large number of tiles (big set-up)

C--   Alternative way of doing global sum without MPI allreduce call
C     but instead, explicit MPI send & recv calls. Expected to be slower.

C--   Alternative way of doing global sum on a single CPU
C     to eliminate tiling-dependent roundoff errors. Note: This is slow.

C=== Other options (to add/remove pieces of code) ===
C--   Flag to turn on checking for errors from all threads and procs
C     (calling S/R STOP_IF_ERROR) before stopping.

C--   Control use of communication with other component:
C     allow to import and export from/to Coupler interface.

C--   Activate some pieces of code for coupling to GEOS AGCM

CBOP
C     !ROUTINE: CPP_EEMACROS.h
C     !INTERFACE:
C     include "CPP_EEMACROS.h"
C     !DESCRIPTION:
C     *==========================================================*
C     | CPP_EEMACROS.h
C     *==========================================================*
C     | C preprocessor "execution environment" supporting
C     | macros. Use this file to define macros for  simplifying
C     | execution environment in which a model runs - as opposed
C     | to the dynamical problem the model solves.
C     *==========================================================*
CEOP

C     In general the following convention applies:
C     ALLOW  - indicates an feature will be included but it may
C     CAN      have a run-time flag to allow it to be switched
C              on and off.
C              If ALLOW or CAN directives are "undef'd" this generally
C              means that the feature will not be available i.e. it
C              will not be included in the compiled code and so no
C              run-time option to use the feature will be available.
C
C     ALWAYS - indicates the choice will be fixed at compile time
C              so no run-time option will be present

C     Flag used to indicate which flavour of multi-threading
C     compiler directives to use. Only set one of these.
C     USE_SOLARIS_THREADING  - Takes directives for SUN Workshop
C                              compiler.
C     USE_KAP_THREADING      - Takes directives for Kuck and
C                              Associates multi-threading compiler
C                              ( used on Digital platforms ).
C     USE_IRIX_THREADING     - Takes directives for SGI MIPS
C                              Pro Fortran compiler.
C     USE_EXEMPLAR_THREADING - Takes directives for HP SPP series
C                              compiler.
C     USE_C90_THREADING      - Takes directives for CRAY/SGI C90
C                              system F90 compiler.

C--   Define the mapping for the _BARRIER macro
C     On some systems low-level hardware support can be accessed through
C     compiler directives here.

C--   Define the mapping for the BEGIN_CRIT() and  END_CRIT() macros.
C     On some systems we simply execute this section only using the
C     master thread i.e. its not really a critical section. We can
C     do this because we do not use critical sections in any critical
C     sections of our code!

C--   Define the mapping for the BEGIN_MASTER_SECTION() and
C     END_MASTER_SECTION() macros. These are generally implemented by
C     simply choosing a particular thread to be "the master" and have
C     it alone execute the BEGIN_MASTER..., END_MASTER.. sections.

CcnhDebugStarts
C      Alternate form to the above macros that increments (decrements) a counter each
C      time a MASTER section is entered (exited). This counter can then be checked in barrier
C      to try and detect calls to BARRIER within single threaded sections.
C      Using these macros requires two changes to Makefile - these changes are written
C      below.
C      1 - add a filter to the CPP command to kill off commented _MASTER lines
C      2 - add a filter to the CPP output the converts the string N EWLINE to an actual newline.
C      The N EWLINE needs to be changes to have no space when this macro and Makefile changes
C      are used. Its in here with a space to stop it getting parsed by the CPP stage in these
C      comments.
C      #define IF ( a .EQ. 1 ) THEN  IF ( a .EQ. 1 ) THEN  N EWLINE      CALL BARRIER_MS(a)
C      #define ENDIF    CALL BARRIER_MU(a) N EWLINE        ENDIF
C      'CPP = cat $< | $(TOOLSDIR)/set64bitConst.sh |  grep -v '^[cC].*_MASTER' | cpp  -traditional -P'
C      .F.f:
C      $(CPP) $(DEFINES) $(INCLUDES) |  sed 's/N EWLINE/\n/' > $@
CcnhDebugEnds

C--   Control storage of floating point operands
C     On many systems it improves performance only to use
C     8-byte precision for time stepped variables.
C     Constant in time terms ( geometric factors etc.. )
C     can use 4-byte precision, reducing memory utilisation and
C     boosting performance because of a smaller working
C     set size. However, on vector CRAY systems this degrades
C     performance.
C- Note: global_sum/max macros were used to switch to  JAM routines (obsolete);
C  in addition, since only the R4 & R8 S/R are coded, GLOBAL RS & RL macros
C  enable to call the corresponding R4 or R8 S/R.

C- Note: a) exch macros were used to switch to  JAM routines (obsolete)
C        b) exch R4 & R8 macros are not practically used ; if needed,
C           will directly call the corrresponding S/R.

C--   Control use of JAM routines for Artic network (no longer supported)
C     These invoke optimized versions of "exchange" and "sum" that
C     utilize the programmable aspect of Artic cards.
CXXX No longer supported ; started to remove JAM routines.
CXXX #ifdef LETS_MAKE_JAM
CXXX #define CALL GLOBAL_SUM_R8 ( a, b) CALL GLOBAL_SUM_R8_JAM ( a, b)
CXXX #define CALL GLOBAL_SUM_R8 ( a, b ) CALL GLOBAL_SUM_R8_JAM ( a, b )
CXXX #define CALL EXCH_XY_RS ( a, b ) CALL EXCH_XY_R8_JAM ( a, b )
CXXX #define CALL EXCH_XY_RL ( a, b ) CALL EXCH_XY_R8_JAM ( a, b )
CXXX #define CALL EXCH_XYZ_RS ( a, b ) CALL EXCH_XYZ_R8_JAM ( a, b )
CXXX #define CALL EXCH_XYZ_RL ( a, b ) CALL EXCH_XYZ_R8_JAM ( a, b )
CXXX #endif

C--   Control use of "double" precision constants.
C     Use d0 where it means REAL*8 but not where it means REAL*16

C--   Substitue for 1.D variables
C     Sun compilers do not use 8-byte precision for literals
C     unless .Dnn is specified. CRAY vector machines use 16-byte
C     precision when they see .Dnn which runs very slowly!

C--   Set the format for writing processor IDs, e.g. in S/R eeset_parms
C     and S/R open_copy_data_file. The default of I9.9 should work for
C     a long time (until we will use 10e10 processors and more)

C o Include/exclude code specific to the ECCO/SEALION version.
C   AUTODIFF or EXF package.
C   Currently controled by a single header file
C   For this to work, PACKAGES_CONFIG.h needs to be included!
cph#if (defined (ALLOW_AUTODIFF) || cph     defined (ALLOW_ECCO) || cph     defined (ALLOW_EXF))
cph# include "ECCO_CPPOPTIONS.h"
cph#endif

C--  File gsw_teos10.F: routines that compute quantities related to seawater.
C--   Contents
C--   TEOS-10 routines (Gibbs seawater, GSW)
C--   o GSW_PT_FROM_CT: function to compute potential temperature
C--              from conservative temperature and absolute salinity
C--   o GSW_CT_FROM_PT: function to compute conservative temperature with
C--              from potential temperature and absolute salinity
C--   o GSW_GIBBS_PT0_PT0: function to compute specific Gibbs free energy
C--              from potential temperature and absolute salinity

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: GSW_PT_FROM_CT
C     !INTERFACE:
      Real*8 FUNCTION GSW_PT_FROM_CT(SA,CT)
C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R  GSW_PT_FROM_CT
C     | o compute potential temperature at reference level 0 dbar
C     |   from conservative temperature (CT) and absolute
C     |    salinity (SA)
C     | o this is a more or less shameless copy of the teos-10 code
C     |   available at http://www.teos-10.org
C     *=============================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     SA :: Absolute Salinity                               (g/kg)
C     CT :: Conservative Temperature                        (deg C)
      Real*8 SA,CT

C     !FUNCTIONS:
C     == Functions ==
CML      Real*8 gsw_gibbs
      Real*8 gsw_gibbs_pt0_pt0
      Real*8 gsw_ct_from_pt
CML      EXTERNAL gsw_gibbs, gsw_gibbs_pt0_pt0, gsw_ct_from_pt
      EXTERNAL gsw_gibbs_pt0_pt0, gsw_ct_from_pt

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER n0, n2
      Real*8 s1, p0, cp0
      Real*8 a0, a1, a2, a3, a4, a5, b0, b1, b2, b3
      Real*8 a5ct, b3ct, ct_factor, pt_num, pt_den, ct_diff
      Real*8 pt, pt_old, ptm, dct_dpt
CEOP

      cp0 = 3991.86795711963D0

      n0 = 0
      n2 = 2

      s1 = SA * 35.D0 / 35.16504D0
      p0 = 0.D0

      a0 = -1.446013646344788D-2
      a1 = -3.305308995852924D-3
      a2 =  1.062415929128982D-4
      a3 =  9.477566673794488D-1
      a4 =  2.166591947736613D-3
      a5 =  3.828842955039902D-3

      b0 =  1.000000000000000D+0
      b1 =  6.506097115635800D-4
      b2 =  3.830289486850898D-3
      b3 =  1.247811760368034D-6

      a5ct = a5*CT
      b3ct = b3*CT

      ct_factor = (a3 + a4*s1 + a5ct)
      pt_num    = a0 + s1*(a1 + a2*s1) + CT*ct_factor
      pt_den    = b0 + b1*s1 + CT*(b2 + b3ct)
      pt        = (pt_num)/(pt_den)

      dct_dpt   = (pt_den)/(ct_factor + a5ct - (b2 + b3ct + b3ct)*pt)

C     start the 1.5 iterations through the modified Newton-Rapshon
C     iterative method.

      ct_diff = gsw_ct_from_pt(sa,pt) - CT
      pt_old  = pt
      pt      = pt_old - (ct_diff)/dct_dpt
      ptm     = 0.5D0*(pt + pt_old)

      dct_dpt = -(ptm + 273.15D0)*gsw_gibbs_pt0_pt0(sa,ptm)/cp0

      pt             = pt_old - (ct_diff)/dct_dpt
      ct_diff        = gsw_ct_from_pt(sa,pt) - CT
      pt_old         = pt
      GSW_PT_FROM_CT = pt_old - (ct_diff)/dct_dpt

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: GSW_CT_FROM_PT
C     !INTERFACE:
      Real*8 FUNCTION GSW_CT_FROM_PT(SA,PT)
C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R  GSW_CT_FROM_PT
C     | o compute conservative temperature from potential
C     |   temperature (PT)  at reference level 0 dbar and absolute
C     |   salinity (SA)
C     | o this is a more or less shameless copy of the teos-10 code
C     |   available at http://www.teos-10.org
C     *=============================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     SA :: Absolute Salinity                               (g/kg)
C     PT :: Potential Temperature                          (deg C)
      Real*8 sa, pt

C     !FUNCTIONS:
C     == Functions ==

C     !LOCAL VARIABLES:
C     == Local variables ==
      Real*8 pot_enthalpy, sfac
      Real*8 x2, x, y, cp0
CEOP

      sfac = 0.0248826675584615D0

      x2   = sfac*sa
      x    = 0.D0
      if (x2.gt.0.D0) x = sqrt(x2)
C     normalize for F03 and F08
      y    = pt*0.025D0

      pot_enthalpy =  61.01362420681071D0 +
     &     y*(168776.46138048015D0 +
     &     y*(-2735.2785605119625D0 + y*(2574.2164453821433D0 +
     &     y*(-1536.6644434977543D0 + y*(545.7340497931629D0 +
     &     (-50.91091728474331D0 - 18.30489878927802D0*y)*y))))) +
     &     x2*(268.5520265845071D0 + y*(-12019.028203559312D0 +
     &     y *(3734.858026725145D0 + y*(-2046.7671145057618D0 +
     &     y*(465.28655623826234D0 + (-0.6370820302376359D0 -
     &     10.650848542359153D0*y)*y)))) +
     &     x*(937.2099110620707D0 + y*(588.1802812170108D0 +
     &     y*(248.39476522971285D0 + (-3.871557904936333D0 -
     &     2.6268019854268356D0*y)*y)) +
     &     x*(-1687.914374187449D0 + x*(246.9598888781377D0 +
     &     x*(123.59576582457964D0 - 48.5891069025409D0*x)) +
     &     y*( 936.3206544460336D0 +
     &     y*(-942.7827304544439D0 + y*(369.4389437509002D0 +
     &     (-33.83664947895248D0 - 9.987880382780322D0*y)*y))))))

      cp0 = 3991.86795711963D0

      gsw_ct_from_pt = pot_enthalpy/cp0

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: GSW_GIBBS_PT0_PT0
C     !INTERFACE:
      Real*8 FUNCTION GSW_GIBBS_PT0_PT0(SA,PT0)
C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R GSW_GIBBS_PT0_PT0
C     | o helper routine that computes the specific Gibbs free
C     |   energy from potential temperature (PT) and absolute
C     |   salinity (SA) for pressure 0 dbar
C     | o this is a more or less shameless copy of the teos-10 code
C     |   available at http://www.teos-10.org
C     *=============================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     SA :: Absolute Salinity                               (g/kg)
C     PT :: Potential Temperature at p = 0 dbar             (deg C)
      Real*8 sa, pt0

C     !FUNCTIONS:
C     == Functions ==

C     !LOCAL VARIABLES:
C     == Local variables ==
      Real*8 sfac, x2, x, y, g03, g08
CEOP

      sfac = 0.0248826675584615

      x2   = sfac*sa
      x    = 0.D0
      if (x2.gt.0.D0) x = sqrt(x2)
      y    = pt0*0.025D0

      g03 = -24715.571866078 +
     &     y*(4420.4472249096725 +
     &     y*(-1778.231237203896 +
     &     y*(1160.5182516851419 +
     &     y*(-569.531539542516  + y*128.13429152494615))))

      g08 = x2*(1760.062705994408  + x*(-86.1329351956084 +
     &     x*( -137.1145018408982  + y*(296.20061691375236 +
     &     y* (-205.67709290374563 + 49.9394019139016*y))) +
     &     y*(  -60.136422517125   + y*10.50720794170734)) +
     &     y*(-1351.605895580406   + y*(1097.1125373015109 +
     &     y*( -433.20648175062206 + 63.905091254154904*y))))

      gsw_gibbs_pt0_pt0 = (g03 + g08)*0.000625

      RETURN
      END

max_stars_repo_name: ElizabethYankovsky/MITgcm

max_issues_repo_name: ElizabethYankovsky/MITgcm

max_forks_repo_name: ElizabethYankovsky/MITgcm


text:
module simplication
	implicit none
	contains
	  ! This subroutine will calculate the x and y values for n number of points
		subroutine calc_xy(x,y,delta_x,n)
		implicit none
		integer, intent(IN)::n
		integer :: i
		real, intent(out)::x(n),y(n) !arrays of points for setting up the polynomial
		real, intent(IN) ::delta_x

		do i=2,n-1
			x(i)=x(i-1)+delta_x !calculatin x coordintes for setting up the polynomial
		end do
		do i=1,n
			y(i)=1.0d0/(1.0d0+25.0d0*x(i)*x(i)) !calculating y coordinates for setting up the polynomical from exact solution
		   !y(i)=cos(x(i)*3.141592654)
		end do

		do i=1,n
			write(25,*) x(i),y(i)
		end do
		end subroutine calc_xy

		!to calculate x and y values for m number of points from the exact solution
		subroutine calc_xf_yf(x_f,y_f,delta_x_f,m)
		implicit none
		integer, intent(IN)::m
		integer :: i
		real, intent(out)::x_f(m),y_f(m) !arrays of points for plotting the polynomial at m points
		real, intent(IN) ::delta_x_f

		do i=2,m-1
			x_f(i)=x_f(i-1)+delta_x_f !calculatin x coordintes for setting up the polynomial
		end do
		!calculating f(x) from the exact solution for the corrosponding values of x(n)
		do i=1,m
			y_f(i)=1.0d0/(1.0d0+25.0d0*x_f(i)*x_f(i)) !calculating y coordinates for setting up the polynomical from exact solution
        !y(i)=cos(x(i)*3.141592654)

		end do

		do i=1,m
			write(16,*) x_f(i),y_f(i)
		end do

		end subroutine calc_xf_yf

		!to calculate the coeficients for newton interpolation
		subroutine calc_alpha(x,y,alpha,n)
		implicit none
		integer,intent(IN):: n
		integer ::i ,j
		real,intent(INOUT) ::x(n),y(n),alpha(n)

		alpha(1)=y(1)
		Do i= 1,n-1
			Do j=1, n-i

				y(j)=(y(j+1)-y(j))/(x(j+i)-x(j))

			end do
			alpha(i+1)=y(1)
		end do
	  end subroutine calc_alpha

		!to calculate y=f(x) from nth value polynomial for x_f(m) number of points
		subroutine calc_poly(x,x_f,alpha,p,n,m)
		implicit none
		integer,intent(IN):: n,m
		integer ::i,j
		real,intent(IN) ::x(n),x_f(m),alpha(n)
		real ::p
		p=alpha(n)

		do j=1,m
		p=alpha(n)

			do i=n,2,-1

				p=(p*(x_f(j)-(x(i-1))))+alpha(i-1)

			end do
			write(43,*) x_f(j),p
		end do
		end subroutine calc_poly

end module simplification

max_stars_repo_name: akkiturb/Ilkay-s-Fortran-Class


program tarefac
    implicit real(8) (a-h,o-z)

    dt = 0.001d0
    tmax = 240d0

    g = 9.8d0
    a_l = 9.8d0
    gama = 0.5d0
    Omega = 2d0/3d0

    write(*,*) 'Entre com o valor de F0:'
    read(*,*) F0

    open(10, file='saida-d-10407962.dat')

    theta_i = 0.2d0
    omega_i = 0d0
    i = 0
    n = 0
    a_lyapunov = 0d0

    do while ( i*dt.le.tmax )

        dw_i = dw(g, a_l, theta_i, gama, omega_i, F0, Omega, dt*i)
        call oscilacao(theta_i, omega_i, dw_i, dt)

        write(10, *) theta_i, omega_i

        i = i+1
    end do

    close(10)

contains
    function dw(g, a_l, theta, gama, dTheta, F0, Omega, t)
        implicit real(8) (a-h,o-z)
        dw = -(g/a_l)*sin(theta) -gama*dTheta +F0*sin(Omega*t)
        return
    end function

    subroutine oscilacao(theta_i, omega_i, dw_i, dt)
        implicit real(8) (a-h,o-z)
        parameter (pi = 4d0*atan(1d0))
        omega_i = omega_i +dw_i*dt
        theta_i = theta_i +omega_i*dt

        do while (abs(theta_i).ge.pi)
            sinal = theta_i/abs(theta_i)
            theta_i = theta_i -sinal*2*pi
        end do
        return
    end subroutine

end program tarefac


text:
! Copyright (c) 1998, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! test that PRIVATE/PUBLIC can be used on names that aren't defined until later

module m1
 real max,min
 public:: max
 private:: min
 public:: sin, setmin
contains
 subroutine setmin(i)
  real i
  min = i
!print *,'min=',min
 end subroutine
 real function sin(x)
  real x
  sin = min*max*x
!print *,'sin=',sin
 end function
end module

program p
 use m1
 real t, u, v
 real expect(3), result(3)
 data expect/5.0,90.0,5.0/
 max = 9.0
 call setmin(2.0)
 t = 5.0
 u = sin(t)
 v = min(t,u)
 result(1) = t
 result(2) = u
 result(3) = v
! print *,t,u,v
 call check(result,expect,3)
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


SUBROUTINE INTH (DP,DQ,DEP,DEQ,DB,DVC,DSAL,DK,DM,DV,DR)
  !
  ! INTEGRATION BY THE METHOD OF ADAMS, A 5 POINTS OF COMPONENT DP GREAT AND SMALL COMPONENT DQ POINT DR, DV IS THE POTENTIAL IN POINT
  !
  ! INTEGRATION PAR LA METHODE DE ADAMS A 5 POINTS DE LA GRANDE COMPOSANTE
  ! DP ET DE LA PETITE COMPOSANTE DQ AU POINT DR,DV ETANT LE POTENTIEL EN
  ! CE POINT
  ! **********************************************************************
  !USE param
  !IMPLICIT REAL*8 (A-H,O-Z)
  IMPLICIT NONE
  REAL*8, intent(inout) :: DP, DQ
  REAL*8, intent(inout) :: DEP(5), DEQ(5)
  REAL*8, intent(in)    :: DB, DVC, DSAL, DK, DM, DV, DR
  !     IMPLICIT REAL (D)
  !COMMON/PS1/DEP(5),DEQ(5),DB,DVC,DSAL,DK,DM
  !save  /PS1/
  !
  ! DEP,DEQ DERIVEES DE DP ET DQ   DB=ENERGIE/DVC    DVC VITESSE DE LA
  ! LUMIERE EN U.A.   DSAL=2.*DVC   DK NOMBRE QUANTIQUE KAPPA
  ! DM=PAS EXPONENTIEL/720.
  ! DKOEF1=405./502., DKOEF2=27./502.
  ! *********************************************************************
  REAL*8 :: dkoef1, dkoef2, dpr, dqr, dsum
  INTEGER:: i
  !
  DATA DKOEF1/.9462151394422310D0/,DKOEF2/.5378486055776890D-1/
  DPR=DP+DM*((251.*DEP(1)+2616.*DEP(3)+1901.*DEP(5))-(1274.*DEP(2)+2774.*DEP(4)))
  DQR=DQ+DM*((251.*DEQ(1)+2616.*DEQ(3)+1901.*DEQ(5))-(1274.*DEQ(2)+2774.*DEQ(4)))
  DO I=2,5
     DEP(I-1)=DEP(I)
     DEQ(I-1)=DEQ(I)
  ENDDO
  DSUM=(DB-DV/DVC)*DR
  DEP(5)=-DK*DPR+(DSAL*DR+DSUM)*DQR
  DEQ(5)=DK*DQR-DSUM*DPR
  DP=DP+DM*((106.*DEP(2)+646.*DEP(4)+251.*DEP(5))-(19.*DEP(1)+264.*DEP(3)))
  DQ=DQ+DM*((106.*DEQ(2)+646.*DEQ(4)+251.*DEQ(5))-(19.*DEQ(1)+264.*DEQ(3)))
  DP=DKOEF1*DP+DKOEF2*DPR
  DQ=DKOEF1*DQ+DKOEF2*DQR
  DEP(5)=-DK*DP+(DSAL*DR+DSUM)*DQ
  DEQ(5)=DK*DQ-DSUM*DP
  RETURN
END SUBROUTINE INTH


text: #redirect wiki:Sacramento:Sacramento Bike Kitchen

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!>
!! Copyright 2016 ARC Centre of Excellence for Climate Systems Science
!!
!! \author  Scott Wales <scott.wales@unimelb.edu.au>
!!
!! Licensed under the Apache License, Version 2.0 (the "License");
!! you may not use this file except in compliance with the License.
!! You may obtain a copy of the License at
!!
!!     http://www.apache.org/licenses/LICENSE-2.0
!!
!! Unless required by applicable law or agreed to in writing, software
!! distributed under the License is distributed on an "AS IS" BASIS,
!! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!! See the License for the specific language governing permissions and
!! limitations under the License.

module farray_mod
    private

    public open_dataset

    ! Represents a netCDF/OPENDAP file
    type, public :: dataset_t
        integer :: ncid
    contains
        procedure :: variable
        procedure :: close => close_dataset
    end type

    ! Represents a variable within the file
    type, public :: variable_t
        integer :: ncid
        integer :: varid
        integer, allocatable :: dimid(:)
    contains
        procedure :: dimension_index => var_dimension_index
        generic :: dimension => dimension_index

        procedure :: shape => var_shape

        procedure :: read_int32
        procedure :: read_int64
        procedure :: read_real32
        procedure :: read_real64
        generic :: read => read_int32, read_int64, read_real32, read_real64

!        procedure :: read_int32_alloc
!        procedure :: read_int64_alloc
!        procedure :: read_real32_alloc
!        procedure :: read_real64_alloc
!        generic :: read_alloc => read_int32_alloc, read_int64_alloc, read_real32_alloc, read_real64_alloc
    end type

    ! Represents a dimension
    type, public :: dimension_t
        integer :: ncid
        integer :: dimid
    contains
        procedure :: size => dimension_size
    end type

    interface open_dataset
        procedure :: open_dataset_f08
        procedure :: open_dataset_f90
    end interface

contains

    function open_dataset_f08(path, comm) result(r)
        use mpi_f08
        use netcdf4_f03
        character(len=*), intent(in) :: path
        type(mpi_comm), intent(in) :: comm
        type(dataset_t) :: r

        r = open_dataset_f90(path, comm%mpi_val)
    end function

    function open_dataset_f90(path, comm) result(r)
        use mpi_f08
        use netcdf4_f03
        character(len=*), intent(in) :: path
        integer, intent(in) :: comm
        type(dataset_t) :: r
        integer :: mode
        integer :: err

        mode = NF_NOWRITE
        err = nf_open_par(path, mode, comm=comm, info=MPI_INFO_NULL%mpi_val, ncid=r%ncid)

        if (err == NF_ENOPAR) then
            ! No parallel netcdf support
            err = nf_open(path, mode, ncid=r%ncid)
        end if
        call check(err)
    end function

    subroutine close_dataset(this)
        use netcdf4_f03
        class(dataset_t), intent(inout) :: this
        integer :: err

        err = nf_close(this%ncid)
        call check(err)
    end subroutine

    function variable(this, name) result(r)
        use netcdf4_f03
        class(dataset_t), intent(in) :: this
        character(len=*), intent(in) :: name
        type(variable_t) :: r
        integer :: ndims
        integer :: err

        r%ncid = this%ncid
        err = nf_inq_varid(r%ncid, name, r%varid)
        call check(err)

        err = nf_inq_varndims(r%ncid, r%varid, ndims)
        call check(err)

        allocate(r%dimid(ndims))

        err = nf_inq_vardimid(r%ncid, r%varid, r%dimid)
        call check(err)
    end function

    function var_shape(this) result(r)
        class(variable_t), intent(in) :: this
        integer, allocatable :: r(:)
        type(dimension_t) :: dd
        integer :: ii

        allocate(r(size(this%dimid)))
        do ii=1,size(this%dimid)
            dd = this%dimension(ii)
            r(ii) = dd%size()
        end do
    end function

    function var_dimension_index(this, idx) result(r)
        use netcdf4_f03
        class(variable_t), intent(in) :: this
        integer, intent(in) :: idx
        type(dimension_t) :: r

        r%ncid = this%ncid
        r%dimid = this%dimid(idx)
    end function

    function dimension_size(this) result(r)
        use netcdf4_f03
        class(dimension_t) :: this
        integer :: r
        integer :: err

        err = nf_inq_dimlen(this%ncid, this%dimid, r)
        call check(err)
    end function

    subroutine read_int32(this, array)
        use netcdf4_f03
        class(variable_t), intent(in) :: this
        integer(kind=4), dimension(..) :: array
        integer :: err

        err = nf_get_vara_int(this%ncid, this%varid, lbound(array), shape(array), array)
        call check(err)
    end subroutine

    subroutine read_int64(this, array)
        use netcdf4_f03
        class(variable_t), intent(in) :: this
        integer(kind=8), dimension(..) :: array
        integer :: err

        err = nf_get_vara_int64(this%ncid, this%varid, lbound(array), shape(array), array)
        call check(err)
    end subroutine

    subroutine read_real32(this, array)
        use netcdf4_f03
        class(variable_t), intent(in) :: this
        real(kind=4), dimension(..) :: array
        integer :: err

        err = nf_get_vara_real(this%ncid, this%varid, lbound(array), shape(array), array)
        call check(err)
    end subroutine

    subroutine read_real64(this, array)
        use netcdf4_f03
        class(variable_t), intent(in) :: this
        real(kind=8), dimension(..) :: array
        integer :: err

        err = nf_get_vara_double(this%ncid, this%varid, lbound(array), shape(array), array)
        call check(err)
    end subroutine

!    subroutine read_int32_alloc(this, array)
!        use netcdf4_f03
!        class(variable_t), intent(in) :: this
!        integer(kind=4), allocatable, dimension(..) :: array
!        integer :: err
!
!        allocate(array(this%shape()))
!        err = nf_get_vara_int(this%ncid, this%varid, lbound(array), shape(array), array)
!        call check(err)
!    end subroutine
!
!    subroutine read_int64_alloc(this, array)
!        use netcdf4_f03
!        class(variable_t), intent(in) :: this
!        integer(kind=8), allocatable, dimension(..) :: array
!        integer :: err
!
!        allocate(array(this%shape()))
!        err = nf_get_vara_int64(this%ncid, this%varid, lbound(array), shape(array), array)
!        call check(err)
!    end subroutine
!
!    subroutine read_real32_alloc(this, array)
!        use netcdf4_f03
!        class(variable_t), intent(in) :: this
!        real(kind=4), allocatable, dimension(..) :: array
!        integer :: err
!
!        allocate(array(this%shape()))
!        err = nf_get_vara_real(this%ncid, this%varid, lbound(array), shape(array), array)
!        call check(err)
!    end subroutine
!
!    subroutine read_real64_alloc(this, array)
!        use netcdf4_f03
!        class(variable_t), intent(in) :: this
!        real(kind=8), allocatable, dimension(..) :: array
!        integer :: err
!
!        allocate(array(this%shape()))
!        err = nf_get_vara_double(this%ncid, this%varid, lbound(array), shape(array), array)
!        call check(err)
!    end subroutine

    subroutine check(err)
        use portable_mod
        use netcdf4_f03
        use mpi_f08
        integer, intent(in) :: err
        integer :: merr

        if (err /= 0) THEN
            write(*,*) nf_strerror(err)
            call traceback()
            call MPI_Abort(MPI_COMM_WORLD,err,merr)
        end if
    end subroutine
end module

max_stars_repo_name: ScottWales/farray

max_issues_repo_name: ScottWales/farray

max_forks_repo_name: ScottWales/farray


text:
c===========================================================================
c
c   This file is part of TISEAN
c
c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
c
c   TISEAN is free software; you can redistribute it and/or modify
c   it under the terms of the GNU General Public License as published by
c   the Free Software Foundation; either version 2 of the License, or
c   (at your option) any later version.
c
c   TISEAN is distributed in the hope that it will be useful,
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c   GNU General Public License for more details.
c
c   You should have received a copy of the GNU General Public License
c   along with TISEAN; if not, write to the Free Software
c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c
c===========================================================================
c   part of the TISEAN randomize package for constraint surrogates
c   cost function
c   spike train power spectrum
c   author T. Schreiber (1999)
c
c-------------------------------------------------------------------
c get cost function specific options
c
      subroutine opts_cost(ncol)
      parameter(mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight

      iweight=ican('W',0)
      fmax=fcan("F",0)
      nfreq=ican("#",0)
      inter=lopt("i",1)
      ncol=1
      end

c-------------------------------------------------------------------
c print version information on cost function
c
      subroutine what_cost()
      call ptext("Cost function: spike train power spectrum")
      end

c-------------------------------------------------------------------
c print cost function specific usage message
c
      subroutine usage_cost()
      call ptext("Cost function options: [-F# -## -w# -i]")
      call popt("W",
     .   "average: 0=max(s) 1=|s|/f 2=(s/f)**2 3=|s| (0)")
      call popt("F","maximal frequency (2*l / total time)")
      call popt("#","number of frequencies (F* total time /2)")
      call popt("w","frequency resolution (0)")
      call popt("i","expect intervals rather than times")
      end

c-------------------------------------------------------------------
c initialise all that is needed for cost function
c
      subroutine cost_init()
      parameter(nx=100000,mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      dimension x(nx)
      common nmax,cost,temp,cmin,rate,x

      if(fmax.le.0.) fmax=2*nmax/(x(nmax)-x(1))
      if(nfreq.le.0) nfreq=fmax*(x(nmax)-x(1))/2
      if(nfreq.gt.mfreq) write(istderr(),'(a)')
     .   "truncated to ", mfreq," frequencies"
      nfreq=min(mfreq,nfreq)
      write(istderr(),*) "randomize_spikespec: total time covered: ",
     .   x(nmax)-x(1)
      write(istderr(),*) "randomize_spikespec: computing ", nfreq,
     .   " frequencies up to ", fmax
      call sspect(nfreq,fmax/nfreq,sp0r,sp0i,sp0)
      end

c-------------------------------------------------------------------
c initial transformation on time series and its inverse
c
      subroutine cost_transform(nmax,mcmax,nxdum,x)
      parameter(mfreq=100000,nx=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      dimension x(nx), lx(nx)

      if(inter.eq.0) goto 1
      do 10 n=2,nmax
 10      x(n)=x(n)+x(n-1)
 1    call sort(nmax,x,lx)
      end

      subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
      parameter(mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      dimension x(nmax), y(nmax)

      do 10 n=1,nmax
 10      y(n)=x(n)
      if(inter.ne.1) return
      do 20 n=nmax,2,-1
 20      y(n)=y(n)-y(n-1)
      end

c-------------------------------------------------------------------
c compute full cost function from scratch
c
      subroutine cost_full(iv)
      parameter(mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      common nmax,cost

      call sspect(nfreq,fmax/nfreq,spr,spi,sp)

      cc=0
      do 10 i=1,nfreq
 10      call aver(cc,sp(i)-sp0(i),i)
      cost=cc
      if(iv.ne.0) call dump()
      end

c-------------------------------------------------------------------
c compute changed cost function on exchange of n1 and n2
c
      subroutine cost_update(n1,n2,cmax,iaccept,iv)
      parameter(mfreq=100000,nx=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      dimension sprcop(mfreq), spicop(mfreq), spcop(mfreq), x(nx)
      common nmax,cost,temp,cmin,rate,x
      data pi/3.1415926/

      comp=0
      iaccept=0
      do 10 i=1,nfreq
         f=i*(fmax/nfreq)
         omega=2*pi*f
         xx=x(n1-1)+x(n1+1)-x(n1)
         sprcop(i)=spr(i)-cos(omega*x(n1))+cos(omega*xx)
         spicop(i)=spi(i)-sin(omega*x(n1))+sin(omega*xx)
         spcop(i)=sprcop(i)**2+spicop(i)**2
         call aver(comp,sp0(i)-spcop(i),i)
 10      if(comp.ge.cmax) return
      cost=comp  ! if got here: accept
      iaccept=1
      call exch(n1,n2)
      if(iv.ne.0) call panic(spcop)
      do 20 i=1,nfreq
         spr(i)=sprcop(i)
         spi(i)=spicop(i)
 20      sp(i)=spcop(i)
      end

c-------------------------------------------------------------------
c compute spectrum from scratch
c
      subroutine sspect(nfreq,fres,spr,spi,sp)
      parameter(nx=100000)
      dimension spr(*), spi(*), sp(*), x(nx)
      common nmax,cost,temp,cmin,rate,x
      data pi/3.1415926/

      do 10 i=1,nfreq
         f=i*fres
         omega=2*pi*f
         sr=0
         si=0
         do 20 n=1,nmax
            sr=sr+cos(omega*x(n))
 20         si=si+sin(omega*x(n))
         spr(i)=sr
         spi(i)=si
 10      sp(i)=sr**2+si**2
      end

c-------------------------------------------------------------------
c weighted average of autocorrelation
c
      subroutine aver(cav,dc,n)
      parameter(mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight

      if(iweight.eq.0) then
         cav=max(cav,abs(dc))    ! max (L-infinity) norm
      else if(iweight.eq.1) then
         cav=cav+abs(dc)/n       ! L-1 norm (sum of moduli), weighted by freq.
      else if(iweight.eq.2) then
         cav=cav+(dc/n)**2       ! L-2 norm (sum of squares), weighted by freq.
      else if(iweight.eq.2) then
         cav=cav+abs(dc)         ! L-1 norm (sum of moduli)
      endif
      end

c-------------------------------------------------------------------
c diagnostic output
c
      subroutine dump()
      parameter(mfreq=100000)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      common nmax

      write(istderr(),'(5hgoal ,4f9.3)') (sp0(n),n=1,min(4,nfreq))
      write(istderr(),'(5his   ,4f9.3)') (sp(n),n=1,min(4,nfreq))
      write(istderr(),'(5hmiss ,4f9.3)')
     .   ((sp0(n)-sp(n)),n=1,min(4,nfreq))
      write(istderr(),'()')
      end

      subroutine panic(spcop)
      parameter(mfreq=100000)
      dimension spcop(*)
      dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
     .   sp0(mfreq), sp(mfreq)
      common /costcom/ nfreq, fmax, inter,
     .    sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
      common nmax

      call cost_full(0)
      write(istderr(),'(7hupdate ,4f9.3)')
     .   (spcop(n),n=1,min(4,nfreq))
      write(istderr(),'(7hfresh  ,4f9.3)') (sp(n),n=1,min(4,nfreq))
      write(istderr(),'(7hdiscr  ,4f9.3)')
     .   ((spcop(n)-sp(n)),n=1,min(4,nfreq))
      write(istderr(),'()')
      end

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# User: Set here the F90 compiler and options
#       Pedefined compilers: INTEL, PGF, HPUX, LAHEY
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

COMPILER = GNU
#COMPILER = LAHEY
#COMPILER = INTEL
#COMPILER = PGF
#COMPILER = HPUX

FC_GNU     = g95
FOPT_GNU   = -cpp -O -pg -fbounds-check
FC_LAHEY   = lf95
FOPT_LAHEY = -Cpp --pca
#FOPT_LAHEY = -Cpp --chk a,e,s,u --pca --ap -O0 -g --trap
FC_INTEL   = ifort
FOPT_INTEL = -cpp -O -mp -pc80 -prec_div -tpp7 -implicitnone
FC_PGF     = pgf90
FOPT_PGF   = -Mpreprocess -O -fast -pc 80 -Kieee
FC_HPUX    = f90
FOPT_HPUX  = -O -u +Oall +check=on

# define FULL_ALGEBRA for non-sparse integration
FC   = $(FC_$(COMPILER))
FOPT = $(FOPT_$(COMPILER)) # -DFULL_ALGEBRA

LIBS =
#LIBS = -llapack -lblas

# Command to create Matlab mex gateway routines
# Note: use $(FC) as the mex Fortran compiler
MEX  = mex

GENSRC = KPP_ROOT_Precision.f90  \
	 KPP_ROOT_Parameters.f90     \
	 KPP_ROOT_Global.f90

GENOBJ = KPP_ROOT_Precision.o    \
	 KPP_ROOT_Parameters.o       \
	 KPP_ROOT_Global.o

FUNSRC = KPP_ROOT_Function.f90
FUNOBJ = KPP_ROOT_Function.o

JACSRC = KPP_ROOT_JacobianSP.f90  KPP_ROOT_Jacobian.f90
JACOBJ = KPP_ROOT_JacobianSP.o    KPP_ROOT_Jacobian.o

HESSRC = KPP_ROOT_HessianSP.f90   KPP_ROOT_Hessian.f90
HESOBJ = KPP_ROOT_HessianSP.o     KPP_ROOT_Hessian.o

STMSRC = KPP_ROOT_StoichiomSP.f90 KPP_ROOT_Stoichiom.f90
STMOBJ = KPP_ROOT_StoichiomSP.o   KPP_ROOT_Stoichiom.o

UTLSRC = KPP_ROOT_Rates.f90 KPP_ROOT_Util.f90 KPP_ROOT_Monitor.f90
UTLOBJ = KPP_ROOT_Rates.o   KPP_ROOT_Util.o   KPP_ROOT_Monitor.o

LASRC  = KPP_ROOT_LinearAlgebra.f90
LAOBJ  = KPP_ROOT_LinearAlgebra.o

STOCHSRC = KPP_ROOT_Stochastic.f90
STOCHOBJ = KPP_ROOT_Stochastic.o

MAINSRC = KPP_ROOT_Main.f90   KPP_ROOT_Initialize.f90   KPP_ROOT_Integrator.f90 KPP_ROOT_Model.f90
MAINOBJ = KPP_ROOT_Main.o     KPP_ROOT_Initialize.o     KPP_ROOT_Integrator.o   KPP_ROOT_Model.o

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# User: modify the line below to include only the
#       objects needed by your application
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ALLOBJ = $(GENOBJ) $(FUNOBJ) $(JACOBJ) $(HESOBJ) $(STMOBJ) \
	 $(UTLOBJ) $(LAOBJ)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# User: modify the line below to include only the
#       executables needed by your application
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
all:    exe

exe:	$(ALLOBJ) $(MAINOBJ)
	$(FC) $(FOPT) $(ALLOBJ) $(MAINOBJ) $(LIBS) -o KPP_ROOT.exe

stochastic:$(ALLOBJ) $(STOCHOBJ) $(MAINOBJ)
	$(FC) $(FOPT) $(ALLOBJ) $(STOCHOBJ) $(MAINOBJ) $(LIBS) \
	-o KPP_ROOT_stochastic.exe

mex:    $(ALLOBJ)
	$(MEX) FC#$(FC) -fortran -O KPP_ROOT_mex_Fun.f90     $(ALLOBJ)
	$(MEX) FC#$(FC) -fortran -O KPP_ROOT_mex_Jac_SP.f90  $(ALLOBJ)
	$(MEX) FC#$(FC) -fortran -O KPP_ROOT_mex_Hessian.f90 $(ALLOBJ)

clean:
	rm -f KPP_ROOT*.o KPP_ROOT*.mod \
	KPP_ROOT*.dat KPP_ROOT.exe KPP_ROOT*.mexglx \
	KPP_ROOT.map

distclean:
	rm -f KPP_ROOT*.o KPP_ROOT*.mod \
	KPP_ROOT*.dat KPP_ROOT.exe KPP_ROOT.map \
	KPP_ROOT*.f90 KPP_ROOT_*.mexglx

KPP_ROOT_Precision.o: KPP_ROOT_Precision.f90
	$(FC) $(FOPT) -c $<

KPP_ROOT_Parameters.o: KPP_ROOT_Parameters.f90 \
	            KPP_ROOT_Precision.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Monitor.o: KPP_ROOT_Monitor.f90 \
	             KPP_ROOT_Precision.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Global.o: KPP_ROOT_Global.f90 \
	            KPP_ROOT_Parameters.o KPP_ROOT_Precision.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Initialize.o: KPP_ROOT_Initialize.f90  $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_Function.o: KPP_ROOT_Function.f90  $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_Stochastic.o: KPP_ROOT_Stochastic.f90  $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_JacobianSP.o: KPP_ROOT_JacobianSP.f90 $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_Jacobian.o: KPP_ROOT_Jacobian.f90  $(GENOBJ) KPP_ROOT_JacobianSP.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_LinearAlgebra.o: KPP_ROOT_LinearAlgebra.f90 $(GENOBJ) KPP_ROOT_JacobianSP.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Rates.o: KPP_ROOT_Rates.f90  $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_HessianSP.o: KPP_ROOT_HessianSP.f90  $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_Hessian.o:  KPP_ROOT_Hessian.f90 $(GENOBJ) KPP_ROOT_HessianSP.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_StoichiomSP.o: KPP_ROOT_StoichiomSP.f90 $(GENOBJ)
	$(FC) $(FOPT) -c $<

KPP_ROOT_Stoichiom.o: KPP_ROOT_Stoichiom.f90  $(GENOBJ) KPP_ROOT_StoichiomSP.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Util.o: KPP_ROOT_Util.f90  $(GENOBJ) KPP_ROOT_Monitor.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Main.o: KPP_ROOT_Main.f90  $(ALLOBJ) KPP_ROOT_Initialize.o KPP_ROOT_Model.o KPP_ROOT_Integrator.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Model.o: KPP_ROOT_Model.f90  $(ALLOBJ) KPP_ROOT_Integrator.o
	$(FC) $(FOPT) -c $<

KPP_ROOT_Integrator.o: KPP_ROOT_Integrator.f90  $(ALLOBJ)
	$(FC) $(FOPT) -c $<

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
{ Operator tests }
"Test 0"
1 1=$["= Test 0 OK"]?~["= Test 0 failed"]?
"Test 1"
1 2=$["= Test 1 failed"]?~["= Test 1 OK"]?
"Test 2"
1 1>$["> Test 0 failed"]?~["> Test 0 OK"]?
2 1>$["> Test 1 OK"]?~["> Test 1 OK"]?
50a:a;0>a;99>~&$["~& Test 0 OK"]?~["~& Test 1 failed"]?
7 8 9 2(
.",".",".",".
1 2 +

max_stars_repo_name: afgbeveridge/EsotericInterpreters

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: afgbeveridge/EsotericInterpreters

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: afgbeveridge/EsotericInterpreters

max_forks_repo_licenses: MIT Unlicense


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine error_norm(rms)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     this function computes the norm of the difference between the
c     computed solution and the exact solution
c---------------------------------------------------------------------

      include 'header.h'

      integer i, j, k, m, d
      double precision xi, eta, zeta, u_exact(5), rms(5), add
      double precision rms_local(5)

      do m = 1, 5
         rms(m) = 0.0d0
      enddo

!$omp parallel default(shared)
!$omp& private(i,j,k,m,zeta,eta,xi,add,u_exact,rms_local)
!$omp&        shared(rms)
      do m = 1, 5
         rms_local(m) = 0.0d0
      enddo
!$omp do
      do k = 0, grid_points(3)-1
         zeta = dble(k) * dnzm1
         do j = 0, grid_points(2)-1
            eta = dble(j) * dnym1
            do i = 0, grid_points(1)-1
               xi = dble(i) * dnxm1
               call exact_solution(xi, eta, zeta, u_exact)

               do m = 1, 5
                  add = u(m,i,j,k)-u_exact(m)
                  rms_local(m) = rms_local(m) + add*add
               enddo
            enddo
         enddo
      enddo
!$omp end do nowait
      do m = 1, 5
!$omp critical
         rms(m) = rms(m) + rms_local(m)
!$omp end critical
      enddo
!$omp end parallel

      do m = 1, 5
         do d = 1, 3
            rms(m) = rms(m) / dble(grid_points(d)-2)
         enddo
         rms(m) = dsqrt(rms(m))
      enddo

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine rhs_norm(rms)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      include 'header.h'

      integer i, j, k, d, m
      double precision rms(5), add
      double precision rms_local(5)

      do m = 1, 5
         rms(m) = 0.0d0
      enddo

!$omp parallel default(shared) private(i,j,k,m,add,rms_local)
!$omp&        shared(rms)
      do m = 1, 5
         rms_local(m) = 0.0d0
      enddo
!$omp do
      do k = 1, grid_points(3)-2
         do j = 1, grid_points(2)-2
            do i = 1, grid_points(1)-2
               do m = 1, 5
                  add = rhs(m,i,j,k)
                  rms_local(m) = rms_local(m) + add*add
               enddo
            enddo
         enddo
      enddo
!$omp end do nowait
      do m = 1, 5
!$omp critical
         rms(m) = rms(m) + rms_local(m)
!$omp end critical
      enddo
!$omp end parallel

      do m = 1, 5
         do d = 1, 3
            rms(m) = rms(m) / dble(grid_points(d)-2)
         enddo
         rms(m) = dsqrt(rms(m))
      enddo

      return
      end


text:
      INCLUDE 'VICMAIN_FOR'
C
C Program signal
C
C 28 JUL 1997 ...RRD...      Ported to UNIX.
C 17 NOV 1995 ...JRY...      Written
C
      SUBROUTINE MAIN44
      INTEGER*2 DN
      INTEGER*4 CNT, DEF, ICNT, IND, ISTAT, NI, NL, NS
      INTEGER*4 EXPOS(50), INSTANCES(30), IUNIT(50), LS(2), PICSCALE(50)
      REAL*4    SIGNAL, REXPOS
      CHARACTER*255 LFNAME, TBL
      CHARACTER*8  TASKS(30)
      CHARACTER*1 TAB

      TAB = CHAR(9)
      CALL IFMESSAGE('signal version 28-Jul-97')
      CALL XVPCNT('INP',NI)
      CALL XVP('LIST',LFNAME,ICNT)
      CALL XVPARM('TBL',TBL,CNT,DEF,0)
      CALL XVPARM('LS',LS,CNT,DEF,0)

      IF (NI .NE. 0) THEN
         DO I=1,NI
            CALL XVUNIT(IUNIT(I),'INP',I,STAT,' ')
            CALL XVOPEN(IUNIT(I),STAT,'OPEN_ACT','SA','IO_ACT','SA',
     &                  'U_FORMAT','HALF',' ')
         ENDDO
      ELSE IF (ICNT .NE. 0) THEN    !input frames are in SRCH-list
     	 OPEN(UNIT=99,FILE=LFNAME,STATUS='OLD',ERR=999)
	 READ(99,FMT=1) LFNAME                      !SKIP FIRST LINE
1        FORMAT(A)
 	 DO I=1,51
	    READ(99,FMT=1,END=11,ERR=999) LFNAME
	    IF (I .GE. 51) GO TO 990
            NI = I
	    CALL XVUNIT(IUNIT(I),'NONE',I,ISTAT,'U_NAME',LFNAME,' ')
	    CALL XVOPEN(IUNIT(I),ISTAT,'OPEN_ACT','SA','IO_ACT','SA',
     &                  'U_FORMAT','HALF',' ')
         ENDDO
11       CONTINUE
      ELSE
         CALL XVMESSAGE('Specify input files as INP or LIST',' ')
         CALL ABEND
      END IF

C     Get PICSCALE from file label (number of frames PICSUM'd for each input)

      CALL XVGET(IUNIT(1),IND,'NL',NL,'NS',NS,' ')
      IF (LS(1) .GT. NL .OR. LS(2) .GT. NS) THEN
         CALL XVMESSAGE('LS EXCEEDS SIZE OF IMAGE', ' ')
         CALL ABEND
      END IF
      DO 5 I=1,NI
         CNT = 30
         CALL XLGET(IUNIT(I),'PROPERTY','EXPOSURE_DURATION',REXPOS,
     &              ISTAT,'PROPERTY','CASSINI-ISS','FORMAT','REAL',' ')
         IF (ISTAT .NE. 1) THEN
            CALL MABEND('ERROR READING LABEL EXPOSURE_DURATION')
         END IF
         EXPOS(I) = IFIX(REXPOS)
         CALL XLHINFO(IUNIT(I),TASKS,INSTANCES,CNT,IND,' ')
         DO J=CNT,1,-1		!Search for last value of picture scale
            CALL XLGET(IUNIT(I),'HISTORY','PICSCALE',PICSCALE(I),IND,
     &            'HIST',TASKS(J),'INSTANCE',INSTANCES(J),
     &            'FORMAT','INT',' ')
            IF (IND.EQ.1) GOTO 5
         ENDDO

C Didn't find a PICSCALE
         PICSCALE(I) = 1		!Default picture scale
5     CONTINUE

C Open output table
      OPEN(15,FILE=TBL,STATUS='UNKNOWN',IOSTAT=JST,ERR=992)
      WRITE(15,100) ' EXPOSURE',TAB,' SIGNAL'
100   FORMAT(A9,A1,A7)

      DO I=1,NI 	!Read data from each flat-field frame...
         CALL XVREAD(IUNIT(I),DN,IND,'LINE',LS(1),'SAMP',LS(2),
     &               'NSAMPS',1,' ')
         SIGNAL = FLOAT(DN)/PICSCALE(I)
         WRITE(15,110) EXPOS(I),TAB,SIGNAL
      ENDDO
110   FORMAT(I10,A1,F8.2)

C Close inputs and output
      DO I=1,NI
        CALL XVCLOSE(IUNIT(I),IND,' ')
      ENDDO
      CLOSE(15)

      CALL XVMESSAGE('signal task completed',' ')
      RETURN

  992 CALL XVMESSAGE('ERROR WRITING OUTPUT TABLE', ' ')
      CALL ABEND
  990 CALL XVMESSAGE('MORE THAN 50 FILENAMES IN LIST', ' ')
      CALL ABEND
  999 CALL XVMESSAGE('ERROR OPENING INPUT SRCH LIST FILE', ' ')
      CALL ABEND
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
! Note that the () will be introduced by the unparser
subroutine where_example()
!    logical :: should_set_cell(nx,ny,nz)
!    logical(logical_kind) :: should_set_cell(nx,ny,nz)
     logical(2) :: i(50)
!    integer :: nx,ny,nz
!    integer :: should_set_cell(nx,ny,nz)
end subroutine where_example

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


subroutine attryn(a, c)
!(LICENSE:PD)
! @(#) test termattrs(3f) and query which attributes this terminal type takes
   use M_ncurses
   implicit none
   integer(C_LONG)       :: a     !chtype a
   integer(C_LONG)       :: c     !chtype c
   integer               :: ierr
   !if( ior(a,c)  .eq. a) then
   !if( iand(a,c)  .ne. 0) then  !! this would get A_NORMAL wrong, because it is zero
   if( iand(a,c) .eq. c) then
      ierr=addstr("Yes "//C_NULL_CHAR)
      ierr=attron(c)
      ierr=addstr("abcdefghijklmnopqrstuvwxyz0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"//C_NULL_CHAR)
      ierr=attroff(c)
   else
      ierr=addstr("No  "//C_NULL_CHAR)
      ierr=attron(c)
      ierr=addstr("abcdefghijklmnopqrstuvwxyz0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"//C_NULL_CHAR)
      ierr=attroff(c)
   endif
   ierr=addch(int(ichar(C_NEW_LINE),C_LONG))
end subroutine attryn

program seeattr
   use M_ncurses
   use ISO_FORTRAN_ENV
   integer(C_LONG)       :: attributes !chtype attributes
   integer(C_INT)        :: y, x
   integer               :: ierr
   integer(C_INT)        :: i
   integer(kind=mmask_t) :: mmask
   stdscr=initscr()
   !if(stdscr == C_NULL_PTR)then
   if (.not.c_associated(stdscr)) then
      ierr=endwin()
      write(*,*)'could not open standard screen. input not a tty? ERRNO=',ierr
      stop
   endif
   ierr=raw()
   ierr=cbreak()    ! Line buffering disabled, Pass on everthing to me
   ierr=noecho()    ! only put on screen what is explicitly placed there, do not echo input
   attributes = termattrs()
   ierr=addstr("This terminal is capable of the following attributes:"//C_NEW_LINE//C_NULL_CHAR)
   ierr=printw("%13s"//C_NULL_CHAR, "AltCharSet: "//C_NULL_CHAR); call attryn(attributes, A_ALTCHARSET)
   ierr=printw("%13s"//C_NULL_CHAR, "Blink: "//C_NULL_CHAR);      call attryn(attributes, A_BLINK)
   ierr=printw("%13s"//C_NULL_CHAR, "Bold: "//C_NULL_CHAR);       call attryn(attributes, A_BOLD)
   ierr=printw("%13s"//C_NULL_CHAR, "Dim: "//C_NULL_CHAR);        call attryn(attributes, A_DIM)
   ierr=printw("%13s"//C_NULL_CHAR, "Invis(ible):"//C_NULL_CHAR); call attryn(attributes, A_INVIS)
   ierr=printw("%13s"//C_NULL_CHAR, "Normal: "//C_NULL_CHAR);     call attryn(attributes, A_NORMAL)
   ierr=printw("%13s"//C_NULL_CHAR, "Reverse: "//C_NULL_CHAR);    call attryn(attributes, A_REVERSE)
   ierr=printw("%13s"//C_NULL_CHAR, "Standout: "//C_NULL_CHAR);   call attryn(attributes, A_STANDOUT)
   ierr=printw("%13s"//C_NULL_CHAR, "Underline: "//C_NULL_CHAR);  call attryn(attributes, A_UNDERLINE)
   ierr=printw("%13s"//C_NULL_CHAR, "Protect: "//C_NULL_CHAR);    call attryn(attributes, A_PROTECT)
   ierr=printw("%13s"//C_NULL_CHAR, "Horizontal: "//C_NULL_CHAR); call attryn(attributes, A_HORIZONTAL)
   ierr=printw("%13s"//C_NULL_CHAR, "Left: "//C_NULL_CHAR);       call attryn(attributes, A_LEFT)
   ierr=printw("%13s"//C_NULL_CHAR, "Low: "//C_NULL_CHAR);        call attryn(attributes, A_LOW)
   ierr=printw("%13s"//C_NULL_CHAR, "Right: "//C_NULL_CHAR);      call attryn(attributes, A_RIGHT)
   ierr=printw("%13s"//C_NULL_CHAR, "Top: "//C_NULL_CHAR);        call attryn(attributes, A_TOP)
   ierr=printw("%13s"//C_NULL_CHAR, "Vertical: "//C_NULL_CHAR);   call attryn(attributes, A_VERTICAL)
   ierr=printw("%13s"//C_NULL_CHAR, "Italic:   "//C_NULL_CHAR);   call attryn(attributes, A_ITALIC)
   !ierr=printw("%13s"//C_NULL_CHAR, "Color:    "//C_NULL_CHAR);   call attryn(attributes, A_COLOR)
   !ierr=printw("%13s"//C_NULL_CHAR, "Attributes:"//C_NULL_CHAR);  call attryn(attributes, A_ATTRIBUTES)
   !ierr=printw("%13s"//C_NULL_CHAR, "CharText:  "//C_NULL_CHAR);  call attryn(attributes, A_CHARTEXT)
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=clear()               ! clear stdscr
   ierr=box(stdscr,0_C_LONG,0_C_LONG)
   ierr=move(0_C_INT,0_C_INT) ! go to top of page
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   call getmaxyx(stdscr,y,x)
!-----------------------------------------------------------------------------------------------------------------------------------
   do i=0,y-1
      ierr=move(i,0)
      ierr=hline(0_C_LONG,x)
   enddo
   ierr=move(0,0)
   ierr=printw("Window size is %d rows, %d columns."//C_NEW_LINE//C_NULL_CHAR,y,x)
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=addstr("This terminal "//C_NULL_CHAR)
   if(has_ic())then
      ierr=addstr("has "//C_NULL_CHAR)
   else
      ierr=addstr("does not have "//C_NULL_CHAR)
   endif
   ierr=addstr("insert/delete character abilities"//C_NEW_LINE//C_NULL_CHAR)
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=addstr("This terminal "//C_NULL_CHAR)
   if(has_il())then
      ierr=addstr("has "//C_NULL_CHAR)
   else
      ierr=addstr("does not have "//C_NULL_CHAR)
   endif
   ierr=addstr("insert/delete line abilities"//C_NEW_LINE//C_NULL_CHAR)
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   if(has_colors())then
      ierr=addstr("This terminal can do colors."//C_NEW_LINE//C_NULL_CHAR)
      if(start_color() .ne. OK) then ! next attempt to initialize curses colors
         ierr=printw("Unable to start colors."//C_NEW_LINE//C_NULL_CHAR)
      else  ! colors are okay; continue
         ierr=printw("Colors have been properly initialized."//C_NEW_LINE//C_NULL_CHAR)
         call getcolor(COLORS,COLOR_PAIRS) !! extension to get these set
         ierr=printw("NCurses reports that you can use %d colors, "//C_NULL_CHAR,COLORS)
         ierr=printw("and %d color pairs."//C_NEW_LINE//C_NULL_CHAR,COLOR_PAIRS)
         ierr=refresh()
         ierr=getch()
         ierr=init_pair(1_C_SHORT, COLOR_BLACK,     COLOR_RED);
         ierr=init_pair(2_C_SHORT, COLOR_RED,       COLOR_GREEN);
         ierr=init_pair(3_C_SHORT, COLOR_GREEN,     COLOR_BLACK);
         ierr=attron(COLOR_PAIR(1))
         ierr=addstr("COLOR_PAIR(1)"//C_NULL_CHAR)
         ierr=attroff(COLOR_PAIR(1))
         ierr=attron(COLOR_PAIR(2))
         ierr=addstr("COLOR_PAIR(2)"//C_NULL_CHAR)
         ierr=attroff(COLOR_PAIR(2))
         ierr=attron(COLOR_PAIR(3))
         ierr=addstr("COLOR_PAIR(3)"//C_NEW_LINE//C_NULL_CHAR)
         ierr=attroff(COLOR_PAIR(3))
      endif
   else
      ierr=addstr("This terminal cannot do colors."//C_NEW_LINE//C_NULL_CHAR)
   endif
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   if(can_change_color())then
      ierr=addstr("This terminal can change the standard colors."//C_NEW_LINE//C_NULL_CHAR)
   else
      ierr=addstr("This terminal cannot change the colors."//C_NEW_LINE//C_NULL_CHAR)
   endif
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=printw("This terminal's baud rate is %d."//C_NEW_LINE//C_NULL_CHAR,baudrate())
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=keypad(stdscr,TRUE)
   mmask = mousemask(ALL_MOUSE_EVENTS,C_NULL_PTR)
   if(mmask == 0)then
      ierr=addstr("Unable to access the mouse on this terminal."//C_NEW_LINE//C_NULL_CHAR)
   else
      ierr=addstr("Mouse events can be captured."//C_NEW_LINE//C_NULL_CHAR)
   endif
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=addstr("VERSION:"//trim(curses_version())//C_NEW_LINE//C_NULL_CHAR)
   ierr=refresh()
   ierr=getch()
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=addstr("TERMINAL TYPE:"//trim(termname())//C_NEW_LINE//C_NULL_CHAR)
   ierr=addstr("LONG TERMINAL TYPE:"//trim(longname())//C_NEW_LINE//C_NULL_CHAR)
   ierr=refresh()
   ierr=getch()
!-------------------------------------------------------------------------------
   ierr=addstr("COMPILER VERSION:   "//compiler_version()//C_NEW_LINE//C_NULL_CHAR)
   ierr=addstr("COMPILER OPTIONS:   "//compiler_options()//C_NEW_LINE//C_NULL_CHAR)
!-------------------------------------------------------------------------------
   ierr=printw("ERASECHAR:          %s"//C_NEW_LINE//C_NULL_CHAR,keyname(ichar(erasechar())))
   ierr=printw("KILLCHAR:           %s"//C_NEW_LINE//C_NULL_CHAR,keyname(ichar(killchar())))
!-----------------------------------------------------------------------------------------------------------------------------------
   do i=1,3
      ierr=getch()
   enddo
!-----------------------------------------------------------------------------------------------------------------------------------
   ierr=endwin()
!-------------------------------------------------------------------------------
! CDEF: char *keybound ( int keycode, int count );
! CDEF: char *slk_label ( int labnum );
! CDEF: char *unctrl ( chtype c );
!-----------------------------------------------------------------------------------------------------------------------------------
end program seeattr


text:
c
c -------------------------------------------------------------
c
       subroutine qad(valbig,mitot,nvar,
     .                svdflx,qc1d,lenbc,lratiox,hx,
     .                maux,aux,auxc1d,delt,mptr)

      use amr_module
       implicit double precision (a-h, o-z)

       logical qprint

       dimension valbig(nvar,mitot)
       dimension qc1d(nvar,lenbc)
       dimension svdflx(nvar,lenbc)
       dimension aux(maux,mitot)
       dimension auxc1d(maux,lenbc)

c
c ::::::::::::::::::::::::::: QAD ::::::::::::::::::::::::::::::::::
c  solve RP between ghost cell value on fine grid and coarse grid
c  value that ghost cell overlaps. The resulting fluctuations
c  are added in to coarse grid value, as a conservation fixup.
c  Done each fine grid time step. If source terms are present, the
c  coarse grid value is advanced by source terms each fine time step too.

c  No change needed in this sub. for spherical mapping: correctly
c  mapped vals already in bcs on this fine grid and coarse saved
c  vals also properly prepared
c
c Side 1 is the left side of the fine grid patch.
c Side 2 is the right side of the fine grid patch.
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c      # local storage
c      # note that dimension here are bigger than dimensions used
c      # in rp1, but shouldn't matter since wave is not used in qad
c      # and for other arrays it is only the last parameter that is wrong
c      #  ok as long as meqn, mwaves < maxvar

       parameter (max1dp1 = max1d+1)
       dimension ql(nvar,max1dp1),    qr(nvar,max1dp1)
       dimension wave(nvar,mwaves,max1dp1), s(mwaves,max1dp1)
       dimension amdq(nvar,max1dp1),  apdq(nvar,max1dp1)
       dimension auxl(maxaux*max1dp1),  auxr(maxaux*max1dp1)
c
c  WARNING: auxl,auxr dimensioned at max possible, but used as if
c  they were dimensioned as the real maux by max1dp1. Would be better
c  of course to dimension by maux by max1dp1 but this wont work if maux=0
c  So need to access using your own indexing into auxl,auxr.
       iaddaux(iaux,i) = iaux + maux*(i-1)

       data qprint/.false./
c
c      aux is auxiliary array with user parameters needed in Riemann solvers
c          on fine grid corresponding to valbig
c      auxc1d is coarse grid stuff from around boundary, same format as qc1d
c      auxl, auxr are work arrays needed to pass stuff to rp1
c      maux is the number of aux variables, which may be zero.
c

       tgrid = rnode(timemult, mptr)
       if (qprint)
     .      write(dbugunit,*)" working on grid ",mptr," time ",tgrid
       level = node(nestlevel, mptr)
       index = 0

c
c--------
c  side 1
c--------
c
       if (maux.gt.0) then
          do 5 ma = 1,maux
             if (auxtype(ma).eq."xleft") then
c                # Assuming velocity at left-face, this fix
c                # preserves conservation in incompressible flow:
                 auxl(iaddaux(ma,2)) = aux(ma,nghost+1)
               else
c                # Normal case -- we set the aux arrays
c                # from the cell corresponding  to q
                 auxl(iaddaux(ma,2)) = aux(ma,nghost)
               endif
  5          continue
          endif
       do 10 ivar = 1, nvar
         ql(ivar,2) = valbig(ivar,nghost)
 10    continue

       index = index + 1
       lind = 1
       if (maux.gt.0) then
         do 24 ma=1,maux
            auxr(iaddaux(ma,lind)) = auxc1d(ma,index)
   24       continue
       endif
       do 25 ivar = 1, nvar
         qr(ivar,lind) = qc1d(ivar,index)
 25    continue

       if (qprint) then
         write(dbugunit,*) 'side 1, ql and qr:'
            write(dbugunit,4101) i,qr(1,1),ql(1,2)
 4101      format(i3,4e16.6)
         if (maux .gt. 0) then
             write(dbugunit,*) 'side 1, auxr:'
             write(dbugunit,4101) i,(auxr(iaddaux(ma,1)),ma=1,maux)
             write(dbugunit,*) 'side 1, auxl:'
             write(dbugunit,4101) i,(auxl(iaddaux(ma,2)),ma=1,maux)
         endif
       endif

       call rp1(max1dp1-2*nghost,nvar,mwaves,maux,nghost,
     .              2-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)

c
c we have the wave. for side 1 add into sdflxm
c
       influx  = 1
       do 40 ivar = 1, nvar
           svdflx(ivar,influx) = svdflx(ivar,influx)
     .                     + amdq(ivar,2) * delt
     .                     + apdq(ivar,2) * delt
 40    continue

c--------
c  side 2
c--------
c
        if (maux.gt.0) then
          do 205 ma = 1,maux
             auxr(iaddaux(ma,1)) = aux(ma,mitot-nghost+1)
 205         continue
          endif
        do 210 ivar = 1, nvar
          qr(ivar,1) = valbig(ivar,mitot-nghost+1)
 210      continue

         index = index + 1
         lind = 1
         if (maux.gt.0) then
            do 224 ma=1,maux
             if (auxtype(ma).eq."xleft") then
c                # Assuming velocity at left-face, this fix
c                # preserves conservation in incompressible flow:
                 auxl(iaddaux(ma,lind+1)) = aux(ma,mitot-nghost+1)
               else
                 auxl(iaddaux(ma,lind+1)) = auxc1d(ma,index)
               endif
  224          continue
            endif
         do 225 ivar = 1, nvar
             ql(ivar,lind+1) = qc1d(ivar,index)
 225     continue

       if (qprint) then
         write(dbugunit,*) 'side 2, ql and qr:'
            write(dbugunit,4101) i,ql(1,2),qr(1,1)
       endif

       call rp1(max1dp1-2*nghost,nvar,mwaves,maux,nghost,
     .              2-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)

c
c we have the wave. for side 2 add into sdflxp
C
          influx  = influx + 1
          do 340 ivar = 1, nvar
              svdflx(ivar,influx) = svdflx(ivar,influx)
     .                     - amdq(ivar,2) * delt
     .                     - apdq(ivar,2) * delt
 340       continue

c      # for source terms:
       if (method(4) .ne. 0) then   ! should I test here if index=0 and all skipped?
           call src1d(nvar,nghost,lenbc,qc1d,maux,auxc1d,tgrid,delt)
c      # how can this be right - where is the integrated src term used?
           endif

       return
       end

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
       subroutine u_ei

! --------------------------------------------------- !
!           F. Califano, February 2006                !
!                                                     !
!     Calcolo delle velocita' degli elettroni e ioni  !
!                                                     !
!     u_e = U - j / n;    u_i = U + de^2 j / n        !
!          ( de^2 = me / mi )             !           !
!                                                     !
! --------------------------------------------------- !

!**************************************************
!  MPI PARALLEL VERSION: VALENTINI-FAGANELLO 2009
!**************************************************

use box_mod
use parameter_mod
use parallel_mod
use dom_distr_mod
use fields_DP_mod
use fields_UJ_mod

IMPLICIT NONE

integer :: ix, iy, iz

! expressions without smoothing
!Uex = Ux - de12inv * Jx * Dinv
!Uey = Uy - de12inv * Jy * Dinv
!Uez = Uz - de12inv * Jz * Dinv
!Uix = Ux + de12inv * de2 * Jx * Dinv
!Uiy = Uy + de12inv * de2 * Jy * Dinv
!Uiz = Uz + de12inv * de2 * Jz * Dinv

! smoothing to ideal-MHD towards boundaries
do iz = 1, nzl
 do iy = 1, nyl
  do ix = 1, nxl
    Uex(ix,iy,iz) = Ux(ix,iy,iz) - de12inv * &
     lambda(ixlg + ix - 1) * Jx(ix,iy,iz) * Dinv(ix,iy,iz)
    Uey(ix,iy,iz) = Uy(ix,iy,iz) - de12inv * &
     lambda(ixlg + ix - 1) * Jy(ix,iy,iz) * Dinv(ix,iy,iz)
    Uez(ix,iy,iz) = Uz(ix,iy,iz) - de12inv * &
     lambda(ixlg + ix - 1) * Jz(ix,iy,iz) * Dinv(ix,iy,iz)
    Uix(ix,iy,iz) = Ux(ix,iy,iz) + de12inv * de2 * &
     lambda(ixlg + ix - 1) * Jx(ix,iy,iz) * Dinv(ix,iy,iz)
    Uiy(ix,iy,iz) = Uy(ix,iy,iz) + de12inv * de2 * &
     lambda(ixlg + ix - 1) * Jy(ix,iy,iz) * Dinv(ix,iy,iz)
    Uiz(ix,iy,iz) = Uz(ix,iy,iz) + de12inv * de2 * &
     lambda(ixlg + ix - 1) * Jz(ix,iy,iz) * Dinv(ix,iy,iz)
  enddo
 enddo
enddo

end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text: 20081228 11:30:37 nbsp  Welcome to the Wiki. Please read Welcome to the Wiki/Business Owner; it will help explain how you can make the wiki a positive experience for you without clashing with established wiki social norms. Please consider the importance of using your RealName, you wouldnt walk into a room an introduce yourself as davistext. Users/JasonAller


! { dg-do compile { target { ! *-*-* } } }
! SKIP THIS FILE
!
! Used by codimension_2.f90
!
! Additional file to check that using the module doesn't generate
! a token symbol. (The module is also used by codimension_2.f90.)
!
subroutine ttest
  use global_coarrays
  implicit none
  b(:) = b(:)[2]
end


text:
C
C	$Id: conbd.f,v 1.6 2008-07-27 00:23:01 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE CONBD
C
C Calling this do-nothing subroutine forces "ld" to load the following
C block data routine (but only if they are in the same ".f" file).
C
        RETURN
C
      END
CNOSPLIT
      BLOCKDATA CONBDX
      COMMON /CONRE1/ IOFFP      ,SPVAL
      COMMON /CONRE4/ ISIZEL     ,ISIZEM     ,ISIZEP     ,NREP       ,
     1                NCRT       ,ILAB       ,NULBLL     ,IOFFD      ,
     2                EXT        ,IOFFM      ,ISOLID     ,NLA        ,
     3                NLM        ,XLT        ,YBT        ,SIDE
      DATA IOFFP,SPVAL/0,0.0/
      DATA ISIZEM,ISIZEP,NLA,NLM,XLT,YBT,SIDE,ISOLID/
     1        2,     1,   16, 40,.05,.05,   .9, 1023/
      DATA ISIZEL,NREP,NCRT/
     1       1,    6,   4 /
      DATA EXT,IOFFD,NULBLL,IOFFM,ILAB/.25,0,3,0,1/
C
C REVISION HISTORY---
C
C JANUARY 1978     DELETED REFERENCES TO THE  *COSY  CARDS AND
C                  ADDED REVISION HISTORY
C FEBURARY 1979    MADE CODE CONFORM TO FORTRAN 66 STANDARDS
C JANUARY 1980     CHANGED LIBRARY NAME FROM NSSL TO
C                  PORTLIB FOR MOVE TO PORTLIB
C MAY 1980         CONRECQCK ON PORTLIB WAS DISCOVERED TO BE ESSENTIALLY
C                  THE WRONG VERSION.  EXTENSIVE MODIFICATIONS TO CLGEN,
C                  MAXMIN, CONREC TO MAKE IT PORTABLE (INVOLVED REMOVAL
C                  OF PWRT CALLS, ENCODE STATEMENTS, INAPPROPRIATE CON-
C                  STANTS, ETC.).
C JUNE 1984        UPDATED TO FORTRAN 77 AND TO GKS
C MARCH 1990       CORRECTED THE USE OF SET CALLS
C-----------------------------------------------------------------------
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      subroutine cudaclaw_bc2_default(maxmx,maxmy,meqn,mbc,
     &      mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt,mthbc)
      implicit none

      integer maxmx, maxmy, meqn, mbc, mx, my, maux, mthbc(4)
      double precision xlower, ylower, dx, dy, t, dt

      double precision q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      double precision aux(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, maux)

      integer m, i, j, ibc, jbc

c
c-------------------------------------------------------
c     # left boundary:
c-------------------------------------------------------
      go to (100,110,120,130) mthbc(1)+1
c     this is how we skip over this side... if (mthbc(1)+1
c     is not 1,2,3 or 4, then the goto above falls through to here...
      goto 199

  100 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(1)=0 and no BCs specified in bc2'
      stop
      go to 199
c
  110 continue
c     # zero-order extrapolation:
      do 115 m=1,meqn
         do 115 ibc=1,mbc
            do 115 j = 1-mbc,my+mbc
               q(1-ibc,j,m) = q(1,j,m)
  115       continue
      go to 199

  120 continue
c     # periodic:
      do 125 m=1,meqn
         do 125 ibc=1,mbc
            do 125 j = 1-mbc, my+mbc
               q(1-ibc,j,m) = q(mx+1-ibc,j,m)
  125       continue
      go to 199

  130 continue
c     # solid wall (assumes 2'nd component is velocity or momentum in x):
      do 135 m=1,meqn
         do 135 ibc=1,mbc
            do 135 j = 1-mbc, my+mbc
               q(1-ibc,j,m) = q(ibc,j,m)
  135       continue
c     # negate the normal velocity:
      do 136 ibc=1,mbc
         do 136 j = 1-mbc, my+mbc
            q(1-ibc,j,2) = -q(ibc,j,2)
  136    continue
      go to 199

  199 continue
c
c-------------------------------------------------------
c     # right boundary:
c-------------------------------------------------------
      go to (200,210,220,230) mthbc(2)+1
      goto 299
c
  200 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2'
      stop
      go to 299

  210 continue
c     # zero-order extrapolation:
      do 215 m=1,meqn
         do 215 ibc=1,mbc
            do 215 j = 1-mbc,my+mbc
               q(mx+ibc,j,m) = q(mx,j,m)
  215       continue
      go to 299

  220 continue
c     # periodic:
      do 225 m=1,meqn
         do 225 ibc=1,mbc
            do 225 j = 1-mbc, my+mbc
               q(mx+ibc,j,m) = q(ibc,j,m)
  225       continue
      go to 299

  230 continue
c     # solid wall (assumes 2'nd component is velocity or momentum in x):
      do 235 m=1,meqn
         do 235 ibc=1,mbc
            do 235 j = 1-mbc, my+mbc
               q(mx+ibc,j,m) = q(mx+1-ibc,j,m)
  235       continue
c     # negate the normal velocity:
      do 236 ibc=1,mbc
         do 236 j = 1-mbc, my+mbc
            q(mx+ibc,j,2) = -q(mx+1-ibc,j,2)
  236    continue
      go to 299

  299 continue
c
c-------------------------------------------------------
c     # bottom boundary:
c-------------------------------------------------------
      go to (300,310,320,330) mthbc(3)+1
      goto 399
c
  300 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(3)=0 and no BCs specified in bc2'
      stop
      go to 399
c
  310 continue
c     # zero-order extrapolation:
      do 315 m=1,meqn
         do 315 jbc=1,mbc
            do 315 i = 1-mbc,mx+mbc
               q(i,1-jbc,m) = q(i,1,m)
  315       continue
      go to 399

  320 continue
c     # periodic:
      do 325 m=1,meqn
         do 325 jbc=1,mbc
            do 325 i = 1-mbc, mx+mbc
               q(i,1-jbc,m) = q(i,my+1-jbc,m)
  325       continue
      go to 399

  330 continue
c     # solid wall (assumes 3'rd component is velocity or momentum in y):
      do 335 m=1,meqn
         do 335 jbc=1,mbc
            do 335 i = 1-mbc,mx+mbc
               q(i,1-jbc,m) = q(i,jbc,m)
  335       continue
c     # negate the normal velocity:
      do 336 jbc=1,mbc
         do 336 i = 1-mbc, mx+mbc
            q(i,1-jbc,3) = -q(i,jbc,3)
  336    continue
      go to 399

  399 continue
c
c-------------------------------------------------------
c     # top boundary:
c-------------------------------------------------------
      go to (400,410,420,430) mthbc(4)+1
      goto 499
c
  400 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(4)=0 and no BCs specified in bc2'
      stop
      go to 499

  410 continue
c     # zero-order extrapolation:
      do 415 m=1,meqn
         do 415 jbc=1,mbc
            do 415 i = 1-mbc,mx+mbc
               q(i,my+jbc,m) = q(i,my,m)
  415       continue
      go to 499

  420 continue
c     # periodic:
      do 425 m=1,meqn
         do 425 jbc=1,mbc
            do 425 i = 1-mbc, mx+mbc
               q(i,my+jbc,m) = q(i,jbc,m)
  425       continue
      go to 499

  430 continue
c     # solid wall (assumes 3'rd component is velocity or momentum in y):
      do 435 m=1,meqn
         do 435 jbc=1,mbc
            do 435 i = 1-mbc, mx+mbc
               q(i,my+jbc,m) = q(i,my+1-jbc,m)
  435       continue
c     # negate the normal velocity:
      do 436 jbc=1,mbc
         do 436 i = 1-mbc, mx+mbc
            q(i,my+jbc,3) = -q(i,my+1-jbc,3)
  436    continue
      go to 499

  499 continue

      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
MODULE ArrMod

  USE MathConstants

  ! Variables for arrival information
  IMPLICIT NONE
  INTEGER, PARAMETER, PRIVATE :: ARRFile = 36
  INTEGER                     :: MaxNArr
  INTEGER, ALLOCATABLE        :: NArr( :, : ), NArr3D( :, :, : )

  TYPE Arrival
     INTEGER :: NTopBnc, NBotBnc
     REAL    :: SrcDeclAngle, SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle, A, Phase
     COMPLEX :: delay
  END TYPE

  TYPE(Arrival), ALLOCATABLE :: Arr( :, :, : ), Arr3D( :, :, :, : )

CONTAINS

  SUBROUTINE AddArr( omega, id, ir, Amp, Phase, delay, SrcDeclAngle, RcvrDeclAngle, NumTopBnc, NumBotBnc )

    ! ADDs the amplitude and delay for an ARRival into a matrix of same.
    ! Extra logic included to keep only the strongest arrivals.

    REAL,      PARAMETER :: PhaseTol = 0.05  ! arrivals with essentially the same phase are grouped into one
    INTEGER,              INTENT( IN ) :: NumTopBnc, NumBotBnc, id, ir
    REAL    ( KIND = 8 ), INTENT( IN ) :: omega, Amp, Phase, SrcDeclAngle, RcvrDeclAngle
    COMPLEX ( KIND = 8 ), INTENT( IN ) :: delay
    LOGICAL              :: NewRay
    INTEGER              :: iArr( 1 ), Nt
    REAL                 :: AmpTot, w1, w2

    Nt     = NArr( id, ir )    ! # of arrivals
    NewRay = .TRUE.

    ! Is this the second bracketting ray of a pair?
    ! If so, we want to combine the arrivals to conserve space.
    ! (test this by seeing if the arrival time is close to the previous one)
    ! (also need that the phase is about the same to make sure surface and direct paths are not joined)

    IF ( Nt >= 1 ) THEN
       IF ( omega * ABS( delay - Arr( id, ir, Nt )%delay ) < PhaseTol .AND. &
           ABS( Arr( id, ir, Nt )%phase - Phase )       < PhaseTol ) NewRay = .FALSE.
    END IF

    IF ( NewRay ) THEN
       IF ( Nt >= MaxNArr ) THEN       ! space available to add an arrival?
          iArr = MINLOC( Arr( id, ir, : )%A )                       ! no: replace weakest arrival
          IF ( Amp > Arr( id, ir, iArr( 1 ) )%A ) THEN
             Arr( id, ir, iArr( 1 ) )%A             = SNGL( Amp )       ! amplitude
             Arr( id, ir, iArr( 1 ) )%Phase         = SNGL( Phase )     ! phase
             Arr( id, ir, iArr( 1 ) )%delay         = CMPLX( delay )    ! delay time
             Arr( id, ir, iArr( 1 ) )%SrcDeclAngle  = SNGL( SrcDeclAngle )  ! angle
             Arr( id, ir, iArr( 1 ) )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle
             Arr( id, ir, iArr( 1 ) )%NTopBnc       = NumTopBnc         ! Number of top     bounces
             Arr( id, ir, iArr( 1 ) )%NBotBnc       = NumBotBnc         !   "       bottom
          ENDIF
       ELSE
          NArr( id, ir         )               = Nt + 1              ! # of arrivals
          Arr(  id, ir, Nt + 1 )%A             = SNGL( Amp )         ! amplitude
          Arr(  id, ir, Nt + 1 )%Phase         = SNGL( Phase )       ! phase
          Arr(  id, ir, Nt + 1 )%delay         = CMPLX( delay )      ! delay time
          Arr(  id, ir, Nt + 1 )%SrcDeclAngle  = SNGL( SrcDeclAngle )    ! angle
          Arr(  id, ir, Nt + 1 )%RcvrDeclAngle = SNGL( RcvrDeclAngle )   ! angle
          Arr(  id, ir, Nt + 1 )%NTopBnc       = NumTopBnc           ! Number of top     bounces
          Arr(  id, ir, Nt + 1 )%NBotBnc       = NumBotBnc           !   "       bottom
       ENDIF
    ELSE      ! not a new ray
       !PhaseArr(   id, ir, Nt ) = PhaseArr( id, ir, Nt )

       ! calculate weightings of old ray information vs. new, based on amplitude of the arrival
       AmpTot = Arr( id, ir, Nt )%A + SNGL( Amp )
       w1     = Arr( id, ir, Nt )%A / AmpTot
       w2     = REAL( Amp ) / AmpTot

       Arr( id, ir, Nt )%delay         = w1 * Arr( id, ir, Nt )%delay         + w2 * CMPLX( delay ) ! weighted sum
       Arr( id, ir, Nt )%A             = AmpTot
       Arr( id, ir, Nt )%SrcDeclAngle  = w1 * Arr( id, ir, Nt )%SrcDeclAngle  + w2 * SNGL( SrcDeclAngle  )
       Arr( id, ir, Nt )%RcvrDeclAngle = w1 * Arr( id, ir, Nt )%RcvrDeclAngle + w2 * SNGL( RcvrDeclAngle )
    ENDIF

    RETURN
  END SUBROUTINE AddArr

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

  SUBROUTINE WriteArrivalsASCII( r, Nrd, Nr, SourceType )

    ! Writes the arrival data (Amplitude, delay for each eigenray)
    ! ASCII output file

    INTEGER,           INTENT( IN ) :: Nrd, Nr
    REAL,              INTENT( IN ) :: r( Nr )
    CHARACTER (LEN=1), INTENT( IN ) :: SourceType
    INTEGER           :: ir, id, iArr
    REAL     (KIND=8) :: factor

    WRITE( ARRFile, * ) MAXVAL( NArr( 1 : Nrd, 1 : Nr ) )

    DO id = 1, Nrd
       DO ir = 1, Nr
          IF ( SourceType == 'X' ) THEN   ! line source
             factor =  4.0 * SQRT( pi )
          ELSE                            ! point source
             IF ( r ( ir ) == 0 ) THEN
                factor = 1e5                   ! avoid /0 at origin
             ELSE
                factor = 1. / SQRT( r( ir ) )  ! cyl. spreading
             END IF
          END IF

          WRITE( ARRFile, * ) NArr( id, ir )
          DO iArr = 1, NArr( id, ir )
             ! You can compress the output file a lot by putting in an explicit format statement here ...
             ! However, you'll need to make sure you keep adequate precision
             WRITE( ARRFile, * ) &
             SNGL( factor ) * Arr( id, ir, iArr )%A,             &
             SNGL( RadDeg ) * Arr( id, ir, iArr )%Phase,         &
                        REAL( Arr( id, ir, iArr )%delay ),       &
                       AIMAG( Arr( id, ir, iArr )%delay ),       &
                              Arr( id, ir, iArr )%SrcDeclAngle,  &
                              Arr( id, ir, iArr )%RcvrDeclAngle, &
                              Arr( id, ir, iArr )%NTopBnc,       &
                              Arr( id, ir, iArr )%NBotBnc
          END DO  ! next arrival
       END DO  ! next receiver depth
    END DO  ! next range

    RETURN
  END SUBROUTINE WriteArrivalsASCII

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

  SUBROUTINE WriteArrivalsBinary( r, Nrd, Nr, SourceType )

    ! Writes the arrival data (amplitude, delay for each eigenray)
    ! Binary output file

    INTEGER,           INTENT( IN ) :: Nrd, Nr
    REAL,              INTENT( IN ) :: r( Nr )
    CHARACTER (LEN=1), INTENT( IN ) :: SourceType
    INTEGER           :: ir, id, iArr
    REAL     (KIND=8) :: factor

    WRITE( ARRFile ) MAXVAL( NArr( 1 : Nrd, 1 : Nr ) )

    DO id = 1, Nrd
       DO ir = 1, Nr
          IF ( SourceType == 'X' ) THEN   ! line source
             factor = 4.0 * SQRT( pi )
          ELSE                            ! point source
             IF ( r ( ir ) == 0 ) THEN
                factor = 1e5                   ! avoid /0 at origin
             ELSE
                factor = 1. / SQRT( r( ir ) )  ! cyl. spreading
             END IF
          END IF

          WRITE( ARRFile ) NArr( id, ir )

          DO iArr = 1, NArr( id, ir )
             ! integers written out as reals below for fast reading in Matlab
             WRITE( ARRFile ) &
            SNGL( factor * Arr( id, ir, iArr )%A ),           &
            SNGL( RadDeg * Arr( id, ir, iArr )%Phase ),       &
                           Arr( id, ir, iArr )%delay,         &
                           Arr( id, ir, iArr )%SrcDeclAngle,  &
                           Arr( id, ir, iArr )%RcvrDeclAngle, &
                     REAL( Arr( id, ir, iArr )%NTopBnc ),     &
                     REAL( Arr( id, ir, iArr )%NBotBnc )

          END DO   ! next arrival
       END DO   ! next receiver depth
    END DO   ! next range

    RETURN
  END SUBROUTINE WriteArrivalsBinary

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

  SUBROUTINE AddArr3D( omega, itheta, id, ir, Amp, Phase, delay, SrcDeclAngle, &
       SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle, NumTopBnc, NumBotBnc )

    ! ADDs the amplitude and delay for an ARRival into a matrix of same.
    ! Extra logic included to keep only the strongest arrivals.

    REAL,                 PARAMETER    :: PhaseTol = 0.5  ! arrivals with essentially the same phase are grouped into one
    INTEGER,              INTENT( IN ) :: itheta, id, ir
    INTEGER,              INTENT( IN ) :: NumTopBnc, NumBotBnc
    REAL    ( KIND = 8 ), INTENT( IN ) :: omega, Amp, Phase, SrcDeclAngle, SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle

    COMPLEX ( KIND = 8 ), INTENT( IN ) :: delay
    LOGICAL                            :: NewRay
    INTEGER                            :: iArr( 1 ), Nt
    REAL                               :: AmpTot, w1, w2

    Nt     = NArr3D( itheta, id, ir )    ! # of arrivals
    NewRay = .TRUE.

    ! Is this the second bracketting ray of a pair?
    ! If so, we want to combine the arrivals to conserve space.
    ! (test this by seeing if the arrival time is close to the previous one)
    ! (also need that the phase is about the same to make sure surface and direct paths are not joined)

    IF ( Nt >= 1 ) THEN
       IF ( omega * ABS( delay - Arr3D( itheta,  id, ir, Nt )%delay ) < PhaseTol .AND. &
           ABS( Arr3D( itheta,  id, ir, Nt )%phase - Phase )       < PhaseTol ) NewRay = .FALSE.
    END IF

    IF ( NewRay ) THEN
       IF ( Nt >= MaxNArr ) THEN       ! space available to add an arrival?
          iArr = MINLOC( Arr3D( itheta,  id, ir, : )%A )                       ! no: replace weakest arrival
          IF ( Amp > Arr3D( itheta,  id, ir, iArr( 1 ) )%A ) THEN
             Arr3D( itheta,  id, ir, iArr( 1 ) )%A             = SNGL( Amp )       ! amplitude
             Arr3D( itheta,  id, ir, iArr( 1 ) )%Phase         = SNGL( Phase )     ! phase
             Arr3D( itheta,  id, ir, iArr( 1 ) )%delay         = CMPLX( delay )    ! delay time
             Arr3D( itheta,  id, ir, iArr( 1 ) )%SrcDeclAngle  = SNGL( SrcDeclAngle  ) ! angle
             Arr3D( itheta,  id, ir, iArr( 1 ) )%SrcAzimAngle  = SNGL( SrcAzimAngle  ) ! angle
             Arr3D( itheta,  id, ir, iArr( 1 ) )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle
             Arr3D( itheta,  id, ir, iArr( 1 ) )%RcvrAzimAngle = SNGL( RcvrAzimAngle ) ! angle
             Arr3D( itheta,  id, ir, iArr( 1 ) )%NTopBnc       = NumTopBnc         ! Number of top     bounces
             Arr3D( itheta,  id, ir, iArr( 1 ) )%NBotBnc       = NumBotBnc         !   "       bottom
          ENDIF
       ELSE
          NArr3D( itheta,  id, ir         )               = Nt + 1              ! # of arrivals
          Arr3D( itheta,   id, ir, Nt + 1 )%A             = SNGL( Amp )         ! amplitude
          Arr3D( itheta,   id, ir, Nt + 1 )%Phase         = SNGL( Phase )       ! phase
          Arr3D( itheta,   id, ir, Nt + 1 )%delay         = CMPLX( delay )      ! delay time
          Arr3D( itheta,   id, ir, Nt + 1 )%SrcDeclAngle  = SNGL( SrcDeclAngle  )   ! angle
          Arr3D( itheta,   id, ir, Nt + 1 )%SrcAzimAngle  = SNGL( SrcAzimAngle  )   ! angle
          Arr3D( itheta,   id, ir, Nt + 1 )%RcvrDeclAngle = SNGL( RcvrDeclAngle )   ! angle
          Arr3D( itheta,   id, ir, Nt + 1 )%RcvrAzimAngle = SNGL( RcvrAzimAngle )   ! angle
          Arr3D( itheta,   id, ir, Nt + 1 )%NTopBnc       = NumTopBnc           ! Number of top     bounces
          Arr3D( itheta,   id, ir, Nt + 1 )%NBotBnc       = NumBotBnc           !   "       bottom
       ENDIF
    ELSE      ! not a new ray
       !PhaseArr(   id, ir, Nt ) = PhaseArr( id, ir, Nt )
       ! calculate weightings of old ray information vs. new, based on amplitude of the arrival
       AmpTot = Arr3D( itheta, id, ir, Nt )%A + SNGL( Amp )
       w1     = Arr3D( itheta, id, ir, Nt )%A / AmpTot
       w2     = REAL( Amp ) / AmpTot

       Arr3D( itheta, id, ir, Nt )%delay     = w1 * Arr3D( itheta, id, ir, Nt )%delay     + w2 * CMPLX( delay ) ! weighted sum
       Arr3D( itheta, id, ir, Nt )%A         = AmpTot
       Arr3D( itheta, id, ir, Nt )%SrcDeclAngle  = w1 * Arr3D( itheta, id, ir, Nt )%SrcDeclAngle  + w2 * SNGL( SrcDeclAngle )
       Arr3D( itheta, id, ir, Nt )%SrcAzimAngle  = w1 * Arr3D( itheta, id, ir, Nt )%SrcAzimAngle  + w2 * SNGL( SrcAzimAngle )
       Arr3D( itheta, id, ir, Nt )%RcvrDeclAngle = w1 * Arr3D( itheta, id, ir, Nt )%RcvrDeclAngle + w2 * SNGL( RcvrDeclAngle )
       Arr3D( itheta, id, ir, Nt )%RcvrAzimAngle = w1 * Arr3D( itheta, id, ir, Nt )%RcvrAzimAngle + w2 * SNGL( RcvrAzimAngle )
    ENDIF

    RETURN
  END SUBROUTINE AddArr3D

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

  SUBROUTINE WriteArrivalsASCII3D( Ntheta, Nrd, Nr )

    ! Writes the arrival data (Amplitude, delay for each eigenray)
    ! ASCII output file

    INTEGER, INTENT( IN ) :: Ntheta, Nrd, Nr
    INTEGER               :: itheta, ir, id, iArr

    WRITE( ARRFile, * ) MAXVAL( NArr3D( 1 : Ntheta,  1 : Nrd, 1 : Nr ) )

    DO itheta = 1, Ntheta
       DO id = 1, Nrd
          DO ir = 1, Nr

             WRITE( ARRFile, * ) NArr3D( itheta,  id, ir )

             DO iArr = 1, NArr3D( itheta,  id, ir )
                ! you can compress the output file a lot by putting in an explicit format statement here ...
                ! However, you'll need to make sure you keep adequate precision
                WRITE( ARRFile, * ) &
                              Arr3D( itheta, id, ir, iArr )%A,             &
                     RadDeg * Arr3D( itheta, id, ir, iArr )%Phase,         &
                        REAL( Arr3D( itheta, id, ir, iArr )%delay ),       &
                       AIMAG( Arr3D( itheta, id, ir, iArr )%delay ),       &
                              Arr3D( itheta, id, ir, iArr )%SrcDeclAngle,  &
                              Arr3D( itheta, id, ir, iArr )%SrcAzimAngle,  &
                              Arr3D( itheta, id, ir, iArr )%RcvrDeclAngle, &
                              Arr3D( itheta, id, ir, iArr )%RcvrAzimAngle, &
                              Arr3D( itheta, id, ir, iArr )%NTopBnc,       &
                              Arr3D( itheta, id, ir, iArr )%NBotBnc
             END DO  ! next arrival
          END DO  ! next receiver depth
       END DO  ! next range
    END DO   ! next angle

    RETURN
  END SUBROUTINE WriteArrivalsASCII3D

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

  SUBROUTINE WriteArrivalsBinary3D( Ntheta, Nrd, Nr )

    ! Writes the arrival data (amplitude, delay for each eigenray)
    ! Binary output file

    INTEGER, INTENT( IN ) :: Ntheta, Nrd, Nr
    INTEGER               :: itheta, ir, id, iArr

    WRITE( ARRFile ) MAXVAL( NArr3D( 1 : Ntheta,  1 : Nrd, 1 : Nr ) )

    DO itheta = 1, Ntheta
       DO id = 1, Nrd
          DO ir = 1, Nr

             WRITE( ARRFile ) NArr3D( itheta,  id, ir )

             DO iArr = 1, NArr3D( itheta,  id, ir )
                ! integers written out as reals below for fast reading in Matlab
                WRITE( ARRFile ) &
                              Arr3D( itheta, id, ir, iArr )%A,             &
               SNGL( RadDeg * Arr3D( itheta, id, ir, iArr )%Phase ),       &
                              Arr3D( itheta, id, ir, iArr )%delay,         &
                              Arr3D( itheta, id, ir, iArr )%SrcDeclAngle,  &
                              Arr3D( itheta, id, ir, iArr )%SrcAzimAngle,  &
                              Arr3D( itheta, id, ir, iArr )%RcvrDeclAngle, &
                              Arr3D( itheta, id, ir, iArr )%RcvrAzimAngle, &
                        REAL( Arr3D( itheta, id, ir, iArr )%NTopBnc ),     &
                        REAL( Arr3D( itheta, id, ir, iArr )%NBotBnc )
             END DO   ! next arrival
          END DO   ! next receiver depth
       END DO   ! next range
    END DO   ! next angle

    RETURN
  END SUBROUTINE WriteArrivalsBinary3D

END MODULE ArrMod

max_stars_repo_name: IvanaEscobar/sandbox

max_issues_repo_name: IvanaEscobar/sandbox

max_forks_repo_name: IvanaEscobar/sandbox


text: Ali Ramezanzadeh is the Campus Outreach Director for the ASUCD Office of External Affairs.