text:
module qsort_struc_m
use structure_c
implicit none
private
public :: qsort

contains

recursive subroutine qsort(struc_rep)
  type(structure), dimension(:) :: struc_rep
  integer :: iq

  if(size(struc_rep) > 1) then
     call Partition(struc_rep,iq)
     call qsort(struc_rep(:iq-1))
     call qsort(struc_rep(iq:))
  end if
end subroutine qsort

subroutine Partition(struc_rep,marker)
  type(structure), dimension(:) :: struc_rep
  type(structure) :: x, temp
  integer, intent(out) :: marker
  integer :: i, j

  x = struc_rep(1)
  i= 0
  j= size(struc_rep) + 1

  do
     j = j-1
     do
        if (struc_rep(j)%energy <= x%energy) exit
        j = j-1
     end do
     i = i+1
     do
        if (struc_rep(i)%energy >= x%energy) exit
        i = i+1
     end do
     if (i < j) then
        ! exchange A(i) and A(j)
        temp = struc_rep(i)
        struc_rep(i) = struc_rep(j)
        struc_rep(j) = temp
     else if (i == j) then
        marker = i+1
        return
     else
        marker = i
        return
     endif
  end do

end subroutine Partition

end module

max_stars_repo_name: frostbaum/RMSDysator

max_issues_repo_name: frostbaum/RMSDysator

max_forks_repo_name: frostbaum/RMSDysator


text:
Unfortunately, Davis has a fair number of methamphetamine addicts, commonly known as tweakers. They are one of the reasons it is always wise to lock up your Bicycling bike.

The term tweaker technically refers to a stage within the methamphetamine abuse cycle.  From the State of California Resource Center, publication No. (ADP) 992455 (http://www.stopdrugs.org/tweaking.html link):

  The most dangerous stage of meth abuse is called tweaking.  A tweaker is an abuser who probably has not slept in 315 days and is irritable and paranoid.  Tweakers often behave or react violently and if a tweaker is using alcohol or another depressant, negative feelings and associated dangers intensify.  The tweaker craves more meth, but no dosage will help recreate the euphoric high, causing frustration, unpredictability, and the potential for Tragic Events violence.

  A tweaker can appear normal: eyes can be clear, speech concise, and movements brisk.  But a closer look will reveal the persons eyes are moving ten times faster than normal, the voice has a slight quiver, and movements are quick and jerky.  Tweakers are often involved in domestic disputes and motor vehicle accidents.  They may hang out at raves and parties and participate in spurofthemoment Crime crimes, such as purse snatchings and assaults, to support their habit.

  Detaining a tweaker alone is not recommended and law enforcement officers should call for backup.

Do we have any stats on this?  I didnt think that meth was a big problem in davis?  Users/GeorgeLewis George Lewis

I spent some time sitting in my fair share of court rooms in Yolo County and lets just say that I heard the word methampetamine come up a lot.  Users/RobRoy

Yes, its hardly unheard of.  Users/MichaelGiardina

Oh, what demons these tweakers are!  In fact, sounds like an undercover cop  Users/ChristopherMckenzie

Some statistics: the Davis Police logs http://cityofdavis.org/police/log/press/ from 04/0109/27/2005 list three meth incidents. However, all three were reported during the past few days (09/2409/26), which seems unlikely. Maybe they drop reports from the online logs? RAND corporation seems to have very detailed stats about drug abuse http://ca.rand.org/stats/health/drug.html  does anyone have access? Users/JoFeuerstein

That log only has data for reports filed in the past week. Users/RyanCastellucci

See Methamphetamine for a list of Meth House busts in California including Davis. They are sorted by county.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c
c     --------------------------------
      character function d2ch ( idig )
c     --------------------------------
c
c     Converts a digit 0 <= idig <= 9 into a corresponding character
c
      implicit none
c
      integer idig
c
      if ( (idig .ge. 0).and.(idig.le.9) ) then
        d2ch = char(48 + idig)
      else
        write(*,*) 'd2ch error : idig =', idig,' is not a digit'
      endif
c
      return
      end
c
c     ------------------------------------------------------------
      subroutine SetViewAngle ( dPhi, dTheta, CTH, STH, CPH, SPH )
c     ------------------------------------------------------------
c
c     Sets sines and cosines of the viewing angle
c
c     input  : dPhi, dTheta - viewing angles in degrees
c     output : CPH, SPH, CTH, STH - cos(Phi),sin(Phi),cos(Theta),sin(Theta)
c
      implicit none
c
      real pi, d2r
      parameter ( pi  = 3.14159265358978245 )
      parameter ( d2r = pi / 180.0 )
c
      real dPhi, dTheta, Phi, Theta, CPH, SPH, CTH, STH
c
      Phi   = dPhi * d2r
      CPH   = cos(Phi)
      SPH   = sin(Phi)
      Theta = dTheta * d2r
      CTH   = cos(Theta)
      STH   = sin(Theta)
c
      return
      end
c
c     ---------------------------------
      subroutine InitGIFDevice ( name )
c     ---------------------------------
c
c     Open GIF file with name 'name'
c
      implicit none
c
      real ScreenSize
      character*15 name
c
      call pgbegin ( 0 , name , 1 , 1 )
      call pgpap ( 7.515 , 0.75 )  ! 640x480
c      call pgpap ( 3.0 , 1.0 )  ! 256x256
      call pgask ( .false. )
      call pgenv ( -0.5 , 0.5 , -0.5 , 0.5 , 1 , -2 )
      call pgslw ( 1 ) ! set minimum possible linewidth
c
      return
      end
c
c     ------------------------------------------
      integer function InitDevice ( name , nch )
c     ------------------------------------------
c
c     Open an arbitrary PGPLOT device
c
c     input : name - the name of PGPLOT device
c             nch  - number of characters in the Device name
c                    (length of the string name)
c
      implicit none
c
      integer nch
      character*50 name
c
      real w,sd
      integer pgopen, istat
      external pgopen
c
      istat = pgopen ( name(1:nch) )
      if ( istat .le. 0 ) then
        write(*,*) 'error: InitDevice: cannot open device',name
        stop
      else
        InitDevice = istat
      endif
c
      w  = 10.0
      sd = 1.0
c
      call PGPAP ( w , sd )
      call PGASK ( .false. )
      call PGENV ( -0.5 , 0.5 , -0.5 , 0.5 , 1 , -2 )
      call PGSLW ( 1 ) ! set minimum possible linewidth
      call PGSCF ( 2 )
c
      return
      end
c

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

      subroutine DrawFrame ( F , psu )
c     --------------------------------
c
c     Draws a wireframe
c
c     input : F(8,3)  - 3D coordinates of the 8 frame corners
c             psu(11) - plot setup
c
      implicit none
c
      real F(8,3), psu(11)
c
      real Fp(8,2) ! array of coordinates of projected corners
c
      integer i
      do i = 1 , 8
        call OnScreen ( F(i,1) , F(i,2) , F(i,3) , psu ,
     &                  Fp(i,1) , Fp(i,2)
     &                 )
      enddo
c
      call PGSCI ( 1 )  ! draw the frame with a color negative to the BG
c
      call Line ( Fp(1,1) , Fp(1,2) , Fp(2,1) , Fp(2,2) )
      call Line ( Fp(1,1) , Fp(1,2) , Fp(4,1) , Fp(4,2) )
      call Line ( Fp(2,1) , Fp(2,2) , Fp(3,1) , Fp(3,2) )
      call Line ( Fp(3,1) , Fp(3,2) , Fp(4,1) , Fp(4,2) )
      call Line ( Fp(5,1) , Fp(5,2) , Fp(6,1) , Fp(6,2) )
      call Line ( Fp(5,1) , Fp(5,2) , Fp(8,1) , Fp(8,2) )
      call Line ( Fp(6,1) , Fp(6,2) , Fp(7,1) , Fp(7,2) )
      call Line ( Fp(7,1) , Fp(7,2) , Fp(8,1) , Fp(8,2) )
      call Line ( Fp(1,1) , Fp(1,2) , Fp(5,1) , Fp(5,2) )
      call Line ( Fp(2,1) , Fp(2,2) , Fp(6,1) , Fp(6,2) )
      call Line ( Fp(3,1) , Fp(3,2) , Fp(7,1) , Fp(7,2) )
      call Line ( Fp(4,1) , Fp(4,2) , Fp(8,1) , Fp(8,2) )
c
      return
      end

c     -------------------------------------
      subroutine Line ( X1 , Y1 , X2 , Y2 )
c     -------------------------------------
c
c     Draws a line between points (X1,Y1) and (X2,Y2)
c     with clipping disabled
c
      implicit none
c
      real X1 , X2 , Y1 , Y2
c
      real  Xl(2) , Yl(2)
c
      call PGSCLP ( 0 )    ! disable clipping at the viewport edges
c
      Xl(1) = X1
      Yl(1) = Y1
      Xl(2) = X2
      Yl(2) = Y2
c
      call PGLine ( 2 , Xl , Yl )
c
      return
      end
c
c     -------------------------------------------------
      subroutine OnScreen ( X , Y , Z , psu , sX , sY )
c     -------------------------------------------------
c
c     From 3D coordinates of a point calculates a 2D screen coordinates
c
c     input  : X , Y , Z  - 3D coordinates of a point
c                           ( in the range [-0.5+Xb,0.5+Xb] by
c                             convention, where Xb is a coordinate
c                             of the box center )
c              psu(1)  = Xb    ! x-coordinate (3D) of the box center
c              psu(2)  = Yb    ! y-coordinate (3D) of the box center
c              psu(3)  = Zb    ! z-coordinate (3D) of the box center
c              psu(4)  = Scale ! scale of the plot for zoom in/out
c              psu(5)  = CTH   ! cosine of the viewing angle Theta
c              psu(6)  = STH   ! sine -"-
c              psu(7)  = CPH   ! cosine of the viewing angle Phi
c              psu(8)  = SPH   ! sine -"-
c              psu(9)  = Xlook ! x-coordinate (3D) of the viewing point
c              psu(10) = Ylook ! y-coordinate (3D) of the viewing point
c              psu(11) = Zlook ! z-coordinate (3D) of the viewing point
c
c     output : sX , sY - 2D screen coordinates
c
      implicit none
c
      real X, Y, Z, psu(11)
      real sX , sY
c
      real Xd, Yd, Zd
c
      Xd = X*psu(5)*psu(7) + Y*psu(5)*psu(8) - Z*psu(6) - psu(1)
      Yd = Y*psu(7) - X*psu(8) - psu(2)
      Zd = X*psu(6)*psu(7) + Y*psu(6)*psu(8) + Z*psu(5) - psu(3)
      sX = psu(4) * psu(11) * Xd / abs(psu(11)-Zd) + psu(9)
      sY = psu(4) * psu(11) * Yd / abs(psu(11)-Zd) + psu(10)
c
      return
      end
c
c     --------------------------
      subroutine InitFrame ( F )
c     --------------------------
c
c     Sets array F(8,3) with coordinates of the box corners
c
      implicit none
c
      real F(8,3)
c
      F(1,1) = -0.5
      F(1,2) = -0.5
      F(1,3) = -0.5
      F(2,1) = -0.5
      F(2,2) =  0.5
      F(2,3) = -0.5
      F(3,1) =  0.5
      F(3,2) =  0.5
      F(3,3) = -0.5
      F(4,1) =  0.5
      F(4,2) = -0.5
      F(4,3) = -0.5
      F(5,1) = -0.5
      F(5,2) = -0.5
      F(5,3) =  0.5
      F(6,1) = -0.5
      F(6,2) =  0.5
      F(6,3) =  0.5
      F(7,1) =  0.5
      F(7,2) =  0.5
      F(7,3) =  0.5
      F(8,1) =  0.5
      F(8,2) = -0.5
      F(8,3) =  0.5
c
      return
      end
c
c     --------------------------------------------
      SUBROUTINE PALETT ( TYPE , CONTRA , BRIGHT )
c     --------------------------------------------
c
c     Sets a "palette" of colors in the range of defined color indices
c     This subroutine is distributed with PGPLOT in one of the demos
c
      INTEGER TYPE
      REAL CONTRA, BRIGHT
C
      REAL GL(2), GR(2), GG(2), GB(2)
      REAL RL(9), RR(9), RG(9), RB(9)
      REAL HL(5), HR(5), HG(5), HB(5)
      REAL CL(5), CR(5), CG(5), CB(5)
      REAL WL(10), WR(10), WG(10), WB(10)
      REAL AL(20), AR(20), AG(20), AB(20)
C
      DATA GL /0.0, 1.0/
      DATA GR /0.0, 1.0/
      DATA GG /0.0, 1.0/
      DATA GB /0.0, 1.0/
C
      DATA RL /-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7/
      DATA RR / 0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0/
      DATA RG / 0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0/
      DATA RB / 0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0/
C
      DATA HL /0.0, 0.2, 0.4, 0.6, 1.0/
      DATA HR /0.0, 0.5, 1.0, 1.0, 1.0/
      DATA HG /0.0, 0.0, 0.5, 1.0, 1.0/
      DATA HB /0.0, 0.0, 0.0, 0.3, 1.0/
C
      DATA CL /0.0, 0.2, 0.4, 0.6, 1.0/
      DATA CR /0.0, 0.0, 0.0, 0.3, 1.0/
      DATA CG /0.0, 0.0, 0.5, 1.0, 1.0/
      DATA CB /0.0, 0.5, 1.0, 1.0, 1.0/
C
      DATA WL /0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0/
      DATA WR /0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0/
      DATA WG /0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0/
      DATA WB /0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0/
C
      DATA AL /0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5,
     :         0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0/
      DATA AR /0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0,
     :         0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
      DATA AG /0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8,
     :         0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0/
      DATA AB /0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9,
     :         0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
C
      IF (TYPE.EQ.1) THEN
C        -- gray scale
         CALL PGCTAB(GL, GR, GG, GB, 2, CONTRA, BRIGHT)
      ELSE IF (TYPE.EQ.2) THEN
C        -- rainbow
         CALL PGCTAB(RL, RR, RG, RB, 9, CONTRA, BRIGHT)
      ELSE IF (TYPE.EQ.3) THEN
C        -- heat
         CALL PGCTAB(HL, HR, HG, HB, 5, CONTRA, BRIGHT)
      ELSE IF (TYPE.EQ.4) THEN
C        -- freeze
         CALL PGCTAB(CL, CR, CG, CB, 5, CONTRA, BRIGHT)
      ELSE IF (TYPE.EQ.5) THEN
C        -- AIPS
         CALL PGCTAB(AL, AR, AG, AB, 20, CONTRA, BRIGHT)
      END IF
      END
c
c     ------------------------------------------------
      SUBROUTINE FIDDLE ( P , SIGN , CONTRA , BRIGHT )
c     ------------------------------------------------

      end
      subroutine SetViewAngle2 ( dPhi, dTheta, CTH, STH, CPH, SPH )
c     ------------------------------------------------------------
c
c     Sets sines and cosines of the viewing angle
c
c     input  : dPhi, dTheta - viewing angles in degrees
c     output : CPH, SPH, CTH, STH - cos(Phi),sin(Phi),cos(Theta),sin(Theta)
c
      implicit none
c
      real pi, d2r
      parameter ( pi  = 3.14159265358978245 )
      parameter ( d2r = pi / 180.0 )
c
      real dPhi, dTheta, Phi, Theta, CPH, SPH, CTH, STH
c
      Phi   = dPhi * d2r
      CPH   = cos(Phi)
      SPH   = sin(Phi)
      Theta = dTheta * d2r
      CTH   = cos(Theta)
      STH   = sin(Theta)
c
      return
      end

max_stars_repo_name: leginee/netbeans

max_issues_repo_name: leginee/netbeans

max_forks_repo_name: leginee/netbeans


text:
Davis Beer Week
http://www.davisbeerweek.com
Davis Beer Week founded by de Vere’s Irish Pub, strives to highlight and promote our community’s thriving craft beer culture. This 7day citywide festival will attract beer tourism, foster knowledge of our regional brewing heritage, and serve as a showcase for great beers, restaurants, pubs, and other businesses with ties to the craft beer community.

Enjoy an entire week of incredible beer events throughout the community, which will be capped off with our secondannual Bike & Brew Fest, an evening of entertainment, food, wine and beer benefiting the Davis United States Bicycling Hall of Fame!

Bike & Brew Fest Beer Vendors
American Craft

 101 North
 21st Amendment
 Anchor
 Angry Orchard
 Anderson Valley
 American River Brewing
 Berryessa Brewery
 Blue Moon
 Coronado Brewing
 Crispin Cider
 Device Brewing Company
 Drakes
 Eel River
 Firestone
 Fox Barrel Cider
 Grand Teton Brewing
 Karl Strauss
 Lagunitas
 Lost Coast
 Mad River Brewing Co.
 Mavericks Beers Half Moon Bay Brewing Co.
 Mendocino Brewing Co.
 North Coast
 Nectar Ales
 New Belgium
 Ninkasi Brewing Company
 Rubicon
 Ruhstaller
 Sam Adams
 Speakeasy
 Sierra Nevada
 Shipyard Brewing
 Six Rivers
 Stone
 Sudwerk
 Track 7
 Trumer Pils
 Old Glory Brewing Company
 Hangar 24 Craft Brewery
 Hoppy Brewing Company
 New Helvetia Brewing

International

 Affligem
 Weihenstephaner
 Stella
 Guinness
 Maes Pils
 Hoegaarden
 Leffe
 Morreti
 Erdinger
 Carlsberg
 Kronenbourg
 Crabbies

Wine

 Artesa
 Michael David
 Helwig
 Lange Twins
 Bogle
 Trinchero Family vineyards
 Da Luca
 Joseph Carr
 Rodney Strong

Food Vendors

 Bistro 33
 de Vere’s Irish Pub
 Dos Coyotes Border Cafe
 Jack’s Urban Eats
 Our House
 Paesanos
 Uncle Vitos
 Sudwerk
 The Hotdogger
 Monticello Seasonal Cuisine
 Thai Recipes

Food Trucks

 Wicked Wich
 California Love Truck
 The Big Red Bus
 Bacon Mania

Entertainment

 Tha Dirt Feeling
 Elements of Soul
 The Silent Comedy

Past Dates


module BoxLibRNGs
   use RNGEngine ! Interfaces for uniform/integer RNG generators written in C and part of BoxLib
   use Random_Numbers ! Used to generate unpredictable seeds or initial seeds for Marsenne twister
   use NonUniformRNGs ! Generates samples from several nonuniform distributions (we use our own for Guassian)
   use iso_c_binding
   implicit none
   public ! Export all names of public entities so this is the one unified Fortran RNG interface for the outside world

   integer, parameter, private :: dp=kind(0.0d0), sp=kind(0.0)

   ! Returns pseudorandom number in interval [0,1).
   ! Donev: The single precision generator should be tested at least once!
   interface UniformRNG
      module procedure UniformRNG_dp
      module procedure UniformRNG_sp
   end interface UniformRNG

   interface NormalRNG  ! The Fortran version of this is below and may be faster
      module procedure NormalRNG_dp
      module procedure NormalRNG_sp
   end interface

   interface UniformRNGVec  ! The Fortran version of this is below and may be faster
      module procedure UniformRNGs_sp
      module procedure UniformRNGs
   end interface

   interface NormalRNGVec  ! The Fortran version of this is below and may be faster
      module procedure NormalRNGs_sp
      module procedure NormalRNGs
   end interface

   interface PoissonRNG ! This is only scalar for now
      module procedure PoissonRNG_sp
      module procedure PoissonRNG_dp
   end interface

   interface BinomialRNG ! This is only scalar for now
      module procedure BinomialRNG_sp
      module procedure BinomialRNG_dp
   end interface

contains ! It is likely that vectorized versions will do better here

!-------------------------------------------
! Uniform RNGs -- interfaces to C Marsenne-Twister codes
!-------------------------------------------

   ! This is in principle callable by C directly, but we go through the wrapper here for safety
   ! void SeedRNG(int *seed); // On output, the actual seed used
   !
   subroutine SeedRNG(seed) BIND(C,NAME="SeedRNG")
      integer, intent(inout) :: seed ! If zero, the clock will be used to find an unpredictable seed

      if(seed==0) then
         call UnpredictableSeeds(seed)
         write(*,*) "SeedRNG @ BoxLibRNGs: Generated unpredictable SEED=", seed
      end if

      call SeedRNG_C(seed)

   end subroutine

   subroutine UniformRNG_dp(number, engine)
      ! Returns pseudorandom number in interval [0,1).
      real(dp), intent(out) :: number
      type(c_ptr), intent(in), optional :: engine
      if (present(engine)) then
         call UniformRNG_C(number, engine)
      else
         call UniformRNG_C(number)
      end if
    end subroutine

   subroutine UniformRNG_sp(number, engine)
      ! Returns pseudorandom number in interval [0,1).
      real(sp), intent(out) :: number
      type(c_ptr), intent(in), optional :: engine
      if (present(engine)) then
         call UniformRNG_C(number, engine)
      else
         call UniformRNG_C(number)
      end if
   end subroutine

  subroutine UniformRNGs(numbers, n_numbers, engine)
    integer, intent(in) :: n_numbers
    real(dp), intent(out) :: numbers(n_numbers)
    type(c_ptr), intent(in), optional :: engine

    integer :: i

    do i=1, n_numbers
       call UniformRNG(numbers(i), engine) ! Marsenne-Twister in C
    end do

  end subroutine

  subroutine UniformRNGs_sp(numbers, n_numbers, engine)
    integer, intent(in) :: n_numbers
    real(sp), intent(out) :: numbers(n_numbers)
    type(c_ptr), intent(in), optional :: engine

    integer :: i

    do i=1, n_numbers
       call UniformRNG(numbers(i), engine) ! Marsenne-Twister in C
    end do

  end subroutine

!-------------------------------------------
! Gaussian/normal RNGs
!-------------------------------------------

   ! Donev: Replaced this with Fortran implementation below instead of using NormalRNG_C
   subroutine NormalRNG_dp_C(number, engine)
      ! Returns normal N(0,1) variate
      real(dp), intent(out) :: number
      type(c_ptr), intent(in), optional :: engine
      if (present(engine)) then
         call NormalRNG_C(number, engine)
      else
         call NormalRNG_C(number)
      end if
   end subroutine

   subroutine NormalRNG_sp(number, engine)
      ! Returns normal N(0,1) variate
      real(sp), intent(out) :: number
      type(c_ptr), intent(in), optional :: engine
      real(dp) :: number_dp
      call NormalRNG_dp(number_dp, engine)
      number=number_dp
   end subroutine

  subroutine NormalRNGs(numbers, n_numbers, engine)
    integer, intent(in) :: n_numbers
    real(dp), intent(out) :: numbers(n_numbers)
    type(c_ptr), intent(in), optional :: engine

    integer :: i

    do i=1, n_numbers
       call NormalRNG(numbers(i), engine)
    end do

  end subroutine

  subroutine NormalRNGs_sp(numbers, n_numbers, engine)
    integer, intent(in) :: n_numbers
    real(sp), intent(out) :: numbers(n_numbers)
    type(c_ptr), intent(in), optional :: engine

    integer :: i

    do i=1, n_numbers
       call NormalRNG(numbers(i), engine)
    end do

  end subroutine

  ! Donev: This should be should be tested for histogram at least once or checked that given the same random integer
  !   it produces the same result as the C code
  subroutine NormalRNG_dp(number, engine)
      ! This is the Fortran equivalent of the C genrandn, it is a bit faster and better for inlining
      ! It also is faster than the normal generator in NURNGs.f90 since it avoids rejection, so we use it here
      ! Note that it is only accurate to 1E-9 (not an exact sampler!) but this is more than sufficient for applications
      real(dp), intent(inout) :: number
      type(c_ptr), intent(in), optional :: engine

      real(dp)     :: p

      real(dp) ::  q,r
      real(dp), parameter :: a1=-39.6968302866538d0
      real(dp), parameter :: a2=220.946098424521d0
      real(dp), parameter :: a3=-275.928510446969d0
      real(dp), parameter :: a4=138.357751867269d0
      real(dp), parameter :: a5=-30.6647980661472d0
      real(dp), parameter :: a6=2.50662827745924d0
      real(dp), parameter :: b1=-54.4760987982241d0
      real(dp), parameter :: b2=161.585836858041d0
      real(dp), parameter :: b3=-155.698979859887d0
      real(dp), parameter :: b4=66.8013118877197d0
      real(dp), parameter :: b5=-13.2806815528857d0
      real(dp), parameter :: c1=-0.00778489400243029d0
      real(dp), parameter :: c2=-0.322396458041136d0
      real(dp), parameter :: c3=-2.40075827716184d0
      real(dp), parameter :: c4=-2.54973253934373d0
      real(dp), parameter :: c5=4.37466414146497d0
      real(dp), parameter :: c6=2.93816398269878d0
      real(dp), parameter :: d1=0.00778469570904146d0
      real(dp), parameter :: d2=0.32246712907004d0
      real(dp), parameter :: d3=2.445134137143d0
      real(dp), parameter :: d4=3.75440866190742d0
      real(dp), parameter :: p_low=0.02425d0
      real(dp), parameter :: p_high=0.9575d0

      do
         call UniformRNG(p, engine)
         if (p .ne. 0.d0) exit
      end do

      if (p.lt.p_low) then
         q=dsqrt(-2.d0*dlog(p))
         number = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/   &
             ((((d1*q+d2)*q+d3)*q+d4)*q+1.d0)
      elseif (p.le.p_high)then
         q=p-0.5d0
         r=q*q
         number = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q/  &
                   (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.d0)
      else
         q=dsqrt(-2.d0*dlog(1.d0-p))
         number = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/   &
               ((((d1*q+d2)*q+d3)*q+d4)*q+1.d0)
      endif

  end subroutine

  subroutine NormalRNGFast(p, engine) ! This has the correct first and second moment only
    ! It is not an actual normal random generator!
    type(c_ptr), intent(in), optional :: engine
    real(dp)     :: u,p
    real(dp), parameter :: f = 3.46410161514d0

    call UniformRNG(u, engine)
    p = f*(u-0.5_dp)

  end subroutine

!-------------------------------------------
! Poisson/binomial RNGs
!-------------------------------------------

 SUBROUTINE PoissonRNG_dp(number,mean,engine)
    INTEGER, INTENT(OUT) :: number
    REAL(dp), INTENT(IN) :: mean
    type(c_ptr), intent(in), optional :: engine

    number=random_Poisson(mu=real(mean), engine=engine)

 END SUBROUTINE

 SUBROUTINE PoissonRNG_sp(number,mean,engine)
    INTEGER, INTENT(OUT) :: number
    REAL(sp), INTENT(IN) :: mean
    type(c_ptr), intent(in), optional :: engine

    number=random_Poisson(mu=mean, engine=engine)

 END SUBROUTINE

 SUBROUTINE BinomialRNG_dp(number,n_trials,success_prob, engine)
    INTEGER, INTENT(OUT) :: number
    INTEGER, INTENT(IN) :: n_trials ! Number of trials
    REAL(dp), INTENT(IN) :: success_prob ! Probability of successful trial
    type(c_ptr), intent(in), optional :: engine

    number=random_binomial(n=n_trials, pp=real(success_prob), engine=engine)

 END SUBROUTINE

 SUBROUTINE BinomialRNG_sp(number,n_trials,success_prob, engine)
    INTEGER, INTENT(OUT) :: number
    INTEGER, INTENT(IN) :: n_trials ! Number of trials
    REAL(sp), INTENT(IN) :: success_prob ! Probability of successful trial
    type(c_ptr), intent(in), optional :: engine

    number=random_binomial(n=n_trials, pp=success_prob, engine=engine)

 END SUBROUTINE

 ! This samples from a multinomial distribution
 ! The last sample is not sampled explicitly since it is just N-sum(samples)
 subroutine MultinomialRNG(samples, n_samples, N, p, engine)
    integer, intent(in) :: n_samples, N
    integer, intent(out) :: samples(n_samples)
    real(dp), intent(in) :: p(n_samples)
    type(c_ptr), intent(in), optional :: engine

    real(dp) :: sum_p
    integer :: sample, sum_n

    ! Donev: We need to give some slack here to allow for roundoff
    ! Not really sure how to best handle this?
    if(sum(p)>(1.0_dp+1e3*epsilon(1.0_dp))) stop "Sum of probabilities must be less than 1"

    sum_p=0
    sum_n=0
    do sample=1, n_samples
       ! One should be careful here to avoid roundoff issues due to some probabilities being zero or summing to not exactly 1
       call BinomialRNG(number=samples(sample), n_trials=N-sum_n, &
               success_prob=min(1.0_dp, p(sample)/max(epsilon(1.0_dp), 1.0_dp-sum_p)), engine=engine)

       ! hack
       do while (samples(sample)<0 .or. samples(sample)>N-sum_n)
          print *, "BinomialRNG ERROR:", N-sum_n, min(1.0_dp, p(sample)/max(epsilon(1.0_dp), 1.0_dp-sum_p)), samples(sample)
          call BinomialRNG(number=samples(sample), n_trials=N-sum_n, &
                  success_prob=min(1.0_dp, p(sample)/max(epsilon(1.0_dp), 1.0_dp-sum_p)), engine=engine)
       end do

       sum_n = sum_n + samples(sample)
       sum_p = sum_p + p(sample)
    end do

 end subroutine

end module


text:
The Davis Legacy Soccer Club (DLSC), formerly known as DYSL, for Davis Youth Soccer League, is a competitive soccer club in Davis. This year, DLSC supports 28 competitive teams from U9 to U18. These teams participate in soccer games and tournaments at the highest level in the US Club League. The DLSC Soccer Fields are located out of town south of Highway 80 to the east of Davis, past the car dealerships. The club is more competitive than Davis AYSO and requires stronger commitments in cost, game play, and parental involvement.

Tournaments
DLSC hosts many renowned tournaments (http://davislegacysoccer.org/tournament/) that attract teams from out of the state, including those from Washington and Canada. The following tournaments are presented this year:
Davis Legacy College Showcase 2010 Girls July 9July 11, 2010
Davis Legacy College Showcase Boys 2010July 1618, 2010
Davis Legacy Super Clasico 2010 August 14/15, 2010
Davis Legacy Premier Cup September 25/26
Davis Legacy Lily Nisen Memorial Tournament October 30/31
Davis Legacy College WINTER Showcase Girls December 11/12, 2010
Davis Legacy College WINTER Showcase Boys December 18/19, 2010

Season
The season depends on the age group. For U9, the season runs from June to early December. For U10 to U13, the season begins in April and runs through January. For the older age groups, U14 and up, season training begins in May continues through April. Tryouts are usually during mid spring for younger age groups and late spring for older. Also, DLSC hosts summer camps open to all youth U5U14, not just for club members, during the summerhttp://davislegacysoccer.org/camps/summercamp/.

Field Renovation
DLSC has built 5 additional fields, making the Legacy Soccer Complex the best soccer venue in northern California. Davis fields have long been recognized as among the best. This expansion helps DLSC grow its programs by sponsoring more teams, hosting larger tournaments and attracting Regional and National events.


text:
C
C $Id: tdsetr.f,v 1.5 2008-07-27 00:17:33 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 TDSETR (PNAM,RVAL)
C
      CHARACTER*(*) PNAM
C
C The subroutine TDSETR may be used to set TDPACK parameters which have
C values of type REAL.
C
C The variables in the following common block define the mapping from
C 3-space to 2-space.
C
        COMMON /TDCOM1/ IH,IT,XM,YM,ZM,XO,YO,ZO,XT,YT,ZT,OE,XE,YE,ZE
        COMMON /TDCOM1/ A1,B1,C1,D1,E1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
        COMMON /TDCOM1/ IS,FV,VL,VR,VB,VT,WL,WR,WB,WT
        SAVE   /TDCOM1/
C
C The variables in the following common block define the parallelogram.
C
        COMMON /TDCOM2/ XACP,YACP,ZACP,XCDX,YCDX,ZCDX,XCDY,YCDY,ZCDY
        SAVE   /TDCOM2/
C
C The variables in the following common block together determine how
C big characters written by TDLBLA will be.
C
        COMMON /TDCOM4/ CSM1,CSM2
        SAVE   /TDCOM4/
C
C The variables in the following common block define the shading type
C and the position of the light source.
C
        COMMON /TDCOM6/ ISHD,XPLS,YPLS,ZPLS
        SAVE   /TDCOM6/
C
C Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
        CALL TDBLDA
C
C Check for an uncleared prior error.
C
        IF (ICFELL('TDSETR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Set the selected parameter.
C
        IF      (PNAM(1:3).EQ.'CS1'.OR.PNAM(1:3).EQ.'cs1') THEN
          CSM1=MAX(0.,RVAL)
        ELSE IF (PNAM(1:3).EQ.'CS2'.OR.PNAM(1:3).EQ.'cs2') THEN
          CSM2=MAX(0.,RVAL)
        ELSE IF (PNAM(1:3).EQ.'FOV'.OR.PNAM(1:3).EQ.'fov') THEN
          FV=RVAL
        ELSE IF (PNAM(1:3).EQ.'HND'.OR.PNAM(1:3).EQ.'hnd') THEN
          IH=MAX(0,MIN(1,INT(RVAL)))
        ELSE IF (PNAM(1:3).EQ.'LSU'.OR.PNAM(1:3).EQ.'lsu') THEN
          XPLS=RVAL
        ELSE IF (PNAM(1:3).EQ.'LSV'.OR.PNAM(1:3).EQ.'lsv') THEN
          YPLS=RVAL
        ELSE IF (PNAM(1:3).EQ.'LSW'.OR.PNAM(1:3).EQ.'lsw') THEN
          ZPLS=RVAL
        ELSE IF (PNAM(1:3).EQ.'SET'.OR.PNAM(1:3).EQ.'set') THEN
          IS=MAX(0,MIN(1,INT(RVAL)))
        ELSE IF (PNAM(1:3).EQ.'SHD'.OR.PNAM(1:3).EQ.'shd') THEN
          ISHD=MAX(0,MIN(1,INT(RVAL)))
        ELSE IF (PNAM(1:3).EQ.'STE'.OR.PNAM(1:3).EQ.'ste') THEN
          IT=MAX(0,MIN(1,INT(RVAL)))
        ELSE IF (PNAM(1:3).EQ.'VPB'.OR.PNAM(1:3).EQ.'vpb') THEN
          VB=MAX(0.,MIN(1.,RVAL))
        ELSE IF (PNAM(1:3).EQ.'VPL'.OR.PNAM(1:3).EQ.'vpl') THEN
          VL=MAX(0.,MIN(1.,RVAL))
        ELSE IF (PNAM(1:3).EQ.'VPR'.OR.PNAM(1:3).EQ.'vpr') THEN
          VR=MAX(0.,MIN(1.,RVAL))
        ELSE IF (PNAM(1:3).EQ.'VPT'.OR.PNAM(1:3).EQ.'vpt') THEN
          VT=MAX(0.,MIN(1.,RVAL))
        ELSE
          CALL SETER ('TDSETR - UNRECOGNIZED PARAMETER NAME',2,1)
          RETURN
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE CLAR2V( N, X, Y, Z, INCX, C, S, INCC )
*
*  -- LAPACK auxiliary routine (version 3.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INCC, INCX, N
*     ..
*     .. Array Arguments ..
      REAL               C( * )
      COMPLEX            S( * ), X( * ), Y( * ), Z( * )
*     ..
*
*  Purpose
*  =======
*
*  CLAR2V applies a vector of complex plane rotations with real cosines
*  from both sides to a sequence of 2-by-2 complex Hermitian matrices,
*  defined by the elements of the vectors x, y and z. For i = 1,2,...,n
*
*     (       x(i)  z(i) ) :=
*     ( conjg(z(i)) y(i) )
*
*       (  c(i) conjg(s(i)) ) (       x(i)  z(i) ) ( c(i) -conjg(s(i)) )
*       ( -s(i)       c(i)  ) ( conjg(z(i)) y(i) ) ( s(i)        c(i)  )
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The number of plane rotations to be applied.
*
*  X       (input/output) COMPLEX array, dimension (1+(N-1)*INCX)
*          The vector x; the elements of x are assumed to be real.
*
*  Y       (input/output) COMPLEX array, dimension (1+(N-1)*INCX)
*          The vector y; the elements of y are assumed to be real.
*
*  Z       (input/output) COMPLEX array, dimension (1+(N-1)*INCX)
*          The vector z.
*
*  INCX    (input) INTEGER
*          The increment between elements of X, Y and Z. INCX > 0.
*
*  C       (input) REAL array, dimension (1+(N-1)*INCC)
*          The cosines of the plane rotations.
*
*  S       (input) COMPLEX array, dimension (1+(N-1)*INCC)
*          The sines of the plane rotations.
*
*  INCC    (input) INTEGER
*          The increment between elements of C and S. INCC > 0.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, IC, IX
      REAL               CI, SII, SIR, T1I, T1R, T5, T6, XI, YI, ZII,
     $                   ZIR
      COMPLEX            SI, T2, T3, T4, ZI
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          AIMAG, CMPLX, CONJG, REAL
*     ..
*     .. Executable Statements ..
*
      IX = 1
      IC = 1
      DO 10 I = 1, N
         XI = REAL( X( IX ) )
         YI = REAL( Y( IX ) )
         ZI = Z( IX )
         ZIR = REAL( ZI )
         ZII = AIMAG( ZI )
         CI = C( IC )
         SI = S( IC )
         SIR = REAL( SI )
         SII = AIMAG( SI )
         T1R = SIR*ZIR - SII*ZII
         T1I = SIR*ZII + SII*ZIR
         T2 = CI*ZI
         T3 = T2 - CONJG( SI )*XI
         T4 = CONJG( T2 ) + SI*YI
         T5 = CI*XI + T1R
         T6 = CI*YI - T1R
         X( IX ) = CI*T5 + ( SIR*REAL( T4 )+SII*AIMAG( T4 ) )
         Y( IX ) = CI*T6 - ( SIR*REAL( T3 )-SII*AIMAG( T3 ) )
         Z( IX ) = CI*T3 + CONJG( SI )*CMPLX( T6, T1I )
         IX = IX + INCX
         IC = IC + INCC
   10 CONTINUE
      RETURN
*
*     End of CLAR2V
*
      END

max_issues_repo_name: cpatrick/calatk

max_forks_repo_name: cpatrick/calatk


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
!
!     determine the critical depth
!
      implicit none
!
      real*8 xflow,rho,b,dg,sqrts0,hk,theta,tth,c1,xflow2,
     &  A,dBBdh,dAdh,BB,dhk
!
      hk=((xflow/(rho*b))**2/(dg*sqrts0))**(1.d0/3.d0)
!
      if(dabs(theta).lt.1.d-10) return
!
!     critical depth for trapezoid, non-rectangular cross section
!
      tth=dtan(theta)
      c1=rho*rho*dg*sqrts0
      xflow2=xflow*xflow
!
      do
         A=hk*(b+hk*tth)
         dBBdh=2.d0*tth
         dAdh=b+hk*dBBdh
         BB=dAdh
         dhk=(xflow2*BB-c1*A**3)/(xflow2*dBBdh-3.d0*c1*A*A*dAdh)
         if(dabs(dhk)/dhk.lt.1.d-3) exit
         hk=hk-dhk
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
!> \file derivnum.f90
!! \BRIEF
!> Module with derivnum subroutine - compute numerical derivatives of carbonate system vars
!> with respect to DIC, Alk, total phosphorus, total silicon, T, S
MODULE mderivnum
CONTAINS

!>    Computes numerical derivatives of standard carbonate system variables
!>    (H+, pCO2, fCO2, CO2*, HCO3- and CO3--, OmegaA, OmegaC) with respect to one given input variable
!!    Input variables are :
!!    temperature, salinity, pressure,
!!    total alkalinity (ALK), dissolved inorganic carbon (DIC),
!!    silica and phosphate concentrations
SUBROUTINE derivnum (dh_dx, dpco2_dx, dfco2_dx, dco2_dx, dhco3_dx,                   &
                     dco3_dx, dOmegaA_dx, dOmegaC_dx,                                &
                     temp, sal, alk, dic, sil, phos, Patm, depth, lat, N, derivar,   &
                     optCON, optT, optP, optB, optK1K2, optKf, optGAS, optS, lon     )
  !   Purpose:
  !     Computes numerical derivatives of standard carbonate system variables
  !     (H+, pCO2, fCO2, CO2*, HCO3- and CO3--, OmegaA, OmegaC) with respect to one given input variable
  !     FROM:
  !     temperature, salinity, pressure,
  !     total alkalinity (ALK), dissolved inorganic carbon (DIC),
  !     silica and phosphate concentrations

  !     INPUT variables:
  !     ================
  !     Patm    = atmospheric pressure [atm]
  !     depth   = depth [m]     (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure)
  !             = pressure [db] (with optP='db')
  !     lat     = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m')
  !             = dummy array (unused when optP='db')
  !     temp    = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not in situ temp)
  !             = in situ   temperature [degrees C] (with optT='Tinsitu', e.g., for data)
  !             = conservative temperature [degrees C] (with optT='Tcsv')
  !     sal     = salinity in [psu]
  !     alk     = total alkalinity in [eq/m^3] with optCON = 'mol/m3'
  !             =               [eq/kg]  with optCON = 'mol/kg'
  !     dic     = dissolved inorganic carbon [mol/m^3] with optCON = 'mol/m3'
  !             =                            [mol/kg]  with optCON = 'mol/kg'
  !     sil     = silica    [mol/m^3] with optCON = 'mol/m3'
  !             =           [mol/kg]  with optCON = 'mol/kg'
  !     phos    = phosphate [mol/m^3] with optCON = 'mol/m3'
  !             =           [mol/kg]  with optCON = 'mol/kg'
  !     derivar = 3-lowercase-character identifier of input variable with respect to which derivative is requested
  !               possibilities are 'alk', 'dic', 'pho', 'sil', 'tem', 'sal'
  !               and dissociation constants 'k0', 'k1', 'k2', 'kw', 'kb', 'ka', 'kc'
  !               and total dissolved inorganic boron 'bt'
  !               The last two of the k's (constants) are the solubility products for Aragonite and Calcite
  !
  !     INPUT options:
  !     ==============
  !     -----------
  !     optCON: choose input & output concentration units - mol/kg (data) vs. mol/m^3 (models)
  !     -----------
  !       -> 'mol/kg' for DIC, ALK, sil, & phos given on mokal scale, i.e., in mol/kg  (std DATA units)
  !       -> 'mol/m3' for DIC, ALK, sil, & phos given in mol/m^3 (std MODEL units)
  !     -----------
  !     optT: choose in situ vs. potential temperature as input
  !     ---------
  !     NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature)
  !       -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed)
  !       -> 'Tcsv' means input is Conservative Temperature (in situ Temp "tempis" is computed)
  !       -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed)
  !     ---------
  !     optP: choose depth (m) vs pressure (db) as input
  !     ---------
  !       -> 'm'  means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed)
  !       -> 'db' means "depth" input is already in situ pressure [db], not m (thus p = depth)
  !     ---------
  !     optB: choose total boron formulation - Uppström (1974) vs. Lee et al. (2010)
  !     ---------
  !       -> 'u74' means use classic formulation of Uppström (1974) for total Boron
  !       -> 'l10' means use newer   formulation of Lee et al. (2010) for total Boron
  !     ---------
  !     optK1K2:
  !     ---------
  !       -> 'l'   means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007)
  !                **** BUT this should only be used when 2 < T < 35 and 19 < S < 43
  !       -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007)
  !                **** Valid for 0 < T < 50°C and 1 < S < 50 psu
  !     ----------
  !     optKf:
  !     ----------
  !       -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007)
  !               **** BUT Valid for  9 < T < 33°C and 10 < S < 40.
  !       -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994)
  !     -----------
  !     optGAS: choose in situ vs. potential fCO2 and pCO2
  !     ---------
  !       PRESSURE corrections for K0 and the fugacity coefficient (Cf)
  !       -> 'Pzero'   = 'zero order' fCO2 and pCO2 (typical approach, which is flawed)
  !                      considers in situ T & only atm pressure (hydrostatic=0)
  !       -> 'Ppot'    = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
  !                      considers potential T & only atm pressure (hydrostatic press = 0)
  !       -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure)
  !                      considers in situ T & total pressure (atm + hydrostatic)
  !     ---------
  !     optS: choose practical [psu] or absolute [g/kg] salinity as input
  !     ----------
  !       -> 'Sprc' means input is practical salinity according to EOS-80 convention
  !       -> 'Sabs' means input is absolute salinity according to TEOS-10 convention (practical sal. will be computed)
  !     ---------
  !     lon:  longitude in degrees East
  !     ----------
  !        Optional, it may be used along with latitude when optS is "Sabs".
  !        Then, they are parameters for conversion from Absolute to Practical Salinity.
  !
  !        When seawater is not of standard composition, Practical Salinity alone is not sufficient
  !        to compute Absolute Salinity and vice-versa. One needs to know the chemical composition,
  !        mainly silicate and nitrate concentration. When those concentrations are unknown and 'lon' and 'lat'
  !        are given, absolute salinity conversion is based on WOA silicate concentration at given location.
  !
  !        Alternative when optS is 'Sabs' :
  !        -------------------------------
  !        When silicate and phosphate concentrations are known, nitrate concentration is inferred from phosphate
  !        (using Redfield ratio), then practical salinity is computed from absolute salinity,
  !        total alcalinity (alk), DIC (dic), silicate (sil) and phosphate (phos).
  !        In that case, do not pass optional parameter 'lon'.
  !
  !        When neither chemical composition nor location are known, an arbitrary geographic point is chosen:
  !        mid equatorial Atlantic. Note that this implies an error on computed practical salinity up to 0.02 psu.
  !        In that case, do pass parameter 'lon' and set each of its elements to 1.e20.
  !     ---------

  !     OUTPUT variables:
  !     =================
  !     dh_dx      = derivative of ion [H+] concentration on total scale in [mol/kg]
  !     dpco2_dx   = derivative of CO2 partial pressure (uatm)
  !     dfco2_dx   = derivative of CO2 fugacity (uatm)
  !     dco2_dx    = derivative of aqueous CO2 concentration in [mol/kg] or [mol/m^3] depending on optCON
  !     dhco3_dx   = derivative of bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] depending on optCON
  !     dco3_dx    = derivative of carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] depending on optCON
  !     dOmegaA_dx = derivative of Omega for aragonite, i.e., the aragonite saturation state
  !     dOmegaC_dx = derivative of Omega for calcite, i.e., the   calcite saturation state

  USE msingledouble
  USE mvars

  IMPLICIT NONE

! Input variables
  !>     number of records
!f2py intent(hide) n
  INTEGER, INTENT(in) :: N
  !> either in situ temperature (when optT='Tinsitu', typical data)
  !! OR potential temperature (when optT='Tpot', typical models) [degree C]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: temp
  !> salinity [psu] or [g/kg]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: sal
  !> total alkalinity in [eq/m^3] (when optCON = 'mol/m3') OR in [eq/kg]  (when optCON = 'mol/kg')
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: alk
  !> dissolved inorganic carbon in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg')
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: dic
  !> SiO2 concentration in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg')
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: sil
  !> phosphate concentration in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg')
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: phos
  !> atmospheric pressure [atm]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: Patm
  !> depth in \b meters (when optP='m') or \b decibars (when optP='db')
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: depth
  !> latitude [degrees north]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: lat
  ! 3-character identifier of input variable with respect to which derivative is requested
  ! ('alk', 'dic', 'pho', 'sil', 'tem', 'sal'
  CHARACTER(3), INTENT(in) ::  derivar

  !> choose either \b 'mol/kg' (std DATA units) or \b 'mol/m3' (std MODEL units) to select
  !! concentration units for input (for alk, dic, sil, phos) & output (co2, hco3, co3)
  CHARACTER(6), INTENT(in) :: optCON
  !> choose \b 'Tinsitu' for in situ temperature or \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models)
  CHARACTER(7), INTENT(in) :: optT
  !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models)
  CHARACTER(2), INTENT(in) :: optP
  !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010).
  !! The 'l10' formulation is based on 139 measurements (instead of 20),
  !! uses a more accurate method, and
  !! generally increases total boron in seawater by 4%
!f2py character*3 optional, intent(in) :: optB='u74'
  CHARACTER(3), OPTIONAL, INTENT(in) :: optB
  !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979)
!f2py character*2 optional, intent(in) :: optKf='pf'
  CHARACTER(2), OPTIONAL, INTENT(in) :: optKf
  !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010)
!f2py character*3 optional, intent(in) :: optK1K2='l'
  CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2
  !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction)
  !! 'Ppot'    - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
  !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf
  !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean
!f2py character*7 optional, intent(in) :: optGAS='Pinsitu'
  CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS
!f2py character*7 optional, intent(in) :: optS='Sprc'
  CHARACTER(4), OPTIONAL, INTENT(in) :: optS
  !> longitude [degrees east]
!!!f2py real(8) optional, intent(in), dimension(n) :: lon = -25.
  REAL(kind=rx), OPTIONAL, INTENT(in), DIMENSION(N) :: lon
!f2py optional :: lon = -25.

! Output variables:
  !> derivative of H on the total scale
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dh_dx
  !> derivative of CO2 partial pressure [uatm]
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dpco2_dx
  !> derivative of CO2 fugacity [uatm]
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dfco2_dx
  !> derivative of aqueous CO2* concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dco2_dx
  !> derivative of (HCO3-) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dhco3_dx
  !> derivative of (CO3--) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dco3_dx
  !> derivative of Omega for aragonite, i.e., the aragonite saturation state
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dOmegaA_dx
  !> derivative of Omega for calcite, i.e., the calcite saturation state
  REAL(kind=rx), INTENT(OUT), DIMENSION(N) :: dOmegaC_dx

! Local variables
  !> pH on the total scale
  REAL(kind=rx), DIMENSION(1,2) :: ph
  !> ion (H+) concentration on the total scale
  REAL(kind=rx), DIMENSION(1,2) :: h
  !> CO2 partial pressure [uatm]
  REAL(kind=rx), DIMENSION(1,2) :: pco2
  !> CO2 fugacity [uatm]
  REAL(kind=rx), DIMENSION(1,2) :: fco2
  !> aqueous CO2* concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), DIMENSION(1,2) :: co2
  !> bicarbonate ion (HCO3-) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), DIMENSION(1,2) :: hco3
  !> carbonate ion (CO3--) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON
  REAL(kind=rx), DIMENSION(1,2) :: co3
  !> Omega for aragonite, i.e., the aragonite saturation state
  REAL(kind=rx), DIMENSION(1,2) :: OmegaA
  !> Omega for calcite, i.e., the calcite saturation state
  REAL(kind=rx), DIMENSION(1,2) :: OmegaC
  !> Revelle factor, i.e., dpCO2/pCO2 / dDIC/DIC
  REAL(kind=rx), DIMENSION(1) :: BetaD
  !> in-situ density of seawater; rhoSW = f(s, t, p) in [kg/m3]
  REAL(kind=rx), DIMENSION(1) :: rhoSW
  !> pressure [decibars]; p = f(depth, latitude) if computed from depth [m] (when optP='m') OR p = depth [db] (when optP='db')
  REAL(kind=rx), DIMENSION(1) :: p
  !> in-situ temperature \b [degrees C]
  REAL(kind=rx), DIMENSION(1) :: tempis

  ! Arrays that are copies of input parameters
  REAL(kind=rx), DIMENSION(1) :: atemp
  REAL(kind=rx), DIMENSION(1) :: asal
  REAL(kind=rx), DIMENSION(1) :: aalk
  REAL(kind=rx), DIMENSION(1) :: adic
  REAL(kind=rx), DIMENSION(1) :: asil
  REAL(kind=rx), DIMENSION(1) :: aphos
  REAL(kind=rx), DIMENSION(1) :: aPatm
  REAL(kind=rx), DIMENSION(1) :: adepth
  REAL(kind=rx), DIMENSION(1) :: alat

  ! value of small delta to apply to input variable when computing numerical derivative
  ! it is actually the ratio of delta relative to input variable value
  REAL(kind=rx) :: rel_delta_x
  ! Value of input variable with respect to which we derive
  REAL(kind=rx) :: input_value
  REAL(kind=rx), DIMENSION(1) :: ainput1, ainput2
  REAL(kind=rx) :: abs_delta, dX

! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007)
  CHARACTER(3) :: opB
  CHARACTER(3) :: opK1K2
  CHARACTER(2) :: opKf
  CHARACTER(7) :: opGAS
  CHARACTER(4) :: opS

  INTEGER :: i
  LOGICAL :: deriv_K    !  if one derives with respect to a dissoc constant
  INTEGER :: var_index  ! index of variable with respect to which one derives

! Approximate values for K0, K1, K2, Kb, Kspa and Kspc
! They will be used to compute an absolute perturbation value on these constants
  REAL,PARAMETER :: K_values(7) = (/0.034, 1.2e-06, 8.3e-10, 2.1e-09, 6.1e-14, 6.7e-07, 4.3e-07/)

! Multiplicative conversion factor to convert from umol/kg to optCON units (mol/kg or mol/m3)
  REAL(kind=rx) :: xfac

! Reference value (typically global surface means, for calculating consistent absolute derivatives)
  REAL(kind=rx) :: ref_value

! Set defaults for optional arguments (in Fortran 90)
! Note:  Optional arguments with f2py (python) are set above with
!        the !f2py statements that precede each type declaraion
  IF (PRESENT(optB)) THEN
    opB = optB
  ELSE
!   Default is Uppstrom (1974) for total boron
    opB = 'u74'
  ENDIF
  IF (PRESENT(optKf)) THEN
    opKf = optKf
  ELSE
!   Default is Perez & Fraga (1987) for Kf
    opKf = 'pf'
  ENDIF
  IF (PRESENT(optK1K2)) THEN
    opK1K2 = optK1K2
  ELSE
!   Default is Lueker et al. 2000) for K1 & K2
    opK1K2 = 'l'
  ENDIF
  IF (PRESENT(optGAS)) THEN
    opGAS = optGAS
  ELSE
    opGAS = 'Pinsitu'
  ENDIF
  IF (PRESENT(optS)) THEN
    opS = optS
  ELSE
    opS = 'Sprc'
  ENDIF

  SELECT CASE (derivar)
      CASE ('alk')
          var_index = 1
          deriv_K   = .FALSE.
      CASE ('dic')
          var_index = 2
          deriv_K   = .FALSE.
      CASE ('pho')
          var_index = 3
          deriv_K   = .FALSE.
      CASE ('sil')
          var_index = 4
          deriv_K   = .FALSE.
      CASE ('tem')
          var_index = 5
          deriv_K   = .FALSE.
      CASE ('sal')
          var_index = 6
          deriv_K   = .FALSE.
      CASE ('k0')
          var_index = 1
          deriv_K   = .TRUE.
      CASE ('k1')
          var_index = 2
          deriv_K   = .TRUE.
      CASE ('k2')
          var_index = 3
          deriv_K   = .TRUE.
      CASE ('kb')
          var_index = 4
          deriv_K   = .TRUE.
      CASE ('kw')
          var_index = 5
          deriv_K   = .TRUE.
      CASE ('ka')
          var_index = 6
          deriv_K   = .TRUE.
      CASE ('kc')
          var_index = 7
          deriv_K   = .TRUE.
      CASE ('bt')
          ! For Bt (propotional to S): call vars_pertK but treat slightly differently than other constants
          var_index = 8
          deriv_K   = .FALSE.
      CASE DEFAULT
          PRINT *,"derivar must be 3-char variable: 'alk', 'dic', 'pho', 'sil', 'tem', or 'sal'"
          PRINT *," or a 2-char dissoc. constant name : 'k0', 'k1', 'k2', 'kb', 'kw', 'ka', 'kc', 'bt'"
          STOP
  END SELECT

  DO i = 1, N

    ! Copy input values
    atemp(1)  = temp(i)
    asal(1)   = sal(i)
    aalk(1)   = alk(i)
    adic(1)   = dic(i)
    asil(1)   = sil(i)
    aphos(1)  = phos(i)
    aPatm(1)  = Patm(i)
    adepth(1) = depth(i)
    alat(1)   = lat(i)

    IF (deriv_K) THEN
        ! Choose value of absolute perturbation
        abs_delta = K_values(var_index) * 1.d-3   ! 0.1 percent of Kx value
        dx = 2 * abs_delta

        ! Call routine vars_pertK() with request to increase one K value
        call vars_pertK(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),      &
                        OmegaA(:,1), OmegaC(:,1),                                          &
                        atemp, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                        var_index, -abs_delta,                                             &
                        optCON, optT, optP, opB, opK1K2, opKf, opGAS,  optS, lon           )
        ! Call routine vars_pertK() with request to decrease the same K value
        call vars_pertK(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),      &
                        OmegaA(:,2), OmegaC(:,2),                                          &
                        atemp, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                        var_index, abs_delta,                                              &
                        optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon            )

    ELSE    ! we do NOT derive w/ respect to one of the  Ks

        ! Simple conversion factors
       IF (optCON .EQ. 'mol/kg') THEN
          xfac = 1.d-6    !Factor to convert umol/kg to mol/kg
       ELSE
          xfac = 1.025d-3 !Factor to convert umol/kg to mol/m3
       ENDIF

        ! Select input value to vary slightly
        ! values for individual rel_delta_x determined by minimizing difference
        ! between numerical derivatives and mocsy-dnad automatic derivatives
        ! (mocsy_dnad), which gives results as good as analytic solutions
        SELECT CASE (var_index)
            CASE (1)
                input_value = alk(i)
                ref_value = 2300.0d0 * xfac ! global surface mean in optCON units
                rel_delta_x = 1.0d-6
            CASE (2)
                input_value = dic(i)
                ref_value = 2000.0d0 * xfac     ! global surface mean in optCON units
                rel_delta_x = 1.0d-6
            CASE (3)
                input_value = phos(i)
                ref_value = 0.5d0 * xfac   ! global surface mean (Orr et al., 2017) in optCON units
                rel_delta_x = 1.0d-3
            CASE (4)
                input_value = sil(i)
                ref_value = 7.5d0 * xfac   ! global surface mean (Orr et al., 2017) in optCON units
                rel_delta_x = 1.0d-3
            CASE (5)
                input_value = temp(i)
                ref_value = 18.0d0 ! global surface mean (C)
                rel_delta_x = 1.0d-4
            CASE (6)
                input_value = sal(i)
                ref_value = 35.0d0 ! global surface mean (C)
                rel_delta_x = 1.0d-4
            CASE (8)
                ! For Bt (propotional to salinity)
                input_value = 0.0002414d0 * (sal(i)  / 1.80655d0) / 10.811d0
                ref_value   = 0.0002414d0 * (35.0d0 / 1.80655d0) / 10.811d0
                rel_delta_x = 1.0d-3
        END SELECT

        ! Determine two slightly different values of selected input value
        abs_delta  = ref_value * rel_delta_x
        ainput1(1) = input_value - abs_delta
        ainput2(1) = input_value + abs_delta
        ! Compute total absolue delta
        dx = ainput2(1) - ainput1(1)

        ! Call routine vars() twice with two slightly different values of selected input
        SELECT CASE (var_index)
            CASE (1)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, ainput1, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, ainput2, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon               )
            CASE (2)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, ainput1, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, ainput2, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon               )
            CASE (3)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, adic, asil, ainput1, aPatm, adepth, alat, 1,       &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, adic, asil, ainput2, aPatm, adepth, alat, 1,       &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon             )
            CASE (4)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, adic, ainput1, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                atemp, asal, aalk, adic, ainput2, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon             )
            CASE (5)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                ainput1, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,       &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                ainput2, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,       &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon             )
            CASE (6)
                call vars(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),               &
                                OmegaA(:,1), OmegaC(:,1), BetaD, rhoSW, p, tempis,                    &
                                atemp, ainput1, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon, verbose=.false. )
                call vars(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),               &
                                OmegaA(:,2), OmegaC(:,2), BetaD, rhoSW, p, tempis,                    &
                                atemp, ainput2, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon             )
            CASE (8)
                ! For Bt (must treat with vars_pertK, but not quite like other 'constants' to account for proportionality to S)
                ! Call routine vars_pertK() with request to increase one Bt value
                call vars_pertK(ph(:,1), pco2(:,1), fco2(:,1), co2(:,1), hco3(:,1), co3(:,1),      &
                                OmegaA(:,1), OmegaC(:,1),                                          &
                                atemp, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                var_index, -abs_delta,                                             &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon            )
                ! Call routine vars_pertK() with request to decrease the same Bt value
                call vars_pertK(ph(:,2), pco2(:,2), fco2(:,2), co2(:,2), hco3(:,2), co3(:,2),      &
                                OmegaA(:,2), OmegaC(:,2),                                          &
                                atemp, asal, aalk, adic, asil, aphos, aPatm, adepth, alat, 1,      &
                                var_index, abs_delta,                                              &
                                optCON, optT, optP, opB, opK1K2, opKf, opGAS, optS, lon            )
        END SELECT
    ENDIF

    ! H+ concentration from ph
    h(1,1) = 10**(-ph(1,1))
    h(1,2) = 10**(-ph(1,2))

    ! Compute derivatives (centered differences)
    dh_dx(i)      = (h(1,2)    - h(1,1))    / dx
    dpco2_dx(i)   = (pco2(1,2) - pco2(1,1)) / dx
    dfco2_dx(i)   = (fco2(1,2) - fco2(1,1)) / dx
    dco2_dx(i)    = (co2(1,2)  - co2(1,1))  / dx
    dhco3_dx(i)   = (hco3(1,2) - hco3(1,1)) / dx
    dco3_dx(i)    = (co3(1,2)  - co3(1,1))  / dx
    dOmegaA_dx(i) = (OmegaA(1,2) - OmegaA(1,1)) / dx
    dOmegaC_dx(i) = (OmegaC(1,2) - OmegaC(1,1)) / dx
  END DO

  RETURN
END SUBROUTINE derivnum
END MODULE mderivnum

max_stars_repo_name: tomaslovato/mocsy

max_issues_repo_name: tomaslovato/mocsy

max_forks_repo_name: tomaslovato/mocsy


text:
program compare_floats ! from urbanjost
implicit none
! below might test as equal
! real :: A(3)=[1.0,2.0,3.0], B(3)=[1.0,2.0,2.9999999999999]
real :: A(3)=[1.0,2.0,3.0], B(3)=[1.0,2.0,2.999999] ! this should be close enough
real :: tolerance=0.00001  ! just a sample tolerance
if (all(A==B)) then  ! testing for exact matches can be problematic
   write(*,*) "A and B are equal"
else if (all(abs(A-B) < tolerance)) then
   write(*,*) "A and B are close enough to equal"
else
   write(*,*) "A and B are NOT equal"
end if
! output: A and B are close enough to equal
end program compare_floats

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
c=======================================================================
c Program to call the  GTN material law using different material
c parameters of VM and GTN to compare its behaviour in a isochoric tension
c-----------------------------------------------------------------------
c In order to call this program in the same folder muste be a csv file of
c at least one file  and maximum 5, with 15 parameters, 3 linear elasticity,
c  2 VM plasticity and 10 fort GTN model. See more information in the
c documentation.
c-----------------------------------------------------------------------
c Version 0.9.1
c Oct 2021
c=======================================================================
        program test_isoc_ten
        use tensor_operations

        implicit none
        double precision :: te_n1
        double precision :: n_ampl=2.0; !2      ! n_ampl>1
        double precision ::  xE, xnu, sigma_y0, kappa
        double precision :: mu
        integer :: stat !Variables t print out information
        integer :: ltype, type_linespace !Selector to compute C and linespace respectively
        integer :: ttype !Selector to compute C
        integer :: n_step !number of steps
        integer,dimension(5) :: posi_last_time=(/2,3,3,4,7/)
        integer :: i !iterators
        integer, dimension(3):: options
        double precision, dimension (7) :: t !Time range
        double precision :: dt!step size for strain ramp
        double precision, dimension (3) :: lam ! Strain range
        double precision, dimension(:), allocatable:: time
        double precision, dimension(:), allocatable:: e11 !11 strain
        double precision, dimension (:,:), allocatable:: sdv, sdv_1
        double precision, dimension (21) ::sdvup, sdvup_1
        double precision, dimension (:), allocatable::s11, s22, s33, !normal stresses
     &       s11_1, s22_1, s33_1
        double precision, dimension (6):: e6, De6, sigma, sigma_1 !Total strain and total stress vectors
        double precision, dimension(6,6) :: A66, A66_1
        double precision:: mat_param(5,15)
        integer:: n1, n2, n_test, matL
        character:: file_name*13, file_name2*30, file_na_ex*27  !******DON'T CHANGE IT *******
c       double precision, dimension(10):: e11

c   xE             : Young's modulus
c   xnu            : Poisson's ratio
c   xsigy0         : initial yield stress
c   xmu            : elastic shear modulus
c   xk             : bulk modulus
c   kappa          : bulk modulus K
c   ltype          : type of load
c   te_n1          : temporal variable 1
c   posi_last_time : vector to help to define the ramp size
c   t              : vector with time range for all the ramp options
c   lamp           : load range
c   n1             : number of test (rows input file)
c   n_test         : number of test
c   n2             : number of parameters (columnss input file)
c   mat_param      : matrix to store imported parameters for test
c   matL           : auxiliary integer to run different material laws
c   type_linespace : Selector functionality of linespace
c   ttype          : Selector tangent stiffness coputation method
c   dt             : delta time increment
c   time            :time vector for the ramp
c   e11            : strain 11 ramp
c   n_step         : number of steps ramp
c   sdv , sdv_1    : tensors to store the state variables history,
c                    of material law 0 and 1 respectively, a row for
c                    each step. The vector form comes already from the
c                    mterial subroutine
c   sdvup          : updated ISV from the material routine
c   A66            : ATS in 6x6 matrix
c   s11, s22, s33  : vetors to stroe the hystorey of the normal stresses
c   e6             : current total strain tensor defined for an specified time
c   De6            : standard strain increment between time increment
c                    defined for this test
c   file_name      : Name of the file to write results,******DON'T CHANGE IT *******
c   stat           : Variable for the state of the opened file to write
c   options        : flags for options in the material routine
c=======================================================================
c                  Define material parameters
c=======================================================================
c-----------Import data with material parameters for test---------------
        n1=size(mat_param,1)
        n2=size(mat_param,2)
        file_name2="4_Test_lin_elas0/dat_LE_01.csv"
        !file_name ="dat_LE_01.csv"
        call read3 (mat_param, n1,n2, file_name2)
        do matL=0,1
        do n_test=1, n1
                print*, "test", n_test
c-------------Compute elastic constant---------------------------------
        xE       = mat_param(n_test,1)
        xnu      = mat_param(n_test,2)
        sigma_y0 = mat_param(n_test,3)
        mu       = xE/(2.0*(1.0+xnu))
        kappa    = xE/(3.0*(1.0-2.0*xnu))
c=======================================================================
c                       define loading ramp
c=======================================================================
c----------------Select the time (load) amplitud------------------------
c     1: linear ramp of load
c     2: Linear load and unload toero
        ltype=1

        te_n1=0.0
        if (ltype==1) then
            t(1:2)=(/0.0, 10.0/) ! 10
c            %This is like a load which is related with sigma_y0 and an integer n_ampl
            lam(1:2)=(/te_n1 , (n_ampl*sigma_y0/2/mu)/2.1 /) !3.0 !total Strain considering the value of mu to onvert from load to strain
        elseif (ltype==2) then
            t(1:3)=(/0.0, 5.0, 10.0/)
            lam(1:2)=(/0.0, 1.59155e-3/3.1 /)
        end if

c-------- computation of number of steps and e11 strain vs time--------
        dt=0.1
        n_step=nint((t(posi_last_time(ltype))-t(1))/dt)
        type_linespace=1
c-------------------Define time ramp------------------------------------
        allocate (time(n_step+1))
        call linespace(t(1),t(posi_last_time(ltype)),time,
     &      type_linespace)

c--------------Define the strain e11 reference load ramp------------------------------
        allocate (e11(n_step+1))
        call loading2 (ltype,posi_last_time ,dt, t, lam, e11)
        !print*, "ramp", e11
c=======================================================================
c       initiate variables for post-processing
c=======================================================================
c-------initiate the tensor sdv according the number of steps-----------
        allocate (sdv(21,n_step+2)) ; allocate (sdv_1(21,n_step+2))
        sdv(20,1)=mat_param(n_test,9); sdv_1(20,1)=mat_param(n_test,9)
        print*,  "initial sdv", sdv(:,1)

c----------Initiate vectors  s11, s22, s33 (normal stress).-------------
        allocate (s11 (n_step+2)); allocate (s11_1 (n_step+2))
        allocate (s22 (n_step+2)); allocate (s22_1 (n_step+2))
        allocate (s33 (n_step+2)); allocate (s33_1 (n_step+2))

c=======================================================================
c   Start loop to compute material stress behavior due to strain ramp
c=======================================================================
        ttype = 0 ! 0: analytical 1: numerical
        do i=1,n_step
c----Define current driven strain at the begining of an iteration-------
            e6=(/1.0, -0.5, -0.5, 0.0, 0.0, 0.0/)*e11(i);
            De6=(/1.0, -0.5, -0.5, 0.0, 0.0, 0.0/)*(e11(i)-e11(i-1))
c---------Call the material law to compute the stress, ATS, and ISV-----
c        outputs
c        sigma: stress n+1 as a vector
c        A66: algritmic tangent (4 grade tensor as a 6x6 matrix) at n+1
c        sdvup: Accumulated Internal state variables at n+1
c
c        Inputs
c        e      : current total strain
c        De6    : strain increment between time increment
c        sdv(i) : current internal state variables
c        options: flags for options in the material routine
c        mat_param: material parameters according to the test
         options(1) = 0 ! ATS selectror
         options(2) = 0 ! Stress Case
         options(3) = matL ! Material Law
         call kGTN (e6,De6,sdv(:,i),options, mat_param(n_test,:), sigma,
     &        A66, sdvup)

c=======================================================================
c store the updated variables to start the next step and for the
c post-processing
c=======================================================================
         sdv(:,i+1) = sdvup

c        We are only interested in the normal stresses
         s11(i+1)=sigma(1)
         s22(i+1)=sigma(2)
         s33(i+1)=sigma(3)

        end do

c=======================================================================
c        Export the results to a cs file for post-processing
c=======================================================================
c+++++++++++++++++++++++Material law 1++++++++++++++++++++++++++++++++++
        if (matL==0) then
            file_na_ex='4_Test_lin_elas0/dVM0'//char(48+n_test)//'.csv'
            !print*, file_na_ex
c---------------Erase old files with the same name-----------------------
            open(unit=1, iostat=stat, file=file_na_ex,
     $             status='old')
                 if (stat == 0) close(1, status='delete')

            ! output data into a file as follow:
                 !e11 = Driving strain
                 !s11 = stress 11
                 !s22 = stress 22
                 !s33 = stress 33
                 !sdv (1)= plastic strain 11
                 !sdv (2)= plastic strain 22
                 !sdv (3)= plastic strain 33
                 !sdv (20)= porosity
                 !sdv (21)= microscopy equivalent plastic strain
c------------------------Export the data--------------------------------
            open(1, file = file_na_ex, status = 'new')
            do i=1,n_step+1
               write(1,*) e11(i), s11(i), s22(i), s33 (i),
     &             sdv(1,i),sdv(2,i),sdv(3,i), sdv(20,i),sdv(21,i)
c     &             ,s11_1(i), s22_1(i), s33_1(i),
c     &             sdv_1(1,i),sdv_1(2,i),sdv_1(3,i), sdv_1(20,i),sdv_1(21,i)
            end do
            close(1)
        elseif (matL==1) then
c+++++++++++++++++++++++Material law 2++++++++++++++++++++++++++++++++++
            file_na_ex='4_Test_lin_elas0/dGN0'//char(48+n_test)//'.csv'
            !print*, file_na_ex
c---------------Erase ol files with the same name-----------------------
            open(unit=1, iostat=stat, file=file_na_ex,
     $             status='old')
                 if (stat == 0) close(1, status='delete')
c------------------------Exprot the data--------------------------------
            open(1, file = file_na_ex, status = 'new')
            do i=1,n_step+1
               write(1,*) e11(i), s11(i), s22(i), s33 (i),
     &             sdv(1,i),sdv(2,i),sdv(3,i), sdv(20,i),sdv(21,i)
            end do
            close(1)
        end if
c-----------------Deallocate loop allocatable variables-----------------
        deallocate (s11); deallocate (s22); deallocate (s33)
        deallocate (s11_1); deallocate (s22_1); deallocate (s33_1)
        deallocate (e11)
        deallocate (sdv); deallocate (sdv_1)
        deallocate (time) !vector with the steps
       end do
       end do
      end program test_isoc_ten

ext: for


text:
module RefinementModule

use NumberKindsModule
use STDIntVectorModule
use OutputWriterModule
use LoggerModule
use ParticlesModule
use FieldModule
use EdgesModule
use FacesModule
use PlaneGeomModule
use SphereGeomModule
use PolyMesh2dModule

implicit none

private
public RefineSetup, New, Delete, FlagFunction
public ScalarIntegralRefinement, ScalarVariationRefinement
public IterateMeshRefinementForFlowMap
public IterateMeshRefinementOneVariable
public IterateMeshRefinementOneVariableAndFlowMap
public IterateMeshRefinementTwoVariables

!----------------
! types and module variables
!----------------
!
type RefineSetup
	logical(klog), pointer :: refineFlag(:) => null()
	contains
		final :: deletePrivate
end type
!
!----------------
! interfaces
!----------------
!
interface New
	module procedure newPrivate
end interface

interface Delete
	module procedure deletePrivate
end interface

interface
	function FlagFunction( mesh, dataField, tolerance, faceIndex )
		use NumberKindsModule
		use PolyMesh2dModule
		use FieldModule
		logical(klog) :: FlagFunction
		type(PolyMesh2d), intent(in) :: mesh
		type(Field), intent(in) :: dataField
		real(kreal), intent(in) :: tolerance
		integer(kint), intent(in) :: faceIndex
	end function
end interface

!
!----------------
! Logging
!----------------
!
logical(klog), save :: logInit = .FALSE.
type(Logger) :: log
character(len=28), save :: logKey = 'Refine'
character(len=MAX_STRING_LENGTH) :: logString
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL

contains

!
!----------------
! public methods
!----------------
!
subroutine newPrivate(self, nMaxFaces )
	type(RefineSetup), intent(out) :: self
	integer(kint), intent(in) :: nMaxFaces

	if (.NOT. logInit ) call InitLogger(log, procRank)

	allocate(self%refineFlag(nMaxFaces))
	self%refineFlag = .FALSE.
end subroutine

subroutine deletePrivate(self)
	type(RefineSetup), intent(inout) :: self
	if ( associated(self%refineFlag)) deallocate(self%refineFlag)
end subroutine

subroutine DivideFlaggedFaces( self, mesh )
	type(RefineSetup), intent(in) :: self
	type(PolyMesh2d), intent(inout) :: mesh
	!
	integer(kint) :: i
	integer(kint) :: nFacesIn
	integer(kint) :: spaceLeft
	logical(klog) :: limitReached
	integer(kint) :: flagCount, refineCount

	nFacesIn = mesh%faces%N
	flagCount = count(self%refineFlag)
	spaceLeft = mesh%faces%N_Max - mesh%faces%N
	if ( spaceLeft / 4 < flagCount ) then
		call LogMessage(log, WARNING_LOGGING_LEVEL, trim(logKey)//" DivideFlaggedFaces WARNING : ", &
			"not enough memory to continue AMR.")
		return
	endif

	call LogMessage(log, DEBUG_LOGGING_LEVEL, trim(logKey)//" entering :", " DivideFlaggedFaces")

	refineCount = 0
	limitReached = .FALSE.
	do i = 1, nFacesIn
		if ( self%refineFlag(i) ) then
			if ( mesh%initNest + mesh%amrLimit > CountParents(mesh%faces, i) ) then
				if ( mesh%faceKind == TRI_PANEL) then
					call DivideTriFace( mesh%faces, i, mesh%particles, mesh%edges)
				elseif ( mesh%faceKind == QUAD_PANEL ) then
					call DivideQuadFace( mesh%faces, i, mesh%particles, mesh%edges)
				endif
				refineCount = refineCount + 1
			else
				limitReached = .TRUE.
			endif
			self%refineFlag(i) = .FALSE.
		endif
	enddo
	if ( limitReached ) then
		call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" DivideFlaggedFaces REMARK : ", &
			"AMR limit reached. Some flagged faces were not divided.")
	endif
	write(logString,'(2(A,I6),A)') "refined ", refineCount, " of ", flagCount, " flagged faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" DivideFlaggedFaces ", logString)
end subroutine

subroutine IterateMeshRefinementOneVariable(self, aMesh, aField, flagFn, tol, desc, nParticlesStart, nParticlesEnd )
	type(RefineSetup), intent(inout) :: self
	type(PolyMesh2d), intent(inout) :: aMesh
	type(Field), intent(in) :: aField
	procedure(FlagFunction) :: flagFn
	real(kreal), intent(in) :: tol
	character(len=*), intent(in) :: desc
	integer(kint), intent(out) :: nParticlesStart
	integer(kint), intent(out) :: nParticlesEnd
	!
	integer(kint) :: counter
	integer(kint) :: i

	counter = 0
	nParticlesStart = aMesh%particles%N
	do i = 1, aMesh%faces%N
		if ( .NOT. aMesh%faces%hasChildren(i) ) then
			if ( flagFn(aMesh, aField, tol, i) ) then
				self%refineFlag(i) = .TRUE.
				counter = counter + 1
			endif
		endif
	enddo

	write(logString, '(2A,I8,A)') trim(desc), ": flagged ", counter, " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	call DivideFlaggedFaces( self, aMesh)
	nParticlesEnd = aMesh%particles%N
end subroutine

subroutine IterateMeshRefinementTwoVariables( self, aMesh, field1, flagFn1, tol1, desc1, field2, flagFn2, tol2, desc2, &
	nParticlesStart, nParticlesEnd )
	type(RefineSetup), intent(inout) :: self
	type(PolyMesh2d), intent(inout) :: aMesh
	type(Field), intent(in) :: field1
	procedure(FlagFunction) :: flagFn1
	real(kreal), intent(in) :: tol1
	character(len=*), intent(in) :: desc1
	type(Field), intent(in) :: field2
	procedure(FlagFunction) :: flagFn2
	real(kreal), intent(in) :: tol2
	character(len=*), intent(in) :: desc2
	integer(kint), intent(out) :: nParticlesStart
	integer(kint), intent(out) :: nParticlesEnd
	!
	integer(kint), dimension(2) :: counters
	integer(kint) :: i

	counters = 0
	nParticlesStart = aMesh%particles%N
	do i = 1, aMesh%faces%N
		if ( .NOT. aMesh%faces%hasChildren(i) ) then
			if ( flagFn1(aMesh, field1, tol1, i) ) then
				self%refineFlag(i) = .TRUE.
				counters(1) = counters(1) + 1
			endif
			if ( flagFn2(aMesh, field2, tol2, i ) ) then
				self%refineFlag(i) = .TRUE.
				counters(2) = counters(2) + 1
			endif
		endif
	enddo

	write(logString,'(2A,I8,A)') trim(desc1), ": flagged ", counters(1), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	write(logString,'(2A,I8,A)') trim(desc2), ": flagged ", counters(2), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	call DivideFlaggedFaces( self, aMesh)
	nParticlesEnd = aMesh%particles%N
end subroutine

subroutine IterateMeshRefinementTwoVariablesAndFlowMap( self, aMesh, field1, flagFn1, tol1, desc1, &
														field2, flagFn2, tol2, desc2, flowMapVarTol, &
														nParticlesStart, nParticlesEnd )
	type(RefineSetup), intent(inout) :: self
	type(PolyMesh2d), intent(inout) :: aMesh
	type(Field), intent(in) :: field1
	procedure(FlagFunction) :: flagFn1
	real(kreal), intent(in) :: tol1
	character(len=*), intent(in) :: desc1
	type(Field), intent(in) :: field2
	procedure(FlagFunction) :: flagFn2
	real(kreal), intent(in) :: tol2
	character(len=*), intent(in) :: desc2
	real(kreal), intent(in) :: flowMapVarTol
	integer(kint), intent(out) :: nParticlesStart
	integer(kint), intent(out) :: nParticlesEnd
	!
	integer(kint), dimension(3) :: counters
	integer(kint) :: i

!	call StartSection(log, trim(logKey)//" meshAMR status")

	counters = 0
	nParticlesStart = aMesh%particles%N
	do i = 1, aMesh%faces%N
		if ( .NOT. aMesh%faces%hasChildren(i) ) then
			if ( flagFn1(aMesh, field1, tol1, i) ) then
				self%refineFlag(i) = .TRUE.
				counters(1) = counters(1) + 1
			endif
			if ( flagFn2(aMesh, field2, tol2, i) ) then
				self%refineFlag(i) = .TRUE.
				counters(2) = counters(2) + 1
			endif
			if ( FlowMapVariationRefinement( aMesh, flowMapVarTol, i ) ) then
				self%refineFlag(i) = .TRUE.
				counters(3) = counters(3) + 1
			endif
		endif
	enddo

	write(logString,'(2A,I8,A)') trim(desc1), ": flagged ", counters(1), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	write(logString,'(2A,I8,A)') trim(desc2), ": flagged ", counters(2), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	write(logString,'(A,I8,A)') "flow map refinement : flagged ", counters(3), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	call DivideFlaggedFaces( self, aMesh)
	nParticlesEnd = aMesh%particles%N

!	call EndSection(log)
end subroutine

subroutine IterateMeshRefinementOneVariableAndFlowMap( self, aMesh, aField, flagFn, tol, desc, flowMapVarTol, &
	nParticlesStart, nParticlesEnd )
	type(RefineSetup), intent(inout) :: self
	type(PolyMesh2d), intent(inout) :: aMesh
	type(Field), intent(in) :: aField
	procedure(FlagFunction) :: flagFn
	real(kreal), intent(in) :: tol
	character(len=*), intent(in) :: desc
	real(kreal), intent(in) :: flowMapVarTol
	integer(kint), intent(out) :: nParticlesStart
	integer(kint), intent(out) :: nParticlesEnd
	!
	integer(kint), dimension(2) :: counters
	integer(kint) :: i

	call LogMessage(log, DEBUG_LOGGING_LEVEL, trim(logKey)//" entering ", "IterateMeshRefinementOneVariableAndFlowMap...")
	counters = 0
	call StartSection(log, trim(logKey)//" mesh status")
	call LogStats(aMesh, log)
	nParticlesStart = aMesh%particles%N
	do i = 1, aMesh%faces%N
		if ( .NOT. aMesh%faces%hasChildren(i) ) then
			!call LogMessage(log, DEBUG_LOGGING_LEVEL, trim(logKey)//" crit 1 face i = ", i)
			if ( flagFn( aMesh, aField, tol, i ) ) then
				self%refineFlag(i) = .TRUE.
				counters(1) = counters(1) + 1
			endif
			!call LogMessage(log, DEBUG_LOGGING_LEVEL, trim(logKey)//" crit 2 face i = ", i)
			if ( FlowMapVariationRefinement( aMesh, flowMapVarTol, i) ) then
				self%refineFlag(i) = .TRUE.
				counters(2) = counters(2) + 1
			endif
		endif
	enddo

	write(logString,'(2A,I8,A)') trim(desc), ": flagged ", counters(1), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	write(logString,'(A,I8,A)') "flow map refinement : flagged ", counters(2), " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	call DivideFlaggedFaces( self, aMesh)
	nParticlesEnd = aMesh%particles%N
	call EndSection(log)
end subroutine

subroutine IterateMeshRefinementForFlowMap( self, aMesh, flowMapVarTol, nParticlesStart, nParticlesEnd)
	type(RefineSetup), intent(inout) :: self
	type(PolyMesh2d), intent(inout) :: aMesh
	real(kreal), intent(in) :: flowMapVarTol
	integer(kint), intent(out) :: nParticlesStart
	integer(kint), intent(out) :: nParticlesEnd
	!
	integer(kint) :: counter
	integer(kint) :: i

	counter = 0
	nParticlesStart = aMesh%particles%N
	do i = 1, aMesh%faces%N
		if ( .NOT. aMesh%faces%hasChildren(i) ) then
			if ( FlowMapVariationRefinement( aMesh, flowMapVarTol, i ) ) then
				self%refineFlag(i) = .TRUE.
				counter = counter + 1
			endif
		endif
	enddo

	write(logString,'(A,I8,A)') "flow map refinement : flagged ", counter, " faces."
	call LogMessage(log, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logString)

	call DivideFlaggedFaces(self, aMesh)
	nParticlesEnd = aMesh%particles%N
end subroutine

function ScalarIntegralRefinement( mesh, dataField, tol, faceIndex )
	logical(klog) :: ScalarIntegralRefinement
	type(PolyMesh2d), intent(in) :: mesh
	type(Field), intent(in) :: dataField
	real(kreal), intent(in) :: tol
	integer(kint), intent(in) :: faceIndex
	!
	integer(kint) :: pIndex
	pIndex = mesh%faces%centerParticle(faceIndex)
	ScalarIntegralRefinement = ( abs(dataField%scalar(pIndex)) * mesh%particles%area(pIndex) > tol )
end function

function ScalarVariationRefinement( mesh, dataField, tol, faceIndex )
	logical(klog) :: ScalarVariationRefinement
	type(PolyMesh2d), intent(in) :: mesh
	type(Field), intent(in) :: dataField
	real(kreal), intent(in) :: tol
	integer(kint), intent(in) :: faceIndex
	!
	integer(kint) :: pIndex, i
	real(kreal) :: maxScalar, minScalar
	type(STDIntVector) :: faceVerts

	pIndex = mesh%faces%centerParticle(faceIndex)
	call CCWVerticesAroundFace( mesh, faceVerts, faceindex)
	maxScalar = dataField%scalar(pIndex)
	minScalar = maxScalar
	do i = 1, faceVerts%N
		if ( dataField%scalar( faceVerts%int(i) ) > maxScalar ) maxScalar = dataField%scalar( faceVerts%int(i))
		if ( dataField%scalar( faceVerts%int(i) ) < minScalar ) minScalar = dataField%scalar( faceVerts%int(i))
	enddo
	ScalarVariationRefinement = ( maxScalar - minScalar > tol )
end function

function FlowMapVariationRefinement( mesh, tol, faceIndex )
	logical(klog) :: FlowMapVariationRefinement
	type(PolyMesh2d), intent(in) :: mesh
	real(kreal), intent(in) :: tol
	integer(kint), intent(in) :: faceIndex
	!
	integer(kint) :: pIndex, i
	real(kreal), dimension(3) :: maxX0, minX0, x0i
	type(STDIntVector) :: faceVerts

	pIndex = mesh%faces%centerParticle(faceIndex)
	!call CCWVerticesAroundFace(mesh, faceVerts, faceIndex)
	maxX0 = LagCoord(mesh%particles, pIndex)
	minX0 = LagCoord(mesh%particles, pIndex)
	do i = 1, mesh%faceKind
		!x0i = LagCoord(mesh%particles, faceVerts%int(i) )
		x0i = LagCoord(mesh%particles, mesh%faces%vertices(i,faceIndex))
		if ( x0i(1) < minX0(1) ) minX0(1) = x0i(1)
		if ( x0i(2) < minX0(2) ) minX0(2) = x0i(2)
		if ( x0i(3) < minX0(3) ) minX0(3) = x0i(3)
		if ( x0i(1) > maxX0(1) ) maxX0(1) = x0i(1)
		if ( x0i(2) > maxX0(2) ) maxX0(2) = x0i(2)
		if ( x0i(3) > maxX0(3) ) maxX0(3) = x0i(3)
	enddo
	FlowMapVariationRefinement = ( sum(maxX0 - minX0) > tol )
end function
!
!----------------
! private methods
!----------------
!
subroutine InitLogger(aLog,rank)
! Initialize a logger for this module and processor
	type(Logger), intent(out) :: aLog
	integer(kint), intent(in) :: rank
	write(logKey,'(A,A,I0.2,A)') trim(logKey),'_',rank,' : '
	if ( rank == 0 ) then
		call New(aLog,logLevel)
	else
		call New(aLog,ERROR_LOGGING_LEVEL)
	endif
	logInit = .TRUE.
end subroutine

end module

max_stars_repo_name: pbosler/LPPM

max_issues_repo_name: pbosler/LPPM

max_forks_repo_name: pbosler/LPPM


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #############################################################
c     ##                                                         ##
c     ##  subroutine prtint  --  output of internal coordinates  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "prtint" writes out a set of Z-matrix internal
c     coordinates to an external disk file
c
c
      subroutine prtint (izmt)
      implicit none
      include 'sizes.i'
      include 'atmtyp.i'
      include 'atoms.i'
      include 'files.i'
      include 'inform.i'
      include 'titles.i'
      include 'zclose.i'
      include 'zcoord.i'
      integer i,k,izmt
      logical opened
      character*2 atmc
      character*5 bndc,angc
      character*43 fstr
      character*120 zmtfile
c
c
c     open the output unit if not already done
c
      inquire (unit=izmt,opened=opened)
      if (.not. opened) then
         zmtfile = filename(1:leng)//'.int'
         call version (zmtfile,'new')
         open (unit=izmt,file=zmtfile,status='new')
      end if
c
c     check for large systems needing extended formatting
c
      atmc = 'i6'
      if (n .ge. 100000)  atmc = 'i7'
      if (n .ge. 1000000)  atmc = 'i8'
      if (digits .le. 6) then
         bndc = 'f10.5'
         angc = 'f10.4'
      else if (digits .le. 8) then
         bndc = 'f12.7'
         angc = 'f12.6'
      else
         bndc = 'f14.9'
         angc = 'f14.8'
      end if
c
c     write out the number of atoms and the title
c
      if (ltitle .eq. 0) then
         fstr = '('//atmc//')'
         write (izmt,fstr(1:4))  n
      else
         fstr = '('//atmc//',2x,a)'
         write (izmt,fstr(1:9))  n,title(1:ltitle)
      end if
c
c     output of first three atoms is handled separately
c
      fstr = '('//atmc//',2x,a3,i6,'//atmc//','//bndc//','//atmc//
     &          ','//angc//','//atmc//','//angc//','//'i6)'
      if (n .ge. 1)
     &   write (izmt,fstr)  1,name(1),type(1)
      if (n .ge. 2)
     &   write (izmt,fstr)  2,name(2),type(2),iz(1,2),zbond(2)
      if (n .ge. 3)
     &   write (izmt,fstr)  3,name(3),type(3),iz(1,3),zbond(3),
     &                      iz(2,3),zang(3)
c
c     convert torsional angles to lie in standard range
c
      do i = 4, n
         if (iz(4,i) .eq. 0) then
            do while (ztors(i) .lt. -180.0d0)
               ztors(i) = ztors(i) + 360.0d0
            end do
            do while (ztors(i) .gt. 180.0d0)
               ztors(i) = ztors(i) - 360.0d0
            end do
         end if
      end do
c
c     output the fourth through final atoms
c
      do i = 4, n
         write (izmt,fstr)  i,name(i),type(i),iz(1,i),zbond(i),
     &                      iz(2,i),zang(i),iz(3,i),ztors(i),iz(4,i)
      end do
c
c     addition and deletion of bonds as required
c
      if (nadd.ne.0 .or. ndel.ne.0) then
         fstr = '(2'//atmc//')'
         write (izmt,'()')
         do i = 1, nadd
            write (izmt,fstr(1:5))  (iadd(k,i),k=1,2)
         end do
         if (ndel .ne. 0)  write (izmt,'()')
         do i = 1, ndel
            write (izmt,fstr(1:5))  (idel(k,i),k=1,2)
         end do
      end if
c
c     close the output unit if opened by this routine
c
      if (.not. opened)  close (unit=izmt)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
      integer Pbb,Nbmax
      parameter(Pbb =6)
      parameter(Nbmax =50)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
module setup_init_mod

use constants_mod

real(dp), save :: nodes(N_NODES),ws(N_NODES)

contains

subroutine hook_init_810

    use share_mod, only: mesh

    use mesh_utils_mod, only: mesh_iterate_quads
    use mesh_utils_mod, only: mesh_refine

    use interpol_mod, only: LegendreGaussNodesAndWeights

    integer :: level

    call LegendreGaussNodesAndWeights(N_NODES-1,nodes,ws)
    nodes = 0.5_dp*nodes
    ws = 0.5_dp*ws

    do level = mesh%minlevel,mesh%maxlevel-1
        call mesh_refine(mesh,setup_init_probe_refine)
    end do
    call mesh_iterate_quads(mesh,fill)

end subroutine

subroutine fill(quad)

    use setup_config_mod
    use equations_const_mod

    use mesh_types_mod, only: quad_t
    use quad_utils_mod, only: quad_get_coords

    use setup_manufacsol_mod

    use share_mod, only: rt => runtime

    type(quad_t), intent(inout) :: quad

    real(dp) :: dens(N_NODES,N_NODES,N_NODES)
    real(dp) :: velx(N_NODES,N_NODES,N_NODES)
    real(dp) :: vely(N_NODES,N_NODES,N_NODES)
    real(dp) :: velz(N_NODES,N_NODES,N_NODES)
    real(dp) :: momx(N_NODES,N_NODES,N_NODES)
    real(dp) :: momy(N_NODES,N_NODES,N_NODES)
    real(dp) :: momz(N_NODES,N_NODES,N_NODES)
    real(dp) :: pres(N_NODES,N_NODES,N_NODES)
    real(dp) :: ener(N_NODES,N_NODES,N_NODES)

    real(dp) :: coords(N_NODES,N_NODES,N_NODES,N_DIMS)
    real(dp) :: radius(N_NODES,N_NODES,N_NODES)

    real(dp) :: u(N_NODES,N_NODES,N_NODES,N_VARS)
    real(dp) :: r(N_NODES,N_NODES,N_NODES,N_VARS)

    real(dp) :: xs(N_NODES)
    real(dp) :: ys(N_NODES)
    real(dp) :: zs(N_NODES)

    real(dp) :: time

    integer :: i,j,k, ii,jj,kk

    coords = quad_get_coords(quad)

    u = 0.0_dp
    r = 0.0_dp

    time = rt%simtime
    !time = 0.0_dp
    ! time = 0.345_dp

    do k = 1,N_NODES; do j = 1,N_NODES; do i = 1,N_NODES;
        xs = coords(i,j,k,1) + quad%delta(1)/REAL(N_NODES,dp)*nodes
        ys = coords(i,j,k,2) + quad%delta(2)/REAL(N_NODES,dp)*nodes
        zs = coords(i,j,k,3) + quad%delta(3)/REAL(N_NODES,dp)*nodes

        do kk = 1,N_NODES; do jj = 1,N_NODES; do ii = 1,N_NODES;
            u(i,j,k,:) = u(i,j,k,:) + ws(ii)*ws(jj)*ws(kk)*init_manufacsol(time,xs(ii),ys(jj),zs(kk))
            r(i,j,k,:) = r(i,j,k,:) + ws(ii)*ws(jj)*ws(kk)*calc_manufacsol(time,xs(ii),ys(jj),zs(kk))
        end do; end do; end do
    end do; end do; end do

    quad%pld%hydro%state(:,:,:,DENS_VAR) = u(:,:,:,DENS_VAR)
    quad%pld%hydro%state(:,:,:,MOMX_VAR) = u(:,:,:,MOMX_VAR)
    quad%pld%hydro%state(:,:,:,MOMY_VAR) = u(:,:,:,MOMY_VAR)
    quad%pld%hydro%state(:,:,:,MOMZ_VAR) = u(:,:,:,MOMZ_VAR)
    quad%pld%hydro%state(:,:,:,ENER_VAR) = u(:,:,:,ENER_VAR)

    quad%aux%manufacsol%u0 = quad%pld%hydro%state
    quad%aux%manufacsol%st = r

end subroutine

function setup_init_probe_refine(quad) result(ok)

    use setup_config_mod
    use mesh_types_mod, only: quad_t
    use mesh_const_mod, only: MESH_CORNERS
    use quad_utils_mod, only: quad_get_center
    use quad_utils_mod, only: quad_get_vertices

    type(quad_t), intent(in)        :: quad
    logical                         :: ok

    real(dp) :: verts(N_DIMS,MESH_CORNERS)
    real(dp) :: center(N_DIMS)

    verts = quad_get_vertices(quad)

    center = quad_get_center(quad)

    ok = .false.

    ! return
    !if (rand() < 0.2) then
    if (any(abs(SQRT(SUM((verts)**2,dim=1))) < 0.5 * 0.5*diameter)) then
    ! if (center(1) < 0) then
    !if (all(center < 0)) then
        ok = .true.
    end if

end function

function setup_init_probe_coarsen(quads) result(ok)

    use mesh_types_mod, only: quad_t
    use mesh_const_mod, only: MESH_CHILDREN
    use setup_config_mod, only: setup_amr_threshold_coarsen

    type(quad_t), intent(in)    :: quads(MESH_CHILDREN)
    logical                     :: ok

    real(dp)                    :: sensor(MESH_CHILDREN)

    integer :: i

    ok = .false.

end function

end module

max_stars_repo_name: jmark/nemo

max_issues_repo_name: jmark/nemo

max_forks_repo_name: jmark/nemo