text:
!> Module for PC-SAFT pure-component parameters and binary interaction
!> parameters. Also contains parameters for the PeTS equation of state.

!> NB: If you want to add new parameters here, beware that different authors are
!> inconsistent wrt. whether beta should be multiplied with pi/6 or not. For
!> this reason you should validate your results before using these parameters.
!> If you get strange results, try multiplying beta with pi/6=0.5236.
module pc_saft_parameters
  use AssocSchemeUtils
  use eosdata, only: eosPC_SAFT, eosPeTS
  use thermopack_constants, only: Rgas => Rgas_default, uid_len, ref_len
  implicit none
  save

  real, parameter :: PI_6 = 4.0*atan(1.0)/6.0

  !> PURE COMPONENT PARAMETERS.
  !> This data structure stores pure component parameters for the PC-SAFT
  !> equation of state.
  ! ---------------------------------------------------------------------------
  type :: pc_saft_data
     integer :: eosidx
     character (len=uid_len) :: compName
     ! Pure component fitted parameters.
     real :: m        !< [-]. Mean number of segments.
     real :: sigma    !< [m]. Temperature-independent segment diameter.
     real :: eps_depth_divk !< [K]. Well depth divided by Boltzmann's constant k.
     ! Association parameters.
     real :: eps  !< [J/mol].
     real :: beta !< [-]. Also known as kappa in SAFT literature.
     ! Association scheme.
     integer :: assoc_scheme
     ! Set number
     character (len=ref_len) :: ref
  end type pc_saft_data

  ! SOURCES:
  ! [0] Regressed by in-house.
  ! [1] Kontogeorgis & Folas, "Thermodynamic Models for Industrial Applications", Wiley 2010 (appendix A).
  ! [2] Tang & Gross. (2010) "Modeling the phase equilibria of hydrogen sulfide
  !     and carbon dioxide in mixtures with hydrocarbons an water using the
  !     PCP-SAFT equation of state". Fluid Phase Equilibria 293.
  ! [3] Smith (2017) MSc thesis. "Measurement and modelling of the vapour--liquid equilibria of binary mixtures of water and alkanols"
  type(pc_saft_data), parameter :: PCcx1 = pc_saft_data(eosPC_SAFT,"C1", &
       1.0000, 3.7039E-10, 150.03, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx2 = pc_saft_data(eosPC_SAFT,"C2", &
       1.6069, 3.5206E-10, 191.42, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx3 = pc_saft_data(eosPC_SAFT,"C3", &
       2.0020, 3.6184E-10, 208.11, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx4 = pc_saft_data(eosPC_SAFT,"NC4", &
       2.3316, 3.7086E-10, 222.88, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx5 = pc_saft_data(eosPC_SAFT,"NC5", &
       2.6896, 3.7729E-10, 231.2, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx6 = pc_saft_data(eosPC_SAFT,"NC6", &
       3.0576, 3.7983E-10, 236.77, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx7 = pc_saft_data(eosPC_SAFT,"NC7", &
       3.4831, 3.8049E-10, 238.4, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx8 = pc_saft_data(eosPC_SAFT,"NC8", &
       3.8176, 3.8373E-10, 242.78, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx9 = pc_saft_data(eosPC_SAFT,"NC9", &
       4.2079, 3.8448E-10, 244.51, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx10 = pc_saft_data(eosPC_SAFT,"NC10", &
       4.6627, 3.8384E-10, 243.87, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx11 = pc_saft_data(eosPC_SAFT,"BENZENE", &
       2.4653, 3.6478E-10, 287.35, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx12 = pc_saft_data(eosPC_SAFT,"TOLU", &
       2.8149, 3.7169E-10, 285.69, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx13 = pc_saft_data(eosPC_SAFT,"N2", &
       1.2053, 3.313E-10,  90.96,  0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx14 = pc_saft_data(eosPC_SAFT,"O2", &
       1.1217, 3.2098E-10, 114.96,  0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx15 = pc_saft_data(eosPC_SAFT,"CYCLOHEX", &
       2.5303, 3.8499E-10, 278.11, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]

  type(pc_saft_data), parameter :: PCcx16 = pc_saft_data(eosPC_SAFT,"H2S", &
       1.6941, 3.0214E-10, 226.79, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1] H2S without association.
  type(pc_saft_data), parameter :: PCcx17 = pc_saft_data(eosPC_SAFT,"H2S", &
       1.355, 3.309E-10, 234.25, 780.8*Rgas, 0.001, assoc_scheme_2B,"Default/Tang_Gross2010") ! [2]

  ! Grenner et al. 2006. doi 10.1021/ie0605332. DeltaP 0.93, DeltaRho=2.62 for T>323K.
  type(pc_saft_data), parameter :: PCcx18 = pc_saft_data(eosPC_SAFT,"H2O", &
       1.50, 2.6273E-10, 180.3, 1804.22*Rgas, 0.0942, assoc_scheme_4C,"Default/Grenner2006")

  type(pc_saft_data), parameter :: PCcx19 = pc_saft_data(eosPC_SAFT,"H2O", &
       1.0656, 3.0007E-10, 366.51, 2500.7*Rgas, 0.034868,  assoc_scheme_2B,"Kontogeorgis_Folas2001") ! [1]

  type(pc_saft_data), parameter :: PCcx20 = pc_saft_data(eosPC_SAFT,"CO2", &
       2.0729, 2.7852E-10, 169.21, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  ! [10]
  type(pc_saft_data), parameter :: PCcx21 = pc_saft_data(eosPC_SAFT,"CO2", &
       2.58239e+00 ,2.56106e-10 ,1.51691e+02 ,0.0 ,0.0,no_assoc,"Tang_Gross2010")

  ! [0] CO2 3B regressed to Tc, Pc and Psat. AAD_P 0.22%, AAD_liqv 3.08%
  type(pc_saft_data), parameter :: PCcx22 = pc_saft_data(eosPC_SAFT,"CO2", &
       1.72897e+00 ,2.76102e-10 ,8.37215e+01 ,9.25650e+03 ,1.71605e-01, assoc_scheme_3B,"Smith2017")

  ! [1]
  type(pc_saft_data), parameter :: PCcx23 = pc_saft_data(eosPC_SAFT,"AR", &
       0.9285, 3.4784E-10, 122.23, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")

  type(pc_saft_data), parameter :: PCcx24 = pc_saft_data(eosPC_SAFT, "NC11", &
       4.8640, 3.8986E-10, 246.9390, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx25 = pc_saft_data(eosPC_SAFT, "NC12", &
       5.2133, 3.9115E-10, 248.0042, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx26 = pc_saft_data(eosPC_SAFT, "NC13", &
       5.5627, 3.9227E-10, 248.9356, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx27 = pc_saft_data(eosPC_SAFT, "NC14", &
       5.9120, 3.9326E-10, 249.7570, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx28 = pc_saft_data(eosPC_SAFT, "NC15", &
       6.2614, 3.9412E-10, 250.4867, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx29 = pc_saft_data(eosPC_SAFT, "NC16", &
       6.6107, 3.9490E-10, 251.1392, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx30 = pc_saft_data(eosPC_SAFT, "NC17", &
       6.9600, 3.9559E-10, 251.7263, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx31 = pc_saft_data(eosPC_SAFT, "NC18", &
       7.3094, 3.9622E-10, 252.2573, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx32 = pc_saft_data(eosPC_SAFT, "NC19", &
       7.6587, 3.9678E-10, 252.7398, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx33 = pc_saft_data(eosPC_SAFT, "NC20", &
       8.0081, 3.9730E-10, 253.1802, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx34 = pc_saft_data(eosPC_SAFT, "NC21", &
       8.3574, 3.9777E-10, 253.5838, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx35 = pc_saft_data(eosPC_SAFT, "NC22", &
       8.7068, 3.9820E-10, 253.9550, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx36 = pc_saft_data(eosPC_SAFT, "NC23", &
       9.0561, 3.9860E-10, 254.2975, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx37 = pc_saft_data(eosPC_SAFT, "NC24", &
       9.4055, 3.9897E-10, 254.6147, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")
  type(pc_saft_data), parameter :: PCcx38 = pc_saft_data(eosPC_SAFT, "NC25", &
       9.7548, 3.9931E-10, 254.9091, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001")

  type(pc_saft_data), parameter :: PCcx39 = pc_saft_data(eosPC_SAFT,"NH3", & ! [0] ammonia (AAD(P,liqv)=(2.0,0.61))
       2.25807 , 2.37802e-10 ,1.26868e+02 , 1.10326e+04, 2.04790e-01, assoc_scheme_3B,"Default/InHouse")

  type(pc_saft_data), parameter :: PCcx40 = pc_saft_data(eosPC_SAFT,"IC4", &
       2.2616, 3.7574E-10, 216.53, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]
  type(pc_saft_data), parameter :: PCcx41 = pc_saft_data(eosPC_SAFT,"IC5", &
       2.562, 3.8296E-10, 230.75, 0.0, 0.0, no_assoc,"Default/Kontogeorgis_Folas2001") ! [1]

  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.52,0.37)
  type(pc_saft_data), parameter :: PCcx42 = pc_saft_data(eosPC_SAFT,"MEOH", &
       2.1, 2.7998E-10, 197.23, 2535.0*Rgas, 0.8230*PI_6, assoc_scheme_2C,"Default/deVilliers2011")

  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.18,0.21)
  type(pc_saft_data), parameter :: PCcx43 = pc_saft_data(eosPC_SAFT,"ETOH", &
       2.3609, 3.1895E-10, 207.56, 2695.69*Rgas, 0.03270*PI_6, assoc_scheme_2C,"Default/deVilliers2011")

  type(pc_saft_data), parameter :: PCcx44 = pc_saft_data(eosPC_SAFT,"ETOH", &
       1.2309, 4.1057E-10, 316.91, 2811.02*Rgas, 0.00633*PI_6, assoc_scheme_2B,"Tang_Gross2010")

  ! Gross and Sadowski 2002, Ind. Eng. Chem. Res. 41 (DeltaP,Deltav)=(0.58,0.71)
  type(pc_saft_data), parameter :: PCcx45 = pc_saft_data(eosPC_SAFT,"ETOH", &
       2.3827, 3.1771E-10, 198.24, 2653.4*Rgas, 0.032384, assoc_scheme_2B,"Gross_Sadowski2002")

  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.27,0.25)
  type(pc_saft_data), parameter :: PCcx46 = pc_saft_data(eosPC_SAFT,"PROP1OL", &
       2.9537, 3.2473E-10, 226.36, 2448.02*Rgas, 0.02280*PI_6, assoc_scheme_2C,"Default/deVilliers2011")
  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.19, 0.33)
  type(pc_saft_data), parameter :: PCcx47 = pc_saft_data(eosPC_SAFT,"BUT1OL", &
       2.9614, 3.5065E-10, 253.29, 2601.00*Rgas, 0.00740*PI_6, assoc_scheme_2C,"Default/deVilliers2011")
  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.09,0.10)
  type(pc_saft_data), parameter :: PCcx48 = pc_saft_data(eosPC_SAFT,"PENT1OL", &
       3.1488, 3.6350E-10, 261.96, 2555.34*Rgas, 0.00750*PI_6, assoc_scheme_2C,"Default/deVilliers2011")
  ! de Villiers 2011. doi 10.1021/ie200521k. 0.5<Tr<0.9, (DeltaP,DeltaRho)=(0.09,0.42)
  type(pc_saft_data), parameter :: PCcx49 = pc_saft_data(eosPC_SAFT,"HEX1OL", &
       3.1542, 3.8188E-10, 277.24, 2946.10*Rgas, 0.00250*PI_6, assoc_scheme_2C,"Default/deVilliers2011")

  ! Gross and Sadowski 2002, Ind. Eng. Chem. Res. 41 (DeltaP,Deltav)=(0.85,1.71)
  type(pc_saft_data), parameter :: PCcx50 = pc_saft_data(eosPC_SAFT,"PROP1OL", &
       2.9997, 3.2522E-10, 233.40, 2276.8*Rgas, 0.015268, assoc_scheme_2B,"Tang_Gross2010")
  ! Gross and Sadowski 2002, Ind. Eng. Chem. Res. 41 (DeltaP,Deltav)=(1.78,1.63)
  type(pc_saft_data), parameter :: PCcx51 = pc_saft_data(eosPC_SAFT,"BUT1OL", &
       2.7515, 3.6189E-10, 259.59, 2544.6*Rgas, 0.006692, assoc_scheme_2B,"Tang_Gross2010")
  ! Gross and Sadowski 2002, Ind. Eng. Chem. Res. 41 (DeltaP,Deltav)=(0.5,0.52)
  type(pc_saft_data), parameter :: PCcx52 = pc_saft_data(eosPC_SAFT,"PENT1OL", &
       3.6260, 3.4508E-10, 247.28, 2252.1*Rgas, 0.010319, assoc_scheme_2B,"Tang_Gross2010")
  ! Gross and Sadowski 2002, Ind. Eng. Chem. Res. 41 (DeltaP,Deltav)=(0.58,0.71)
  type(pc_saft_data), parameter :: PCcx53 = pc_saft_data(eosPC_SAFT,"HEX1OL", &
       3.5146, 3.6735E-10, 262.32, 2538.9*Rgas, 0.005747, assoc_scheme_2B,"Tang_Gross2010")

  type(pc_saft_data), parameter :: PCcx54 = pc_saft_data(eosPeTS,"AR", &
       1.000, 3.42E-10, 136.0, 0.0, 0.0, no_assoc,"Default")

  integer, parameter :: nPCmodels = 54
  type(pc_saft_data), dimension(nPCmodels), parameter :: PCarray = (/ &
       PCcx1, &
       PCcx2, &
       PCcx3, &
       PCcx4, &
       PCcx5, &
       PCcx6, &
       PCcx7, &
       PCcx8, &
       PCcx9, &
       PCcx10, &
       PCcx11, &
       PCcx12, &
       PCcx13, &
       PCcx14, &
       PCcx15, &
       PCcx16, &
       PCcx17, &
       PCcx18, &
       PCcx19, &
       PCcx20, &
       PCcx21, &
       PCcx22, &
       PCcx23, &
       PCcx24, &
       PCcx25, &
       PCcx26, &
       PCcx27, &
       PCcx28, &
       PCcx29, &
       PCcx30, &
       PCcx31, &
       PCcx32, &
       PCcx33, &
       PCcx34, &
       PCcx35, &
       PCcx36, &
       PCcx37, &
       PCcx38, &
       PCcx39, &
       PCcx40, &
       PCcx41, &
       PCcx42, &
       PCcx43, &
       PCcx44, &
       PCcx45, &
       PCcx46, &
       PCcx47, &
       PCcx48, &
       PCcx49, &
       PCcx50, &
       PCcx51, &
       PCcx52, &
       PCcx53, &
       PCcx54/)

  !> TEMPERATURE-INDEPENDENT INTERACTION PARAMETERS FOR PC-SAFT DISPERSION TERM.
  ! ----------------------------------------------------------------------------
  type :: PCkijdata
     integer:: eosidx
     character (len=uid_len) :: uid1, uid2
     character(len=ref_len) :: ref
     real :: kijvalue
  end type PCkijdata

  ! All interaction parameters from the following reference:
  ! [2] Tang & Gross. (2010) "Modeling the phase equilibria of hydrogen sulfide
  !     and carbon dioxide in mixtures with hydrocarbons and water using the
  !     PCP-SAFT equation of state". Fluid Phase Equilibria 293.
  integer, parameter :: PCmaxkij = 31
  type(PCkijdata), dimension(PCmaxkij), parameter :: PCkijdb = (/ &
       PCkijdata(eosPC_SAFT,"H2S","C1","Default/Tang_Gross2010",0.0425), &
       PCkijdata(eosPC_SAFT,"H2S","C2","Default/Tang_Gross2010",0.072), &
       PCkijdata(eosPC_SAFT,"H2S","C3","Default/Tang_Gross2010",0.069), &
       PCkijdata(eosPC_SAFT,"H2S","NC4","Default/Tang_Gross2010",0.067), &
       PCkijdata(eosPC_SAFT,"H2S","NC5","Default/Tang_Gross2010",0.073), &
       PCkijdata(eosPC_SAFT,"H2S","NC6","Default/Tang_Gross2010",0.073), &
       PCkijdata(eosPC_SAFT,"H2S","NC7","Default/Tang_Gross2010",0.078), &
       PCkijdata(eosPC_SAFT,"H2S","NC9","Default/Tang_Gross2010",0.086), &
       PCkijdata(eosPC_SAFT,"H2S","NC10","Default/Tang_Gross2010",0.077), &
       PCkijdata(eosPC_SAFT,"H2S","IC4","Default/Tang_Gross2010",0.060), &
       PCkijdata(eosPC_SAFT,"H2S","IC5","Default/Tang_Gross2010",0.076), &
       PCkijdata(eosPC_SAFT,"H2S","CYCLOHEX","Default/Tang_Gross2010",0.082), &
       PCkijdata(eosPC_SAFT,"H2S","NC8","Default/Tang_Gross2010",0.0), &
       PCkijdata(eosPC_SAFT,"CO2","H2S","Default/Tang_Gross2010",0.0223), &
       PCkijdata(eosPC_SAFT,"CO2","BENZENE","Default/Tang_Gross2010",0.025), &
       PCkijdata(eosPC_SAFT,"CO2","TOLU","Default/Tang_Gross2010",0.026), &
       PCkijdata(eosPC_SAFT,"CO2","C1","Default/Gross_Sadowski2001",0.0650), & ! Value from Gross and Sadowski 2001.
       PCkijdata(eosPC_SAFT,"CO2","C2","Default/Tang_Gross2010",0.102), &
       PCkijdata(eosPC_SAFT,"CO2","C3","Default/Tang_Gross2010",0.0107), &
       PCkijdata(eosPC_SAFT,"CO2","NC4","Default/Tang_Gross2010",0.109), &
       PCkijdata(eosPC_SAFT,"CO2","NC5","Default/Tang_Gross2010",0.120), &
       PCkijdata(eosPC_SAFT,"CO2","NC6","Default/Tang_Gross2010",0.123), &
       PCkijdata(eosPC_SAFT,"CO2","NC7","Default/Tang_Gross2010",0.115), &
       PCkijdata(eosPC_SAFT,"CO2","NC8","Default/Tang_Gross2010",0.132), &
       PCkijdata(eosPC_SAFT,"CO2","NC9","Default/Tang_Gross2010",0.122), &
       PCkijdata(eosPC_SAFT,"CO2","NC10","Default/Tang_Gross2010",0.133), &
       PCkijdata(eosPC_SAFT,"CO2","IC4","Default/Tang_Gross2010",0.112), &
       PCkijdata(eosPC_SAFT,"CO2","IC5","Default/Tang_Gross2010",0.116), &
       PCkijdata(eosPC_SAFT,"CO2","CYCLOHEX","Default/Tang_Gross2010",0.125), &
       PCkijdata(eosPC_SAFT,"CO2","BENZENE","Default/Tang_Gross2010",0.087), &
       PCkijdata(eosPC_SAFT,"CO2","TOLU","Default/Tang_Gross2010",0.108) &
       /)

contains

  !> Checks that we use the correct gas constant when initializing the
  !> pc_saft_data parameters.
  logical function Rgas_is_correct()
    use thermopack_constants, only: Rgas_default
    Rgas_is_correct = (Rgas == Rgas_default)
  end function Rgas_is_correct

  !> Get the index in the PCarray of the component having uid given by
  !> compName. idx=0 if component isn't in database.
  function getPCdataIdx(eosidx,compName,param_ref) result(idx)
    use stringmod, only: str_eq, string_match
    integer, intent(in) :: eosidx
    character(len=*), intent(in) :: compName, param_ref
    integer :: idx, idx_default
    logical :: found

    if (.not. Rgas_is_correct()) call stoperror("Rgas_default must be default Rgas parameter.")

    found = .false.
    idx = 1
    idx_default = -1
    do while (idx <= nPCmodels)
      if ((eosidx==PCarray(idx)%eosidx) .and. &
           str_eq(compName, PCarray(idx)%compName)) then
        if (string_match(param_ref,PCarray(idx)%ref)) then
          found = .true.
          exit
        elseif (string_match("DEFAULT", PCarray(idx)%ref)) then
          idx_default = idx
        endif
      endif
      idx = idx + 1
    enddo

    if (.not. found) then
      if (verbose) then
        print *, "No PC-SAFT parameters for compName, ref ", compName, trim(param_ref)
      endif
      if (idx_default > 0) then
        idx = idx_default
        print *, "Using default parameter set for "//trim(compName)
      else
        print *, "ERROR FOR COMPONENT ", compname
        call stoperror("The PC-SAFT parameters don't exist.")
      endif
    end if

  end function getPCdataIdx

  !> Retrieve binary interaction parameter for components uid1 and uid2.
  !> If no kij is stored in the database PCkijdb, it returns 0.0.
  function getPCkij (eosidx,uid1,uid2,param_ref) result(kijvalue)
    use stringmod, only: str_eq, string_match
    integer, intent(in) :: eosidx
    character(len=*), intent(in) :: uid1, uid2, param_ref
    real :: kijvalue
    integer :: idx, idx_default
    logical :: match_11_22, match_12_21, found

    kijvalue = 0.0 ! default value if the binary is not in PCkijdb.
    idx = 1
    idx_default = -1
    found = .false.
    do idx = 1,PCmaxkij
       match_11_22 = str_eq(uid1,PCkijdb(idx)%uid1) .and. str_eq(uid2,PCkijdb(idx)%uid2)
       match_12_21 = str_eq(uid1,PCkijdb(idx)%uid2) .and. str_eq(uid2,PCkijdb(idx)%uid1)

       if ( eosidx==PCkijdb(idx)%eosidx .and. (match_11_22 .or. match_12_21)) then
         if (string_match(param_ref,PCkijdb(idx)%ref)) then
           kijvalue = PCkijdb(idx)%kijvalue
           found = .true.
           exit
         elseif (string_match("DEFAULT", PCkijdb(idx)%ref)) then
           idx_default = idx
         endif
       endif
     end do

     if (.not. found .and. idx_default > 0) then
       kijvalue = PCkijdb(idx_default)%kijvalue
     endif
  end function getPCkij

  subroutine getPcSaftKij_allComps(nc,comp,eosidx,kij)
    use compdata, only: gendata_pointer
    ! Input
    integer, intent(in) :: nc
    type(gendata_pointer), intent(in) :: comp(nc)
    integer, intent(in) :: eosidx
    ! Output
    real, intent(out) :: kij(nc,nc)
    ! Locals
    integer :: ic,jc

    kij = 0.0
    do ic=1,nc
       do jc=ic+1,nc
          kij(ic,jc) = getPCkij(eosidx,comp(ic)%p_comp%ident,comp(jc)%p_comp%ident,"DEFAULT")
          kij(jc,ic) = kij(ic,jc)
       end do
    end do

  end subroutine getPcSaftKij_allComps

  subroutine getPcSaftPureParams_allComps(nc,comp,eosidx,param_ref,found,m,sigma,eps_depth_divk,eps,beta,scheme)
    use compdata, only: gendata_pointer
    ! Input
    integer, intent(in) :: nc
    type(gendata_pointer), intent(in) :: comp(nc)
    integer, intent(in) :: eosidx
    character(len=*), intent(in) :: param_ref
    ! Output
    logical, intent(out) :: found(nc)
    real, intent(out) :: m(nc), sigma(nc), eps_depth_divk(nc), eps(nc), beta(nc)
    integer, intent(out) :: scheme(nc)
    ! Locals
    integer :: ic

    do ic=1,nc
       call getPcSaftpureParams_singleComp(comp(ic)%p_comp%ident,eosidx,param_ref,found(ic),m(ic),&
            sigma(ic),eps_depth_divk(ic),eps(ic),beta(ic),scheme(ic))
    end do

  end subroutine getPcSaftPureParams_allComps

  subroutine getPcSaftpureParams_singleComp(compName,eosidx,param_ref,found,m,sigma,eps_depth_divk,eps,beta,scheme)
    ! Input
    character(len=*), intent(in) :: compName, param_ref
    integer, intent(in) :: eosidx
    ! Output
    logical, intent(out) :: found
    real, intent(out) :: m, sigma, eps_depth_divk, eps, beta
    integer, intent(out) :: scheme
    ! Locals
    integer :: idx

    idx = getPCdataIdx(eosidx,compName,param_ref)
    if ( idx == 0 ) then
       found = .false.
       return
    end if

    found = .true.
    m = PCarray(idx)%m
    scheme = PCarray(idx)%assoc_scheme
    sigma = PCarray(idx)%sigma
    eps_depth_divk = PCarray(idx)%eps_depth_divk
    eps = PCarray(idx)%eps
    beta = PCarray(idx)%beta
    scheme = PCarray(idx)%assoc_scheme
  end subroutine getPcSaftpureParams_singleComp

end module pc_saft_parameters

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


text:
\ color.f       - x4 console colour handling
\ ------------------------------------------------------------------------

  .( color.f )

\ ------------------------------------------------------------------------

  headers>

enum: colours
  := black
  := red
  := green
  := brown
  := blue
  := magenta
  := cyan
  := white
  := gray
  := lt_red
  := lt_green
  := yellow
  := lt_blue
  := lt_magenta
  := lt_cyan
  := lt_white
;enum

\ ------------------------------------------------------------------------
\ character attributes

enum: attrs
  1 /= :standout:
  2 /= :underline:
  4 /= :reverse:
  8 /= :bold:
 16 /= :alt:             \ alt charset (ibm box charset)
;enum

\ ------------------------------------------------------------------------
\ character attributes

enum: attrs
  1 /= :standout:
  2 /= :underline:
  4 /= :reverse:
  8 /= :bold:
 16 /= :alt:             \ alt charset (ibm box charset)
;enum

\ ------------------------------------------------------------------------

  0 var attrib              \ current color atributes + bold/standout etc

\ ------------------------------------------------------------------------

  <headers

: (>attr)       ( n1 mask --- )
  attrib and or
  !> attrib ;

\ ------------------------------------------------------------------------

: (>fg)         ( n1 --- )      $fff0 (>attr) ;
: (>bg)         ( n1 --- ) 4 << $ff0f (>attr) ;

\ ------------------------------------------------------------------------
\ set new foreground or background colors

  headers>

: >fg           ( n1 --- ) dup (>fg) setaf ;
: >bg           ( n1 --- ) dup (>bg) setab ;

\ ------------------------------------------------------------------------

: fg            ( --- n1 ) attrib $f and ;
: bg            ( --- n1 ) attrib $f0 and 4 >> ;

\ ------------------------------------------------------------------------

: >attrib       ( n1 --- )
  $ff00 (>attr)
  fg setaf
  bg setab ;

\ ------------------------------------------------------------------------
\ combine fg and bg colors into a single 8 bit attribute

: >color        ( fg bg --- color ) 4 << or ;

\ ------------------------------------------------------------------------
\ turn various attributes on and off

\ the terminfo for some terminals doesnt list an sgr format string even
\ though the terminals themselves support one.  (rxvt types only show
\ the sgr0 format string - eterm, linux terminal have both sgr and sgr0)
\
\ some of these attributes have two associated escape sequences.  one to
\ turn the attribute on, the other to turn it off.  someone in their
\ infinite moronic stupidity neglected to include a format string to turn
\ bold off. the only way to turn bold off reliably over all terminal types
\ is to turn off 'all' attributes by using an sgr0 and to then
\ re-establish the ones you realy didnt want to lose.
\
\ the book termcap and terminfo rationalizes this somewhat but its a
\ complete CROCK OF SHIT!  terminfo is BRAINDEAD here and the people who
\ put together incomplete terminfo file are ALSO braindead!!!!!

  <headers

: (>pref)       ( c1 --- )
  sgr0                      \ clear all attribs
  dup :standout:  and ?: smso noop
  dup :underline: and ?: smul noop
  dup :reverse:   and ?: rev  noop
      :bold:      and ?: bold noop ;

\ ------------------------------------------------------------------------

\ the above sgr0 very kindly resets the forground and background colours
\ to their defaults.  nice of the terminal to force me into changing them
\ back again, specially when changing attributes is one of the slowest
\ things you can do in a terminal

  headers>

: >pref     ( c1 --- )
  dup 8 << $ff (>attr)
  (>pref)
  fg setaf bg setab ;          \ restore the damned colors

\ ------------------------------------------------------------------------

  <headers

: +pref         ( bit --- ) attrib 8 >> or >pref ;
: -pref         ( bit --- ) not attrib 8 >> and >pref ;

\ ------------------------------------------------------------------------
\ fetch an attribute mask. see if this attribute is set or not

: (aset/aclr)   ( a1 --- mask f1 )
  @ dup 8 <<                \ fetch attribute mask for set/clear
  attrib and ;              \ is this atribute set?

\ ------------------------------------------------------------------------
\ create word to set an attribute

: aset          ( bit --- )  create, does> (aset/aclr) ?: drop +pref ;
: aclr          ( bit --- )  create, does> (aset/aclr) ?: -pref drop ;

\ ------------------------------------------------------------------------
\ these words will not set/clr any attribute that is already set/clr

  headers>

 :standout:  aset >so       :standout:  aclr <so
 :underline: aset >ul       :underline: aclr <ul
 :reverse:   aset >rev      :reverse:   aclr <rev
 :bold:      aset >bold     :bold:      aclr <bold

\ ------------------------------------------------------------------------
\ this does not reset colors (no thanx to terminfo)

: >norm 0 >pref ;

\ ========================================================================

ext: f


module potentialBCs_mumps

use, intrinsic :: iso_fortran_env, only : stderr=>error_unit
use, intrinsic :: ieee_arithmetic, only : ieee_is_finite

use mpimod, only : mpi_cfg
use phys_consts, only: wp, pi, Re, debug
use grid, only: lx1, lx2, lx2all, lx3all, gridflag
use meshobj, only: curvmesh
use interpolation, only : interp1,interp2
use timeutils, only : dateinc, date_filename, find_lastdate
use reader, only : get_grid2, get_simsize2, get_Efield
use config, only: gemini_cfg
use efielddataobj, only: efielddata

implicit none (type, external)
private
public :: potentialbcs2D, potentialbcs2D_fileinput, init_Efieldinput, &
            compute_rootBGEfields

type(efielddata) :: efield

integer, private :: ix1ref,ix2ref,ix3ref    ! reference locations for field line mapping
integer, private :: ix1eq=-1                ! index for the equatorial location in terms of index into the x%x1 array; used by default boundary conditions

contains

subroutine init_Efieldinput(dt,t,cfg,ymd,UTsec,x)
  !> Initialize variables to hold electric field input file data, can be called by any worker but only root does anything
  real(wp), intent(in) :: dt,t
  type(gemini_cfg), intent(in) :: cfg
  integer, dimension(3), intent(in) :: ymd
  real(wp), intent(in) :: UTsec
  class(curvmesh), intent(in) :: x

  !> initializes the auroral electric field/current and particle inputs to read in a file corresponding to the first time step
  if (mpi_cfg%myid==0 .and. cfg%flagE0file==1) then    !only root needs these...
    call efield%init(cfg,cfg%E0dir,x,dt,cfg%dtE0,ymd,UTsec)
  end if
end subroutine init_Efieldinput

subroutine potentialBCs2D_fileinput(dtmodel,t,ymd,UTsec,cfg,x,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3, &
                                  Vmaxx3,E01all,E02all,E03all,flagdirich)
  !! A FILE INPUT BASED BOUNDARY CONDITIONS FOR ELECTRIC POTENTIAL OR
  !! FIELD-ALIGNED CURRENT.
  !! NOTE: THIS IS ONLY CALLED BY THE ROOT PROCESS
  real(wp), intent(in) :: dtmodel
  real(wp), intent(in) :: t
  integer, dimension(3), intent(in) :: ymd    !date for which we wish to calculate perturbations
  real(wp), intent(in) :: UTsec
  type(gemini_cfg), intent(in) :: cfg
  class(curvmesh), intent(in) :: x

  real(wp), dimension(:,:), intent(inout), target :: Vminx1,Vmaxx1
  !! intent(out)
  real(wp), dimension(:,:), intent(inout) :: Vminx2,Vmaxx2
  !! intent(out)
  real(wp), dimension(:,:), intent(inout) :: Vminx3,Vmaxx3
  !! intent(out)
  real(wp), dimension(:,:,:), intent(inout) :: E01all,E02all,E03all
  !! intent(out)
  integer, intent(out) :: flagdirich
  integer :: ix1,ix2,ix3

  !> COMPUTE SOURCE/FORCING TERMS FROM BACKGROUND FIELDS, ETC.
  E01all = 0
  !! do not allow a background parallel field

  !! tell the efield data object to update
  call efield%update(cfg,dtmodel,t,x,ymd,UTsec)

  !! now we need to take data in the efield object and map along the field lines
  call compute_rootBGEfields(x,E02all,E03all)

  !! object double flag to int
  flagdirich=nint(efield%flagdirich)

  !! set boundary condition output arguments based on object data
  !!   Note that we are effectively making copies of object properties; however, additional processing is
  !!   being done here that is specific to the potential solver so that may be justified extra memory use...
  do ix3=1,lx3all
    do ix2=1,lx2all
      Vminx1(ix2,ix3)=efield%Vminx1inow(ix2,ix3)
      Vmaxx1(ix2,ix3)=efield%Vmaxx1inow(ix2,ix3)
    end do
  end do

  !! This forces certain patterns to the boundary conditions to make sure solvers don't get garbage data...
  if (lx2all/=1 .and. lx3all/=1) then
    ! full 3D grid
    do ix3=1,lx3all
      do ix1=1,lx1
        Vminx2(ix1,ix3)=efield%Vminx2isnow(ix3)
        Vmaxx2(ix1,ix3)=efield%Vmaxx2isnow(ix3)
      end do
    end do

    do ix2=1,lx2all
      do ix1=1,lx1
        Vminx3(ix1,ix2)=efield%Vminx3isnow(ix2)
        Vmaxx3(ix1,ix2)=efield%Vmaxx3isnow(ix2)
      end do
    end do
  else
    ! some type of 2D grid, lateral boundary will be overwritten
    Vminx2 = 0
    Vmaxx2 = 0
    if (flagdirich==1) then
      ! Dirichlet:  needs to be the same as the physical top corner grid points
      if (gridflag/=1) then   !non-inverted so logical end is the max alt.
        if (lx2all==1) then
          do ix1=1,lx1
            Vminx3(ix1,:)=Vmaxx1(:,1)
            Vmaxx3(ix1,:)=Vmaxx1(:,lx3all)
          end do
          Vminx2=0
          Vmaxx2=0
        else
          do ix1=1,lx1
            Vminx2(ix1,:)=Vmaxx1(1,:)
            Vmaxx2(ix1,:)=Vmaxx1(lx2all,:)
          end do
          Vminx3=0
          Vmaxx3=0
        end if
      else                    !inverted so logical beginning is max alt.
        if (lx2all==1) then
          do ix1=1,lx1
            Vminx3(ix1,:)=Vminx1(:,1)
            Vmaxx3(ix1,:)=Vminx1(:,lx3all)
          end do
          Vminx2=0
          Vmaxx2=0
        else
          do ix1=1,lx1
            Vminx2(ix1,:)=Vminx1(1,:)
            Vmaxx2(ix1,:)=Vminx1(lx2all,:)
          end do
          Vminx3=0
          Vmaxx3=0
        end if
      end if
    else
      ! Neumann in x1:  sides are grounded...
      if (lx2all==1) then
        Vminx3 = 0
        Vmaxx3 = 0
      else
        Vminx2=0
        Vmaxx2=0
      end if
    end if
  end if
end subroutine potentialBCs2D_fileinput

subroutine compute_rootBGEfields(x,E02all,E03all)
  !> Returns a background electric field calculation for use by external program units.
  !   This requires that all necessary files, etc. have already been loaded into module
  !   variables.  This is only to be called by a root process as it deals with fullgrid
  !   data.  An interface for workers and root is in the top-level potential module.  This
  !   particular bit of code is needed both when setting boundary conditions and also when
  !   initializing background electric field; hence it is a subroutine as opposed to block of code
  class(curvmesh), intent(in) :: x
  real(wp), dimension(:,:,:), intent(inout) :: E02all,E03all
  !! intent(out)
  integer :: ix1,ix2,ix3
  real(wp) :: h2ref,h3ref

  !! the only danger here is that this routine could be called before any module data are loaded
  !   so check just to make sure it isn't being misused in this way
  !!    FIXME: does this accomplish anything???
  if (.not. associated(efield%E0xinow)) error stop  &
        'potentialBCs:compute_rootBGEfields is trying to access unallocated module data'

  !! recompute reference locations here (also computed in object)
  if (lx2all > 1 .and. lx3all>1) then ! 3D sim
    ix2ref = lx2all/2      !note integer division
    ix3ref = lx3all/2
  else if (lx2all==1 .and. lx3all>1) then
    ix2ref = 1
    ix3ref=lx3all/2
  else if (lx2all>1 .and. lx3all==1) then
    ix2ref=lx2all/2
    ix3ref=1
  else
    error stop 'Unable to orient boundary conditions for electric potential'
  endif

  !! by default the code uses 300km altitude as a reference location, using the center x2,x3 point
  !! These are the coordinates for inputs varying along axes 2,3
  ix1ref = minloc(abs(x%rall(:,ix2ref,ix3ref) - Re - 300e3_wp), dim=1)

  !! scale electric fields at some reference point into the full grid
  do ix3=1,lx3all
    do ix2=1,lx2all
      h2ref=x%h2all(ix1ref,ix2,ix3)
      !! define a reference metric factor for a given field line
      h3ref=x%h3all(ix1ref,ix2,ix3)
      do ix1=1,lx1
        E02all(ix1,ix2,ix3)=efield%E0xinow(ix2,ix3)*h2ref/x%h2all(ix1,ix2,ix3)
        E03all(ix1,ix2,ix3)=efield%E0yinow(ix2,ix3)*h3ref/x%h3all(ix1,ix2,ix3)
      end do
    end do
  end do
end subroutine compute_rootBGEfields

subroutine potentialBCs2D(UTsec,cfg,x,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3, &
                                      Vmaxx3,E01all,E02all,E03all,flagdirich)
  ! This is a default routine for setting electromagnetic boundary conditions in cases where user file input is not specified.  It also computes the equatorial vertical drift for the EIA if requested by the user. This routine *could* be modified to hard-code specific conditions in if needed but we really recommend using file input for that.
  real(wp), intent(in) :: UTsec
  type(gemini_cfg), intent(in) :: cfg
  class(curvmesh), intent(in) :: x
  real(wp), dimension(:,:), intent(inout), target :: Vminx1,Vmaxx1
  !! intent(out)
  real(wp), dimension(:,:), intent(inout) :: Vminx2,Vmaxx2
  !! intent(out)
  real(wp), dimension(:,:), intent(inout) :: Vminx3,Vmaxx3
  !! intent(out)
  real(wp), dimension(:,:,:), intent(inout) :: E01all,E02all,E03all
  !! intent(out)
  integer, intent(out) :: flagdirich
  real(wp), dimension(1:size(Vmaxx1,1),1:size(Vmaxx1,2)) :: Emaxx1    !pseudo-electric field
  real(wp) :: Phipk
  integer :: ix1,ix2,ix3    !grid sizes are borrowed from grid module
  integer :: im
  !    integer, parameter :: lmodes=8; this type of thing done from input scripts now...
  real(wp) :: phase
  real(wp), dimension(1:size(Vmaxx1,1)) :: x3dev
  real(wp) :: meanx2,sigx2,meanx3,sigx3,meant,sigt,sigcurv,x30amp,varc    !for setting background field
  real(wp), dimension(:,:), pointer :: Vtopalt,Vbotalt
  real(wp) :: vamp,LThrs,veltime,z,glonmer

  !CALCULATE/SET TOP BOUNDARY CONDITIONS
  sigx2=1/20._wp*(x%x2all(lx2all)-x%x2all(1))
  meanx2=0.5_wp*(x%x2all(1)+x%x2all(lx2all))
  sigx3=1/20._wp*(x%x3all(lx3all)-x%x3all(1))    !this requires that all workers have a copy of x3all!!!!
  meanx3=0.5_wp*(x%x3all(1)+x%x3all(lx3all))

  ! FIXME: the pointer swapping to deal with top vs. bottom here is confusing/superfluous; it may be better simply to have the input preparation scripts assign things accordingly, which is what we will do from now on.  For this routine right now it doesn't matter since both just zeroed out anyway...
  if (gridflag/=2) then
    Vtopalt=>Vminx1
    Vbotalt=>Vmaxx1
  else
    Vtopalt=>Vmaxx1
    Vbotalt=>Vminx1
  end if

  Phipk = 0      !pk current density
  flagdirich = 0    !Neumann conditions
  do ix3=1,lx3all
    do ix2=1,lx2all
      Vtopalt(ix2,ix3) = 0
    end do
  end do

  !SOME USER INFO
  if (debug) print *, 'At time (UT seconds):  ',UTsec,'  Max FAC set to be:  ',maxval(abs(Vtopalt))

  !BOTTOM BOUNDARY IS ALWAYS ZERO CURRENT - SIDES ARE JUST GROUNDED
  Vbotalt = 0   !since we need to have no current through bottom boundary
  Vminx2 = 0
  Vmaxx2 = 0
  Vminx3 = 0
  Vmaxx3 = 0

  !PI's EIA code COMPUTE SOURCE/FORCING TERMS FROM BACKGROUND FIELDS, ETC.
  if (cfg%flagEIA) then
    if (ix1eq<=0) then   !recompute the position of the equator in terms of the x1 variable
      ix1eq = minloc(abs(x%x1), dim=1)    !equator location is that closest to zero in the x1 (q) variable
      if (debug) print*, 'equator ix1:  ',ix1eq,x%x1(ix1eq)
    end if

    vamp=cfg%v0equator    !amplitude of vertical drift at equator from input config.nml file

    E01all=0
    E02all=0

    do ix2=1,lx2all    !for a swapped grid this is longitude
      !for each meridional slice define a local time
      glonmer=x%glonall(ix1eq,ix2,lx3all/2)     !just use halfway up in altitude at the magnetic equator
      do while (glonmer<0)
        glonmer=glonmer+360
      end do

      LThrs=UTsec/3600+glonmer/15                 !Local time of center of meridian
      veltime = sin(2*pi*(LThrs-7)/24)    ! Huba's formulate for velocity amplitude vs. time

      do ix3=1,lx3all     !here this is L-shell
        z = x%altall(ix1eq,ix2,ix3)  !Current altitude of center of this flux tube
        do ix1=1,lx1
          if (z<=150e3_wp) then
            E03all(ix1,ix2,ix3) = 0
          elseif ((z>=150e3_wp) .and. (z<=300e3_wp)) then
            E03all(ix1,ix2,ix3) = -(veltime*vamp*(z-150e3_wp)/150e3_wp)*x%Bmagall(ix1eq,ix2,ix3)    !minus sign to deal with permuted dimensions
          elseif (z>300e3_wp) then
            E03all(ix1,ix2,ix3) = -veltime*vamp*x%Bmagall(ix1eq,ix2,ix3)
          end if
        end do
      end do
    end do

  !  print*, '  Applied EIA perturbation to background electric field...'
  !  print*, '    ',minval(E03all),maxval(E03all)
  else
    E01all = 0
    E02all = 0
    E03all = 0
  end if
end subroutine potentialBCs2D

end module potentialBCs_mumps


! RUN: %S/test_errors.sh %s %t %flang_fc1
subroutine s1(x, y)
  !ERROR: Array pointer 'x' must have deferred shape or assumed rank
  real, pointer :: x(1:)  ! C832
  !ERROR: Allocatable array 'y' must have deferred shape or assumed rank
  real, dimension(1:,1:), allocatable :: y  ! C832
end

subroutine s2(a, b, c)
  real :: a(:,1:)
  real :: b(10,*)
  real :: c(..)
  !ERROR: Array pointer 'd' must have deferred shape or assumed rank
  real, pointer :: d(:,1:)  ! C832
  !ERROR: Allocatable array 'e' must have deferred shape or assumed rank
  real, allocatable :: e(10,*)  ! C832
  !ERROR: Assumed-rank array 'f' must be a dummy argument
  real, pointer :: f(..)  ! C837
  !ERROR: Assumed-shape array 'g' must be a dummy argument
  real :: g(:,1:)
  !ERROR: Assumed-size array 'h' must be a dummy argument
  real :: h(10,*)  ! C833
  !ERROR: Assumed-rank array 'i' must be a dummy argument
  real :: i(..)  ! C837
end

subroutine s3(a, b)
  real :: a(*)
  !ERROR: Dummy array argument 'b' may not have implied shape
  real :: b(*,*)  ! C835, C836
  !ERROR: Implied-shape array 'c' must be a named constant or a dummy argument
  real :: c(*)  ! C836
  !ERROR: Named constant 'd' array must have constant or implied shape
  integer, parameter :: d(:) = [1, 2, 3]
end

subroutine s4()
  type :: t
    integer, allocatable :: a(:)
    !ERROR: Component array 'b' without ALLOCATABLE or POINTER attribute must have explicit shape
    integer :: b(:)  ! C749
    real, dimension(1:10) :: c
    !ERROR: Array pointer component 'd' must have deferred shape
    real, pointer, dimension(1:10) :: d  ! C745
  end type
end

function f()
  !ERROR: Array 'f' without ALLOCATABLE or POINTER attribute must have explicit shape
  real, dimension(:) :: f  ! C832
end

subroutine s5()
  !ERROR: Allocatable array 'a' must have deferred shape or assumed rank
  integer :: a(10), b(:)
  allocatable :: a
  allocatable :: b
end subroutine

subroutine s6()
!C835   An object whose array bounds are specified by an
!  implied-shape-or-assumed-size-spec shall be a dummy data object or a named
!  constant.
!
!C843   An entity with the INTENT attribute shall be a dummy data object or a
!  dummy procedure pointer.
!
!C849   An entity with the OPTIONAL attribute shall be a dummy argument.

  !ERROR: Implied-shape array 'local1' must be a named constant or a dummy argument
  real, dimension (*) :: local1
  !ERROR: INTENT attributes may apply only to a dummy argument
  real, intent(in) :: local2
  !ERROR: INTENT attributes may apply only to a dummy argument
  procedure(), intent(in) :: p1
  !ERROR: OPTIONAL attribute may apply only to a dummy argument
  real, optional :: local3
  !ERROR: OPTIONAL attribute may apply only to a dummy argument
  procedure(), optional :: p2
end subroutine


text:
subroutine plcd_EndsteElmope(b)
   use Mod_plcd_BaseElmope
   use Mod_PLCD
   use Mod_plcd_BMatrix
   use Mod_plcd_BMatrixFactory
   use Mod_plcd_ExacsoError
   implicit none
   class(PLCDProblem), target :: b

   real(rp) :: ielem2

   a=>b

   !ToNullsub
   call SetPointersAndHooksToNULLSUB
   !ResetProcedureComposition
   call ResetProcedureComposition

   !SetPointers
   !Initialize
   call SetPointersExacsoError%Initialize

   !Set
   call SetPointersExacsoError%Set

   !Finalize
   call SetPointersExacsoError%Finalize

   call a%Mesh%ElementAlloc(e,a%Memor,'DefaultRule','plcd_EnditeElmope')
   call a%Mesh%GetNelem(nelem)

   call ProcHook%Initializations
   elements : do ielem = 1,nelem
      call a%Mesh%ElementLoad(ielem,e)

      ielem2 = ielem

      !Move history data to converged
      ElementMatData => a%ElementMaterialsData(ielem)%p
      call ElementMatData%MoveHistoryVariablesToConverged

      !Derivatives at the center of gravity, detjm etc
      call e%elmdcg

      call ProcHook%PreGauss

      do igaus = 1,e%pgaus
         e%igaus = igaus

         call ProcHook%InGaussElmats

      enddo

   enddo elements

   call ProcHook%Finalizations

   !DeallocateElement
   call a%Mesh%ElementDeAlloc(e,a%Memor,'DefaultRule','plcd_Elmope')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
-- Counter is a recursively defined _record_.
Counter = Rec C. { get: Nat, inc: Unit -> C, dec: Unit -> C };

-- An instance c of Counter.
c = let create = fix (lambda f: { x: Nat } -> Counter. lambda s: { x: Nat }.
                        { get = s.x,
                          inc = lambda _: Unit. f { x = succ (s.x) },
                          dec = lambda _: Unit. f { x = pred (s.x) } })
    in (create { x = 0 }) as Counter;

-- The Counters defined here are purely functional...
c1 = (c.inc unit) as Counter;
c2 = (c1.inc unit) as Counter;
c2.get;

max_stars_repo_name: brunoflores/fulluntyped

max_issues_repo_name: brunoflores/fulluntyped

max_forks_repo_name: brunoflores/fulluntyped


text:
      SUBROUTINE RWQBEN2 (TIMTMP)
C
C M. MORTON 01/30/98: CHANGED CODE TO ALLOW FOR TEMPORALLY
C   VARYING BENTHIC FLUXES IN THE BENFN FILE.  PREVIOUS VERSION ONLY
C   PROVIDED SPATIALLY VARYING FLUX (NO PROVISION FOR TIME VARYING).
C CHANGE RECORD
C READ IN SPATIALLY AND/OR TEMPORALLY VARYING PARAMETERS FOR BENTHIC
C FLUXES OF PO4D, NH4, NO3, SAD, COD, O2
C FORMAT OF BENFN FILE IS:
C    TITLE 1
C    TITLE 2
C    TITLE 3
C  270.00000  <-- DAY AT WHICH FOLLOWING FLUXES BECOME ACTIVE
C  350.00000  <-- DAY AT WHICH FOLLOWING FLUXES BECOME ACTIVE
C 9999.99999 <-- ENTER LARGE DAY AT END OF FILE
C
      USE GLOBAL

      CHARACTER TITLE(3)*79, CCMRM*1
      INTEGER,SAVE,ALLOCATABLE,DIMENSION(:)::IZONE
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFCOD
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFNH4
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFNO3
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFO2
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFPO4D
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XBFSAD

      IF(.NOT.ALLOCATED(IZONE  ))THEN
		ALLOCATE(IZONE(NSMZM))
		ALLOCATE(XBFCOD(NSMZM))
		ALLOCATE(XBFNH4(NSMZM))
		ALLOCATE(XBFNO3(NSMZM))
		ALLOCATE(XBFO2(NSMZM))
		ALLOCATE(XBFPO4D(NSMZM))
		ALLOCATE(XBFSAD(NSMZM))
	    IZONE=0
	    XBFCOD=0.0
	    XBFNH4=0.0
	    XBFNO3=0.0
	    XBFO2=0.0
	    XBFPO4D=0.0
	    XBSFAD=0.0
	ENDIF
C
      OPEN(1,FILE=BENFN,STATUS='UNKNOWN')
      OPEN(2,FILE='WQ3D'//ans(partid2)//'.OUT',STATUS='UNKNOWN',POSITION='APPEND')
C
C SKIP OVER THREE HEADER RECORDS:
C
      READ(1,50) (TITLE(M),M=1,3)
      WRITE(2,999)
      WRITE(2,50) (TITLE(M),M=1,3)
C
C SKIP OVER ALL COMMENT CARDS AT BEGINNING OF FILE:
C
      REWIND(1)
      CCMRM = '#'
      CALL SKIPCOMM(1, CCMRM)
      READ(1, *) IBENZ
      WRITE(2, 65) TIMTMP, IBENZ
   65 FORMAT(' * BENTHIC FLUXES AT     ', F10.5,' DAYS OF MODEL RUN',/,
     &    ' NUMBER OF BENTHIC FLUX ZONES = ', I4)
C
C SEQUENTIALLY READ THROUGH BENTHIC FLUX FILE UNTIL THE APPROPRIATE
C TIME IS FOUND:
C   BDAY   = CURRENT DAY AT WHICH BENTHIC FLUX IS IN EFFECT
C   BENDAY = NEXT DAY AT WHICH BENTHIC FLUX CHANGES (PASSED TO MAIN PROG
C
   10 READ(1, *, END=15) BENDAY
      IF(BENDAY .GT. TIMTMP) GOTO 20
      BDAY = BENDAY
      DO I=1,IBENZ
        READ(1,*,END=15) MM, XBFPO4D(MM), XBFNH4(MM), XBFNO3(MM),
     &      XBFSAD(MM), XBFCOD(MM), XBFO2(MM)
        IZONE(I) = MM
      ENDDO
      GOTO 10
C
C UNEXPECTED END-OF-FILE ENCOUNTERED:
C
   15 WRITE(2,16) BENFN
   16 FORMAT(//,' ************* WARNING *************',/,
     &    ' END-OF-FILE ENCOUNTERED IN FILE: ', A20,/,/
     &    ' BENTHIC FLUXES SET TO VALUES CORRESPONDING ',
     &    ' TO LAST DAY IN FILE.',/)
      BENDAY=(TCON*TBEGIN + NTC*TIDALP)/86400.0  ! *** PMC SINGLE LINE
   20 CONTINUE
      WRITE(2, 48) BDAY
   48 FORMAT(/,' DAY IN BENTHIC FLUX FILE: ',F10.5,/,
     &    '    ZONE    FPO4    FNH4    FNO3    FSAD    FCOD    FSOD')
      DO I=1,IBENZ
        MM = IZONE(I)
        WRITE(2,51) MM,XBFPO4D(MM),XBFNH4(MM),XBFNO3(MM),XBFSAD(MM),
     &      XBFCOD(MM),XBFO2(MM)
      ENDDO
C
C DETERMINE BENTHIC FLUX FOR EACH CELL (L) BY INTERPOLATING BETWEEN
C THE MUD AND SAND FLUXES.  XBENMUD(L) IS THE PERCENT MUD FOR EACH
C CELL.
C
      DO L=2,LA
        IZM = IBENMAP(L,1)
        IZS = IBENMAP(L,2)
        XM = XBENMUD(L)
        WQBFPO4D(L) = XM*XBFPO4D(IZM) + (1.0-XM)*XBFPO4D(IZS)
        WQBFNH4(L)  = XM*XBFNH4(IZM)  + (1.0-XM)*XBFNH4(IZS)
        WQBFNO3(L)  = XM*XBFNO3(IZM)  + (1.0-XM)*XBFNO3(IZS)
        WQBFSAD(L)  = XM*XBFSAD(IZM)  + (1.0-XM)*XBFSAD(IZS)
        WQBFCOD(L)  = XM*XBFCOD(IZM)  + (1.0-XM)*XBFCOD(IZS)
        WQBFO2(L)   = XM*XBFO2(IZM)   + (1.0-XM)*XBFO2(IZS)
      ENDDO
      CLOSE(1)
      CLOSE(2)
  999 FORMAT(1X)
   50 FORMAT(A79)
   51 FORMAT(I8, 10F8.3)
   52 FORMAT(I7, 1X, A3)
   60 FORMAT(/, A24, I5, A24)
      RETURN
      END

max_stars_repo_name: fearghalodonncha/DeepCurrent

max_issues_repo_name: fearghalodonncha/DeepCurrent

max_forks_repo_name: fearghalodonncha/DeepCurrent


text:
        subroutine itrBCElas (umesh, disp, iBC, BC, iper, ilwork)
c
c----------------------------------------------------------------------
c
c This program satisfies the boundary conditions on the mesh-elastic BCs.
c
c input:
c  disp   (nshg,nelas)  : mesh-elastic displacement
c  iBC    (nshg)        : Boundary Condition Code
c  BC     (nshg,4)      : boundary condition constraint parameters
c                         pass only 4 componenets to this subroutine
c
c output:
c  disp   (nshg,nelas)  : Adjusted mesh displacement value(s) corresponding
c                          to a constraint d.o.f.
c
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension umesh(nshg,nelas),        disp(nshg,nelas),
     &            iBC(nshg),                BC(nshg,4)
c
        dimension ilwork(nlwork),           iper(nshg)
c
c.... -----------------> Mesh Elastic Velocity <-----------------------
c.... 3D
c
c.... x1-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 1)
            disp(:,1)  = BC(:,1) - BC(:,2) * disp(:,2)
     &                           - BC(:,3) * disp(:,3)
            umesh(:,1) = disp(:,1) / Delt(1)
          endwhere
c
c.... x2-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 2)
            disp(:,2)  = BC(:,1) - BC(:,2) * disp(:,1)
     &                           - BC(:,3) * disp(:,3)
            umesh(:,2) = disp(:,2) / Delt(1)
          endwhere
c
c.... x1-velocity and x2-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 3)
            disp(:,1)  = BC(:,1) - BC(:,2) * disp(:,3)
            disp(:,2)  = BC(:,3) - BC(:,4) * disp(:,3)
            umesh(:,1) = disp(:,1) / Delt(1)
            umesh(:,2) = disp(:,2) / Delt(1)
          endwhere
c
c.... x3-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 4)
            disp(:,3)  = BC(:,1) - BC(:,2) * disp(:,1)
     &                           - BC(:,3) * disp(:,2)
            umesh(:,3) = disp(:,3) / Delt(1)
          endwhere
c
c.... x1-velocity and x3-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 5)
            disp(:,1)  = BC(:,1) - BC(:,2) * disp(:,2)
            disp(:,3)  = BC(:,3) - BC(:,4) * disp(:,2)
            umesh(:,1) = disp(:,1) / Delt(1)
            umesh(:,3) = disp(:,3) / Delt(1)
          endwhere
c
c.... x2-velocity and x3-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 6)
            disp(:,2)  = BC(:,1) - BC(:,2) * disp(:,1)
            disp(:,3)  = BC(:,3) - BC(:,4) * disp(:,1)
            umesh(:,2) = disp(:,2) / Delt(1)
            umesh(:,3) = disp(:,3) / Delt(1)
          endwhere
c
c.... x1-velocity, x2-velocity and x3-velocity, 3D
c
          where (ibits(iBC,14,3) .eq. 7)
            disp(:,1)  = BC(:,1)
            disp(:,2)  = BC(:,2)
            disp(:,3)  = BC(:,3)
            umesh(:,1) = disp(:,1) / Delt(1)
            umesh(:,2) = disp(:,2) / Delt(1)
            umesh(:,3) = disp(:,3) / Delt(1)
          endwhere
c
c       endif
c
c.... end of velocity
c
c.... communications
c
        if (numpe > 1) then
           call commu (disp, ilwork, nelas, 'out')
        endif
c
c       slave has masters value, for abc we need to rotate it
c
        if(iabc==1) then        !are there any axisym bc's
           call rotabc(disp, iBC, 'out')
        endif
c
c.... return
c
        return
        end

max_stars_repo_name: SCOREC/phasta

max_issues_repo_name: SCOREC/phasta

max_forks_repo_name: SCOREC/phasta


!*==zlarfgp.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZLARFGP + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
!
!       .. Scalar Arguments ..
!       INTEGER            INCX, N
!       COMPLEX*16         ALPHA, TAU
!       ..
!       .. Array Arguments ..
!       COMPLEX*16         X( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZLARFGP generates a complex elementary reflector H of order n, such
!> that
!>
!>       H**H * ( alpha ) = ( beta ),   H**H * H = I.
!>              (   x   )   (   0  )
!>
!> where alpha and beta are scalars, beta is real and non-negative, and
!> x is an (n-1)-element complex vector.  H is represented in the form
!>
!>       H = I - tau * ( 1 ) * ( 1 v**H ) ,
!>                     ( v )
!>
!> where tau is a complex scalar and v is a complex (n-1)-element
!> vector. Note that H is not hermitian.
!>
!> If the elements of x are all zero and alpha is real, then tau = 0
!> and H is taken to be the unit matrix.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the elementary reflector.
!> \endverbatim
!>
!> \param[in,out] ALPHA
!> \verbatim
!>          ALPHA is COMPLEX*16
!>          On entry, the value alpha.
!>          On exit, it is overwritten with the value beta.
!> \endverbatim
!>
!> \param[in,out] X
!> \verbatim
!>          X is COMPLEX*16 array, dimension
!>                         (1+(N-2)*abs(INCX))
!>          On entry, the vector x.
!>          On exit, it is overwritten with the vector v.
!> \endverbatim
!>
!> \param[in] INCX
!> \verbatim
!>          INCX is INTEGER
!>          The increment between elements of X. INCX > 0.
!> \endverbatim
!>
!> \param[out] TAU
!> \verbatim
!>          TAU is COMPLEX*16
!>          The value tau.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date November 2017
!
!> \ingroup complex16OTHERauxiliary
!
!  =====================================================================
      SUBROUTINE ZLARFGP(N,Alpha,X,Incx,Tau)
      IMPLICIT NONE
!*--ZLARFGP108
!
!  -- LAPACK auxiliary routine (version 3.8.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     November 2017
!
!     .. Scalar Arguments ..
      INTEGER Incx , N
      COMPLEX*16 Alpha , Tau
!     ..
!     .. Array Arguments ..
      COMPLEX*16 X(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION TWO , ONE , ZERO
      PARAMETER (TWO=2.0D+0,ONE=1.0D+0,ZERO=0.0D+0)
!     ..
!     .. Local Scalars ..
      INTEGER j , knt
      DOUBLE PRECISION alphi , alphr , beta , bignum , smlnum , xnorm
      COMPLEX*16 savealpha
!     ..
!     .. External Functions ..
      DOUBLE PRECISION DLAMCH , DLAPY3 , DLAPY2 , DZNRM2
      COMPLEX*16 ZLADIV
      EXTERNAL DLAMCH , DLAPY3 , DLAPY2 , DZNRM2 , ZLADIV
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , DCMPLX , DIMAG , SIGN
!     ..
!     .. External Subroutines ..
      EXTERNAL ZDSCAL , ZSCAL
!     ..
!     .. Executable Statements ..
!
      IF ( N<=0 ) THEN
         Tau = ZERO
         RETURN
      ENDIF
!
      xnorm = DZNRM2(N-1,X,Incx)
      alphr = DBLE(Alpha)
      alphi = DIMAG(Alpha)
!
      IF ( xnorm/=ZERO ) THEN
!
!        general case
!
         beta = SIGN(DLAPY3(alphr,alphi,xnorm),alphr)
         smlnum = DLAMCH('S')/DLAMCH('E')
         bignum = ONE/smlnum
!
         knt = 0
         IF ( ABS(beta)<smlnum ) THEN
            DO
!
!           XNORM, BETA may be inaccurate; scale X and recompute them
!
               knt = knt + 1
               CALL ZDSCAL(N-1,bignum,X,Incx)
               beta = beta*bignum
               alphi = alphi*bignum
               alphr = alphr*bignum
               IF ( (ABS(beta)>=smlnum) .OR. (knt>=20) ) THEN
!
!           New BETA is at most 1, at least SMLNUM
!
                  xnorm = DZNRM2(N-1,X,Incx)
                  Alpha = DCMPLX(alphr,alphi)
                  beta = SIGN(DLAPY3(alphr,alphi,xnorm),alphr)
                  EXIT
               ENDIF
            ENDDO
         ENDIF
         savealpha = Alpha
         Alpha = Alpha + beta
         IF ( beta<ZERO ) THEN
            beta = -beta
            Tau = -Alpha/beta
         ELSE
            alphr = alphi*(alphi/DBLE(Alpha))
            alphr = alphr + xnorm*(xnorm/DBLE(Alpha))
            Tau = DCMPLX(alphr/beta,-alphi/beta)
            Alpha = DCMPLX(-alphr,alphi)
         ENDIF
         Alpha = ZLADIV(DCMPLX(ONE),Alpha)
!
         IF ( ABS(Tau)<=smlnum ) THEN
!
!           In the case where the computed TAU ends up being a denormalized number,
!           it loses relative accuracy. This is a BIG problem. Solution: flush TAU
!           to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
!
!           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
!           (Thanks Pat. Thanks MathWorks.)
!
            alphr = DBLE(savealpha)
            alphi = DIMAG(savealpha)
            IF ( alphi/=ZERO ) THEN
               xnorm = DLAPY2(alphr,alphi)
               Tau = DCMPLX(ONE-alphr/xnorm,-alphi/xnorm)
               DO j = 1 , N - 1
                  X(1+(j-1)*Incx) = ZERO
               ENDDO
               beta = xnorm
            ELSEIF ( alphr>=ZERO ) THEN
               Tau = ZERO
            ELSE
               Tau = TWO
               DO j = 1 , N - 1
                  X(1+(j-1)*Incx) = ZERO
               ENDDO
               beta = -savealpha
            ENDIF
!
         ELSE
!
!           This is the general case.
!
            CALL ZSCAL(N-1,Alpha,X,Incx)
!
         ENDIF
!
!        If BETA is subnormal, it may lose relative accuracy
!
         DO j = 1 , knt
            beta = beta*smlnum
         ENDDO
         Alpha = beta
!
!        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
!
      ELSEIF ( alphi/=ZERO ) THEN
!           Only "reflecting" the diagonal entry to be real and non-negative.
         xnorm = DLAPY2(alphr,alphi)
         Tau = DCMPLX(ONE-alphr/xnorm,-alphi/xnorm)
         DO j = 1 , N - 1
            X(1+(j-1)*Incx) = ZERO
         ENDDO
         Alpha = xnorm
      ELSEIF ( alphr>=ZERO ) THEN
!              When TAU.eq.ZERO, the vector is special-cased to be
!              all zeros in the application routines.  We do not need
!              to clear it.
         Tau = ZERO
      ELSE
!              However, the application routines rely on explicit
!              zero checks when TAU.ne.ZERO, and we must clear X.
         Tau = TWO
         DO j = 1 , N - 1
            X(1+(j-1)*Incx) = ZERO
         ENDDO
         Alpha = -Alpha
      ENDIF
!
!
!     End of ZLARFGP
!
      END SUBROUTINE ZLARFGP


    module polygon_writer
    use polygon_data

    implicit none

    contains

    !*********************************************************
    subroutine print_polygon(poly, message)
    type(Polygon), intent(in)               :: poly
    character(*), intent(in), optional      :: message

    integer         :: i

    write(*,'(A)') '---------------------'
    if (present(message)) then
        write(*,'(A)') trim(adjustl(message))
    end if

    do i = 1, poly%nsides
        write(*,'(i0,a,f8.3,a,f8.3,a,f8.3,a,f8.3,a)') i, ' : (', poly%edges(i)%start%x, ', ', poly%edges(i)%start%y, ')  -->  (',&
             poly%edges(i)%finish%x, ', ', poly%edges(i)%finish%y, ')'
    end do

    write(*,'(A)') '---------------------'
    end subroutine print_polygon

    !*********************************************************
    subroutine print_list_of_points(points, message)
    type(Vertex), allocatable, intent(in)       :: points(:)
    character(*), intent(in), optional          :: message

    integer                     :: i, n

    n = size(points, 1)

    write(*,'(A)') '---------------------'
    if (present(message)) then
        write(*,'(A)') trim(adjustl(message))
    end if

    do i = 1, n
        write(*,'(f8.3,2X,f8.3)') points(i)%x, points(i)%y
    end do

    write(*,'(A)') '---------------------'
    end subroutine print_list_of_points

    !*********************************************************
    subroutine save_polygon_to_file(fname, poly)
    character(*), intent(in)            :: fname
    type(Polygon), intent(in)           :: poly

    integer             :: nsides, i
    integer             :: fid

    nsides = poly%nsides

    open(newunit=fid, file=fname)
    write(fid, '(i0)') nsides
    do i = 1, nsides
        write(fid, '(i0,f12.4,f12.4)') i, poly%edges(i)%start%x, poly%edges(i)%start%y
    end do

    close(fid)
    end subroutine save_polygon_to_file

    end module polygon_writer


!***********************************************************************************************
      SUBROUTINE ATMOS_DEPOSITION(jyear,redn_Part,redn_TP,redn_TN,N_Atm,P_Atm,Part_Atm)
!***********************************************************************************************
!           -------- Atmospheric Deposition --------
!    Nitrogen forms:   NH4, NO3, DON
!    Phosphorus forms: DOP, DRP, SRP, POP, THP (=BAP)
!    Constant rates are based on LTADS (2005) report and TERC DATA
!************************************************************************
      USE DLMWQ_VARIABLES
      USE DLMWQ_SERVICES
		IMPLICIT NONE
      SAVE

		INTEGER*4 i,j,k,jyear
!
!		Nutrient Budget Control
!
		REAL*8 N_Atm,P_Atm, Part_Atm(7)
!
!		Atmospheric Neutrients deposition rates for Lake Tahoe (mg/m2/day)
!		Units should be in mg/m2/day so that it will match with mg/m3 or ug/L
  		REAL*8 surfacelayer_depth
!---------Nitrogen Neutrients-----------------------
		REAL*8 rate_NO3, rate_NH4, rate_PON, rate_DON
		REAL*8 rate0_NO3, rate0_NH4, rate0_PON, rate0_DON
		REAL*8 rate1_NO3, rate1_NH4, rate1_PON, rate1_DON
!---------Phophorous Neutrients--------------------
      REAL*8	rate_SRP,rate_THP, rate_POP, rate_DOP
		REAL*8 	rate0_SRP,rate0_POP, rate0_DOP
		REAL*8 	rate1_SRP,rate1_POP, rate1_DOP
!---------Details in season-wise--------------------
!************START DRY Nitrogen**********************************************
!------------Winter Nitrogen Dry---------------------------------------------
      REAL*8 winter_NO3_DRY,winter_NH4_DRY,winter_DON_DRY,winter_PN_DRY

      PARAMETER (winter_NO3_DRY = 0.17884, winter_NH4_DRY = 0.53652,		&
     				  winter_DON_DRY = 0.40060, winter_PN_DRY  = 0.04650)
!-------------Spring Nitrogen Dry --------------------------------------------
		REAL*8 spring_NO3_DRY,spring_NH4_DRY,spring_DON_DRY,spring_PN_DRY

      PARAMETER (spring_NO3_DRY = 0.13757, spring_NH4_DRY = 0.41270,		&
     	           spring_DON_DRY = 0.22891, spring_PN_DRY  = 0.03247)

!-------------Summer Nitrogen Dry --------------------------------------------
		REAL*8 summer_NO3_DRY, summer_NH4_DRY,summer_DON_DRY,summer_PN_DRY

      PARAMETER (summer_NO3_DRY = 0.22723, summer_NH4_DRY = 0.53020,		&
     	           summer_DON_DRY = 0.15593, summer_PN_DRY  = 0.03595)

!-------------Fall Nitrogen Dry --------------------------------------------
		REAL*8 fall_NO3_DRY,fall_NH4_DRY,fall_DON_DRY,fall_PN_DRY

      PARAMETER (fall_NO3_DRY = 0.23201, fall_NH4_DRY = 0.66702,			&
     	           fall_DON_DRY = 0.12041, fall_PN_DRY  = 0.04489)
!********************END DRY Nitrogen***************************************
!********************START DRY Phosphrous***********************************
!--------------Winter (December to June) Phophrous Dry --------------------
		REAL*8	winter_SRP_DRY,winter_POP_DRY,winter_DOP_DRY

		PARAMETER (winter_SRP_DRY = 0.00593,winter_POP_DRY = 0.01423,winter_DOP_DRY = 0.00474)

!-------------Summer (July to November) Phosphorous Dry--------------------
		REAL*8 summer_SRP_DRY,summer_POP_DRY,summer_DOP_DRY

		PARAMETER (summer_SRP_DRY = 0.01341,summer_POP_DRY = 0.02850,summer_DOP_DRY = 0.01174)
!*******************END DRY Phophrous*************************************
!*******************START of WET Nitrogen and Phosphrous*****************
!--------------Wet deposition same for all season---------------------
		REAL*8 NO3_WET,NH4_WET,DON_WET,PN_WET,SRP_WET,POP_WET,DOP_WET

		PARAMETER (NO3_WET = 0.68980618, NH4_WET = 0.65148361,				&
     				  DON_WET = 0.82930032, PN_WET = 0.09044125)

		PARAMETER (SRP_WET = 0.03832257,POP_WET = 0.03832257,DOP_WET = 0.03065805)
!****************END of WET Nitrogen and Phosphrous**********************
!-----------------------------------------------------------------------
!
!    Atmospheric Paricles deposition rates for Lake Tahoe (# particles/day/m2)
		REAL(8):: rate_Particles(7)
		REAL(8):: winter_dry_Particles_rate(7), winter_wet_Particles_rate(7),	&
					 spring_dry_Particles_rate(7), spring_wet_Particles_rate(7),	&
					 summer_dry_Particles_rate(7), summer_wet_Particles_rate(7),	&
					   fall_dry_Particles_rate(7),   fall_wet_Particles_rate(7)

		DATA winter_dry_Particles_rate /9.072598E+07,2.996198E+07,					&
			 3.932766E+06,3.894573E+06,6.112662E+05,1.180736E+05,0.00E+00/

		DATA winter_wet_Particles_rate /8.116250E+08,2.680367E+08,					&
			 8.154126E+06,8.072247E+06,2.389807E+05,4.776825E+04,0.00E+00/

		DATA spring_dry_Particles_rate /5.418927E+07,1.789584E+07,					&
			 3.467331E+06,3.429486E+06,5.013876E+05,9.704511E+04,0.00E+00/

		DATA spring_wet_Particles_rate /8.113205E+08,2.679361E+08,					&
			 1.936294E+07,1.918008E+07,2.583086E+05,5.241607E+04,0.00E+00/

		DATA summer_dry_Particles_rate /6.780904E+07,2.239373E+07,					&
			 3.037390E+06,2.990081E+06,6.304551E+05,1.246261E+05,0.00E+00/

		DATA summer_wet_Particles_rate /3.868652E+08,1.277610E+08,					&
			 5.181073E+06,4.890468E+06,2.202288E+05,5.227175E+04,0.00E+00/

		DATA fall_dry_Particles_rate /8.827393E+07,2.915219E+07,						&
			 3.913478E+06,3.702312E+06,4.622116E+05,8.972125E+04,0.00E+00/

		DATA fall_wet_Particles_rate /3.167317E+07,1.045997E+07,						&
			 5.675661E+05,5.514941E+05,2.555857E+04,5.609794E+03,0.00E+00/
!-------------------------------------------------------------------------
		REAL*8 redn, redn_part, redn_TP, redn_TN

!-----------------------------------------------------------------
!-------Initialize to zero to avoid cummulating the values in every time step-----
		rate_NO3 = 0.0
		rate_NH4 = 0.0
		rate_DON = 0.0
		rate_PON = 0.0
		rate_SRP = 0.0
		rate_POP = 0.0
		rate_DOP = 0.0
		rate_THP = 0.0

		rate0_NO3 = 0.0
		rate0_NH4 = 0.0
		rate0_DON = 0.0
		rate0_PON = 0.0
		rate0_SRP = 0.0
		rate0_POP = 0.0
		rate0_DOP = 0.0

		rate1_NO3 = 0.0
		rate1_NH4 = 0.0
		rate1_DON = 0.0
		rate1_PON = 0.0
		rate1_SRP = 0.0
		rate1_POP = 0.0
		rate1_DOP = 0.0

		DO i = 1,7
			rate_Particles(i)=0.0
		ENDDO
!		WRITE(*,10)jday,jyear, redn_part,redn_TP,redn_TN
!10	FORMAT(2i8, 3f10.2, 4x, 'Atmospher')
!		WRITE(1000,*) jyear,jday,redn,'Atmos'
!-------End of initialization----------------------------------------------------
!		redn_part1=1.0d0-0.00d0/100.0d0
!		redn_TP1  =1.0d0-0.00d0/100.0d0
!		redn_TN1  =1.0d0-0.00d0/100.0d0
		WRITE(*,10)jday,jyear, redn_Part,redn_TP,redn_TN
10		FORMAT(2i8, 3f10.2, 4x, 'Atmospher')

!    Depth (m) of the surface layer
!
      surfacelayer_depth = depth(ns)-depth(ns-1)
!***************************************************************************
!***********Inert Atmospheric Particles load into Tahoe Lake****************
!***************************************************************************
!           Winter: (Jan, Feb, March)  90 days   1-90
!           Spring: (April, May, June) 91 days  91-181
!           Summer: (July, Aug, Sep)   92 days 182-273
!           Fall  : (Oct, Nov, Dec)    92 days 274-365
!----------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------
		DO i =1,7
			IF(jday.ge.1.and.jday.le.90) THEN  !Winter
				IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
					rate_Particles(i)=1.95*winter_wet_Particles_rate(i)
				ELSE
					rate_Particles(i)=1.0*winter_dry_Particles_rate(i)
				ENDIF
			ELSEIF(jday.ge.91.and.jday.le.181) THEN	  !Spring
	  			IF(Rain.ge. 2.54)THEN
	    			rate_Particles(i)=1.95*spring_wet_Particles_rate(i)
				ELSE
					rate_Particles(i)=1.0*spring_dry_Particles_rate(i)
				ENDIF
			ELSEIF(jday.ge.182.and.jday.le.273) THEN   !Summer
	  			IF(Rain.ge. 2.54)THEN
	       		rate_Particles(i)=1.95*summer_wet_Particles_rate(i)
				ELSE
					rate_Particles(i)=1.0*summer_dry_Particles_rate(i)
				ENDIF
			ELSE	 !Fall
				IF(Rain.ge. 2.54)THEN
					rate_Particles(i)=1.95*fall_wet_Particles_rate(i)
				ELSE
					rate_Particles(i)=1.0*fall_dry_Particles_rate(i)
				ENDIF
			ENDIF
		ENDDO

		DO i = 1,7
			cf(ns,i) = cf(ns,i) + redn_Part*rate_Particles(i)/(surfacelayer_depth)		   !redn_Part
			Part_Atm(i) = Part_Atm(i) + redn_Part*													&
							 (rate_Particles(i)/(surfacelayer_depth))*vol(ns)*1000.0d0			!redn_Part
	ENDDO

!********************************************************************************
!***********Water Soluble Atmospheric Nutrients load into Tahoe Lake*************
!********************************************************************************
		IF(jday.ge.1.and.jday.le.90) THEN  !Winter  (January to March)
			IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	winter_NO3_DRY
				rate0_NH4 = winter_NH4_DRY
				rate0_DON =	winter_DON_DRY
				rate0_PON =	winter_PN_DRY
				rate0_SRP = winter_SRP_DRY
				rate0_POP = winter_POP_DRY
				rate0_DOP =	winter_DOP_DRY
			ENDIF
		ELSEIF(jday.ge.91.and.jday.le.181) THEN	  !Spring  (April to June)
			IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	spring_NO3_DRY
				rate0_NH4 = spring_NH4_DRY
				rate0_DON =	spring_DON_DRY
				rate0_PON =	spring_PN_DRY
				rate0_SRP = winter_SRP_DRY
				rate0_POP = winter_POP_DRY
				rate0_DOP =	winter_DOP_DRY
			ENDIF
		ELSEIF(jday.ge.182.and.jday.le.273) THEN   !Summer (July to September)
			IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	summer_NO3_DRY
				rate0_NH4 = summer_NH4_DRY
				rate0_DON =	summer_DON_DRY
				rate0_PON =	summer_PN_DRY
				rate0_SRP = summer_SRP_DRY
				rate0_POP = summer_POP_DRY
				rate0_DOP =	summer_DOP_DRY
			ENDIF
		ELSEIF(jday.ge.274.and.jday.le.310) THEN  !Fall
   			IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	fall_NO3_DRY
				rate0_NH4 = fall_NH4_DRY
				rate0_DON =	fall_DON_DRY
				rate0_PON =	fall_PN_DRY
				rate0_SRP = summer_SRP_DRY
				rate0_POP = summer_POP_DRY
				rate0_DOP =	summer_DOP_DRY
			ENDIF

		ELSEIF(jday.ge.311.and.jday.le.334) THEN 	 !Fall	EXCLDUING DECEMBER FOR DRY P (October to November)
			IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	fall_NO3_DRY
				rate0_NH4 = fall_NH4_DRY
				rate0_DON =	fall_DON_DRY
				rate0_PON =	fall_PN_DRY
				rate0_SRP = summer_SRP_DRY
				rate0_POP = summer_POP_DRY
				rate0_DOP =	summer_DOP_DRY
			ENDIF
		ELSE	!Only December
				IF(Rain.ge. 2.54)THEN	 !Wet days for rainfall>=0.1 inch or 2.54 mm
				rate0_NO3 =	NO3_WET
				rate0_NH4 = NH4_WET
				rate0_DON =	DON_WET
				rate0_PON =	PN_WET
				rate0_SRP = SRP_WET
				rate0_POP = POP_WET
				rate0_DOP =	DOP_WET
			ELSE
				rate0_NO3 =	fall_NO3_DRY
				rate0_NH4 = fall_NH4_DRY
				rate0_DON =	fall_DON_DRY
				rate0_PON =	fall_PN_DRY
				rate0_SRP = winter_SRP_DRY
				rate0_POP = winter_POP_DRY
				rate0_DOP =	winter_DOP_DRY
			ENDIF
		ENDIF
!********************Reduction **********************************************
		rate1_SRP = redn_TP*rate0_SRP	!redn_TP
		rate1_POP = redn_TP*rate0_POP	!redn_TP
		rate1_DOP =	redn_TP*rate0_DOP	!redn_TP

		rate1_NO3 =	redn_TN*rate0_NO3	!redn_TN
		rate1_NH4 = redn_TN*rate0_NH4	!redn_TN
		rate1_DON =	redn_TN*rate0_DON	!redn_TN
		rate1_PON =	redn_TN*rate0_PON	!redn_TN
!************CONVERSION FOR BAP AND BAN********************************
!   -------Bioavalable nitrogen and phosphorous----------------------
!    Refractive phosphrous is not bioavailable Quim's Dissertation
!   -------------BAP-------------------------
		rate_THP = 0.95*rate1_SRP + 0.30*rate1_POP + 0.10*rate1_DOP

		rate_SRP = rate1_SRP - 0.95*rate1_SRP
		rate_POP = rate1_POP - 0.30*rate1_POP	 !30%
		rate_DOP = rate1_DOP - 0.10*rate1_DOP	 !10%

!    so 10% of DON and PN is biavailable
		rate_NO3 = rate1_NO3 + 0.55*rate1_PON + 0.55*rate1_DON
		rate_NH4 = rate1_NH4
		rate_PON = rate1_PON - 0.55*rate1_PON
		rate_DON = rate1_DON - 0.55*rate1_DON
!****************END CONVERSION****************************************

		wqual(ns,16) = wqual(ns,16) + rate_THP/(surfacelayer_depth)		!    THP
		wqual(ns,15) = wqual(ns,15) + rate_SRP/(surfacelayer_depth)		!    SRP (wqual (i,12)
      wqual(ns,19) = wqual(ns,19) + rate_POP/(surfacelayer_depth)		!    POP
      wqual(ns,20) = wqual(ns,20) + rate_DOP/(surfacelayer_depth)		!    DOP
      wqual(ns,21) = wqual(ns,21) + rate_NO3/(surfacelayer_depth)		!    NO3
	   wqual(ns,22) = wqual(ns,22) + rate_NH4/(surfacelayer_depth)		!    NH4
		wqual(ns,25) = wqual(ns,25) + rate_PON/(surfacelayer_depth)		!    PON
		wqual(ns,26) = wqual(ns,26) + rate_DON/(surfacelayer_depth)		!    DON
!
!    Nutrient Budget Control
!
      N_Atm=(wqual(ns,21)+wqual(ns,22)+wqual(ns,25)+wqual(ns,26))*vol(ns)
      P_Atm=(wqual(ns,15)+wqual(ns,16)+wqual(ns,19)+wqual(ns,20))*vol(ns)

    	RETURN
		END SUBROUTINE ATMOS_DEPOSITION


text:
        subroutine genBC1 (BCtmp,  iBC,  BC)
c
c----------------------------------------------------------------------
c
c This subroutine adjusts the boundary conditions to accommodate for
c the velocity constraints in the non-axes directions. It copies the
c reduced constraint parameters in BC.
c
c input:
c  BCtmp (nshg,6+5*I3nsd) : input BC parameters (density, temperature,
c                             pressure, (nsd-1)(nsd+1) velocity params,
c                             upto 4 scalar params)
c  iBC   (nshg)           : boundary condition code
c
c output:
c  BC    (nshg,ndofBC)    : the constraint eq's parameters
c
c
c Farzin Shakib, Winter 1986.
c Zdenek Johan,  Winter 1991.  (Fortran 90)
c----------------------------------------------------------------------
c
      include "common.h"

c
        dimension BCtmp(nshg,ndofBC2),    iBC(nshg),
     &            BC(nshg,ndofBC),        tmpbc(4)
c
        dimension tmp(nshg)
c
        integer elas3, elas4, elas5, elas6
c
c.... scalars
c
        do isclr=1,nsclr
           where (btest(iBC,5+isclr)) BC(:,6+isclr) = BCtmp(:,12+isclr)
        enddo
c
c.... set up the thermodynamic properties
c
        where (btest(iBC,0)) BC(:,1) = BCtmp(:,1) ! density
        where (btest(iBC,1)) BC(:,2) = BCtmp(:,2) ! temperature
        where (btest(iBC,2)) BC(:,1) = BCtmp(:,3) ! pressure
c
c.... if the velocity in the x1-direction is specified
c
        where (ibits(iBC,3,3) .eq. 1)
          tmp     = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2
          BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,4)
          BC(:,4) =       BCtmp(:,5) / BCtmp(:,4)
          BC(:,5) =       BCtmp(:,6) / BCtmp(:,4)
        endwhere
c
c.... if the velocity in the x2-direction is specified
c
        where (ibits(iBC,3,3) .eq. 2)
          tmp     = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2
          BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,5)
          BC(:,4) =       BCtmp(:,4) / BCtmp(:,5)
          BC(:,5) =       BCtmp(:,6) / BCtmp(:,5)
        endwhere
c
c.... if the two velocities are specified (x1 & x2-direction)
c
c
c  Protect against user flipping the order of x1 and x2 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),3,3) .eq. 3 .and.
     &       (BCtmp(i,4).eq.0 .or. BCtmp(i,9).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,4:7)
              BCtmp(i,4:7)=BCtmp(i,8:11)
              BCtmp(i,8:11)=tmpbc(1:4)
          endif
        enddo
        where (ibits(iBC,3,3) .eq. 3)
          tmp         = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2
     &                                       + BCtmp(:, 6)**2)
          BCtmp(:, 4) = BCtmp(:, 4) / tmp
          BCtmp(:, 5) = BCtmp(:, 5) / tmp
          BCtmp(:, 6) = BCtmp(:, 6) / tmp
          BCtmp(:, 7) = BCtmp(:, 7) * tmp
c
          tmp         = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2
     &                                       + BCtmp(:,10)**2)
          BCtmp(:, 8) = BCtmp(:, 8) / tmp
          BCtmp(:, 9) = BCtmp(:, 9) / tmp
          BCtmp(:,10) = BCtmp(:,10) / tmp
          BCtmp(:,11) = BCtmp(:,11) * tmp
c
          BCtmp(:, 4) = BCtmp(:, 9) * BCtmp(:, 4)
     &                - BCtmp(:, 5) * BCtmp(:, 8)
          BCtmp(:, 6) = BCtmp(:, 9) * BCtmp(:, 6)
     &                - BCtmp(:, 5) * BCtmp(:,10)
          BCtmp(:, 7) = BCtmp(:, 9) * BCtmp(:, 7)
     &                - BCtmp(:, 5) * BCtmp(:,11)
          BC(:,3)     = BCtmp(:, 7) / BCtmp(:, 4)
          BC(:,4)     = BCtmp(:, 6) / BCtmp(:, 4)
c
          BCtmp(:, 9) = BCtmp(:, 4) * BCtmp(:, 9)
          BCtmp(:,10) = BCtmp(:, 4) * BCtmp(:,10)
     &                - BCtmp(:, 8) * BCtmp(:, 6)
          BCtmp(:,11) = BCtmp(:, 4) * BCtmp(:,11)
     &                - BCtmp(:, 8) * BCtmp(:, 7)
          BC(:,5)     = BCtmp(:,11) / BCtmp(:, 9)
          BC(:,6)     = BCtmp(:,10) / BCtmp(:, 9)
        endwhere
c
c.... if the velocity in the x3-direction is specified
c
        if (nsd .eq. 3) then
        where (ibits(iBC,3,3) .eq. 4)
          tmp     = BCtmp(:,4)**2 + BCtmp(:,5)**2 + BCtmp(:,6)**2
          BC(:,3) = tmp * BCtmp(:,7) / BCtmp(:,6)
          BC(:,4) =       BCtmp(:,4) / BCtmp(:,6)
          BC(:,5) =       BCtmp(:,5) / BCtmp(:,6)
        endwhere
        endif
c
c.... if two velocities are specified (x1 & x3-direction)
c
        if (nsd .eq. 3) then
c
c  Protect against user flipping the order of x1 and x3 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),3,3) .eq. 5 .and.
     &       (BCtmp(i,4).eq.0 .or. BCtmp(i,10).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,4:7)
              BCtmp(i,4:7)=BCtmp(i,8:11)
              BCtmp(i,8:11)=tmpbc(1:4)
           endif
        enddo
        where (ibits(iBC,3,3) .eq. 5)
          tmp         = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2
     &                                       + BCtmp(:, 6)**2)
          BCtmp(:, 4) = BCtmp(:, 4) / tmp
          BCtmp(:, 5) = BCtmp(:, 5) / tmp
          BCtmp(:, 6) = BCtmp(:, 6) / tmp
          BCtmp(:, 7) = BCtmp(:, 7) * tmp
c
          tmp         = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2
     &                                       + BCtmp(:,10)**2)
          BCtmp(:, 8) = BCtmp(:, 8) / tmp
          BCtmp(:, 9) = BCtmp(:, 9) / tmp
          BCtmp(:,10) = BCtmp(:,10) / tmp
          BCtmp(:,11) = BCtmp(:,11) * tmp
c
          BCtmp(:, 4) = BCtmp(:,10) * BCtmp(:, 4)
     &                - BCtmp(:, 6) * BCtmp(:, 8)
          BCtmp(:, 5) = BCtmp(:,10) * BCtmp(:, 5)
     &                - BCtmp(:, 6) * BCtmp(:, 9)
          BCtmp(:, 7) = BCtmp(:,10) * BCtmp(:, 7)
     &                - BCtmp(:, 6) * BCtmp(:,11)
          BC(:,3)     = BCtmp(:, 7) / BCtmp(:, 4)
          BC(:,4)     = BCtmp(:, 5) / BCtmp(:, 4)
c
          BCtmp(:, 9) = BCtmp(:, 4) * BCtmp(:, 9)
     &                - BCtmp(:, 8) * BCtmp(:, 5)
          BCtmp(:,10) = BCtmp(:, 4) * BCtmp(:,10)
          BCtmp(:,11) = BCtmp(:, 4) * BCtmp(:,11)
     &                - BCtmp(:, 8) * BCtmp(:, 7)
          BC(:,5)     = BCtmp(:,11) / BCtmp(:,10)
          BC(:,6)     = BCtmp(:, 9) / BCtmp(:,10)
        endwhere
        endif
c
c.... if two velocities are specified (x2 & x3-direction)
c
        if (nsd .eq. 3) then
c
c  Protect against user flipping the order of x2 and x3 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),3,3) .eq. 6 .and. (
     &       BCtmp(i,5).eq.0 .or. BCtmp(i,10).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,4:7)
              BCtmp(i,4:7)=BCtmp(i,8:11)
              BCtmp(i,8:11)=tmpbc(1:4)
           endif
        enddo
        where (ibits(iBC,3,3) .eq. 6)
          tmp         = sqrt (BCtmp(:, 4)**2 + BCtmp(:, 5)**2
     &                                       + BCtmp(:, 6)**2)
          BCtmp(:, 4) = BCtmp(:, 4) / tmp
          BCtmp(:, 5) = BCtmp(:, 5) / tmp
          BCtmp(:, 6) = BCtmp(:, 6) / tmp
          BCtmp(:, 7) = BCtmp(:, 7) * tmp
c
          tmp         = sqrt (BCtmp(:, 8)**2 + BCtmp(:, 9)**2
     &                                       + BCtmp(:,10)**2)
          BCtmp(:, 8) = BCtmp(:, 8) / tmp
          BCtmp(:, 9) = BCtmp(:, 9) / tmp
          BCtmp(:,10) = BCtmp(:,10) / tmp
          BCtmp(:,11) = BCtmp(:,11) * tmp
c
          BCtmp(:, 4) = BCtmp(:,10) * BCtmp(:, 4)
     &                - BCtmp(:, 6) * BCtmp(:, 8)
          BCtmp(:, 5) = BCtmp(:,10) * BCtmp(:, 5)
     &                - BCtmp(:, 6) * BCtmp(:, 9)
          BCtmp(:, 7) = BCtmp(:,10) * BCtmp(:, 7)
     &                - BCtmp(:, 6) * BCtmp(:,11)
          BC(:,3)     = BCtmp(:, 7) / BCtmp(:, 5)
          BC(:,4)     = BCtmp(:, 4) / BCtmp(:, 5)
c
          BCtmp(:, 8) = BCtmp(:, 5) * BCtmp(:, 8)
     &                - BCtmp(:, 9) * BCtmp(:, 4)
          BCtmp(:,10) = BCtmp(:, 5) * BCtmp(:,10)
          BCtmp(:,11) = BCtmp(:, 5) * BCtmp(:,11)
     &                - BCtmp(:, 9) * BCtmp(:, 7)
          BC(:,5)     = BCtmp(:,11) / BCtmp(:,10)
          BC(:,6)     = BCtmp(:, 8) / BCtmp(:,10)
        endwhere
        endif
c
c.... if all velocities are specified
c

        if (nsd .eq. 3) then
        where (ibits(iBC,3,3) .eq. 7)
          BC(:,3) = BCtmp(:,7) * BCtmp(:,4)
          BC(:,4) = BCtmp(:,7) * BCtmp(:,5)
          BC(:,5) = BCtmp(:,7) * BCtmp(:,6)
        endwhere
        endif

c      write(*,*) ' In genbc1: ibc, bc'
c      do i = 1,nshg
c        if (ibits(iBC(i),3,3) .eq. 7) then
c          write(*,'(a,x,i4,2x,12f6.2)')'comp3:', iBC(i),BCtmp(i,:)
c          write(*,'(a,x,i4,2x,6f10.2)')'comp3:', iBC(i),BC(i,:)
c        endif
c      enddo
c
c.... end
c

c.... if we have the mesh elastic solver
      if (iALE .eq. 2) then
c
c----------------Re-organize Mesh Elastic Velocities-----------------
c.... convert mesh velocity to mesh displacement
c
        elas3 = ndof + 2
        elas4 = ndof + 3
        elas5 = ndof + 4
        elas6 = ndof + 5
c
c.... if the velocity in the x1-direction is specified
c
        where (ibits(iBC,14,3) .eq. 1)
          tmp     = BCtmp(:,17)**2 + BCtmp(:,18)**2 + BCtmp(:,19)**2
          BC(:,elas3) = tmp * BCtmp(:,20) / BCtmp(:,17) * Delt(1)
          BC(:,elas4) =       BCtmp(:,18) / BCtmp(:,17)
          BC(:,elas5) =       BCtmp(:,19) / BCtmp(:,17)
        endwhere
c
c.... if the velocity in the x2-direction is specified
c
        where (ibits(iBC,14,3) .eq. 2)
          tmp     = BCtmp(:,17)**2 + BCtmp(:,18)**2 + BCtmp(:,19)**2
          BC(:,elas3) = tmp * BCtmp(:,20) / BCtmp(:,18) * Delt(1)
          BC(:,elas4) =       BCtmp(:,17) / BCtmp(:,18)
          BC(:,elas5) =       BCtmp(:,19) / BCtmp(:,18)
        endwhere
c
c.... if the two velocities are specified (x1 & x2-direction)
c
c
c  Protect against user flipping the order of x1 and x2 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),14,3) .eq. 3 .and.
     &       (BCtmp(i,17).eq.0 .or. BCtmp(i,22).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,17:20)
              BCtmp(i,17:20)=BCtmp(i,21:24)
              BCtmp(i,21:24)=tmpbc(1:4)
          endif
        enddo
        where (ibits(iBC,14,3) .eq. 3)
          tmp         = sqrt (BCtmp(:, 17)**2 + BCtmp(:, 18)**2
     &                                       + BCtmp(:, 19)**2)
          BCtmp(:, 17) = BCtmp(:, 17) / tmp
          BCtmp(:, 18) = BCtmp(:, 18) / tmp
          BCtmp(:, 19) = BCtmp(:, 19) / tmp
          BCtmp(:, 20) = BCtmp(:, 20) * tmp
c
          tmp         = sqrt (BCtmp(:, 21)**2 + BCtmp(:, 22)**2
     &                                       + BCtmp(:,23)**2)
          BCtmp(:, 21) = BCtmp(:, 21) / tmp
          BCtmp(:, 22) = BCtmp(:, 22) / tmp
          BCtmp(:, 23) = BCtmp(:, 23) / tmp
          BCtmp(:, 24) = BCtmp(:, 24) * tmp
c
          BCtmp(:, 17) = BCtmp(:, 22) * BCtmp(:, 17)
     &                 - BCtmp(:, 18) * BCtmp(:, 21)
          BCtmp(:, 19) = BCtmp(:, 22) * BCtmp(:, 19)
     &                 - BCtmp(:, 18) * BCtmp(:, 23)
          BCtmp(:, 20) = BCtmp(:, 22) * BCtmp(:, 20)
     &                 - BCtmp(:, 18) * BCtmp(:, 23)
          BC(:,elas3)  = BCtmp(:, 20) / BCtmp(:, 17) * Delt(1)
          BC(:,elas4)  = BCtmp(:, 19) / BCtmp(:, 17)
c
          BCtmp(:, 22) = BCtmp(:, 17) * BCtmp(:, 22)
          BCtmp(:, 23) = BCtmp(:, 17) * BCtmp(:, 23)
     &                 - BCtmp(:, 21) * BCtmp(:, 19)
          BCtmp(:, 24) = BCtmp(:, 17) * BCtmp(:, 24)
     &                 - BCtmp(:, 21) * BCtmp(:, 20)
          BC(:,elas5)  = BCtmp(:, 24) / BCtmp(:, 22) * Delt(1)
          BC(:,elas6)  = BCtmp(:, 23) / BCtmp(:, 22)
        endwhere
c
c.... if the velocity in the x3-direction is specified
c
        if (nsd .eq. 3) then
        where (ibits(iBC,3,3) .eq. 4)
          tmp     = BCtmp(:,17)**2 + BCtmp(:,18)**2 + BCtmp(:,19)**2
          BC(:,elas3) = tmp * BCtmp(:,20) / BCtmp(:,19) * Delt(1)
          BC(:,elas4) =       BCtmp(:,17) / BCtmp(:,19)
          BC(:,elas5) =       BCtmp(:,18) / BCtmp(:,19)
        endwhere
        endif
c
c.... if two velocities are specified (x1 & x3-direction)
c
        if (nsd .eq. 3) then
c
c  Protect against user flipping the order of x1 and x3 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),14,3) .eq. 5 .and.
     &       (BCtmp(i,17).eq.0 .or. BCtmp(i,23).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,17:20)
              BCtmp(i,17:20)=BCtmp(i,21:24)
              BCtmp(i,21:24)=tmpbc(1:4)
           endif
        enddo
        where (ibits(iBC,14,3) .eq. 5)
          tmp         = sqrt (BCtmp(:, 17)**2 + BCtmp(:, 18)**2
     &                                       + BCtmp(:, 19)**2)
          BCtmp(:, 17) = BCtmp(:, 17) / tmp
          BCtmp(:, 18) = BCtmp(:, 18) / tmp
          BCtmp(:, 19) = BCtmp(:, 19) / tmp
          BCtmp(:, 20) = BCtmp(:, 20) * tmp
c
          tmp         = sqrt (BCtmp(:, 21)**2 + BCtmp(:, 22)**2
     &                                       + BCtmp(:, 23)**2)
          BCtmp(:, 21) = BCtmp(:, 21) / tmp
          BCtmp(:, 22) = BCtmp(:, 22) / tmp
          BCtmp(:, 23) = BCtmp(:, 23) / tmp
          BCtmp(:, 24) = BCtmp(:, 24) * tmp
c
          BCtmp(:, 17) = BCtmp(:, 23) * BCtmp(:, 17)
     &                 - BCtmp(:, 19) * BCtmp(:, 21)
          BCtmp(:, 18) = BCtmp(:, 23) * BCtmp(:, 18)
     &                 - BCtmp(:, 19) * BCtmp(:, 22)
          BCtmp(:, 20) = BCtmp(:, 23) * BCtmp(:, 20)
     &                 - BCtmp(:, 19) * BCtmp(:, 24)
          BC(:,elas3)  = BCtmp(:, 20) / BCtmp(:, 17) * Delt(1)
          BC(:,elas4)  = BCtmp(:, 18) / BCtmp(:, 17)
c
          BCtmp(:, 22) = BCtmp(:, 17) * BCtmp(:, 22)
     &                 - BCtmp(:, 21) * BCtmp(:, 18)
          BCtmp(:, 23) = BCtmp(:, 17) * BCtmp(:, 23)
          BCtmp(:, 24) = BCtmp(:, 17) * BCtmp(:, 24)
     &                 - BCtmp(:, 21) * BCtmp(:, 20)
          BC(:,elas5)  = BCtmp(:, 24) / BCtmp(:, 23) * Delt(1)
          BC(:,elas6)  = BCtmp(:, 22) / BCtmp(:, 23)
        endwhere
        endif
c
c.... if two velocities are specified (x2 & x3-direction)
c
        if (nsd .eq. 3) then
c
c  Protect against user flipping the order of x2 and x3 in
c  the vector 1 and vector 2.  Without this it will blow up.
c
        do i=1,nshg
          if(ibits(iBC(i),14,3) .eq. 6 .and. (
     &       BCtmp(i,18).eq.0 .or. BCtmp(i,23).eq.0)) then !flip them
              tmpbc(1:4)=BCtmp(i,17:20)
              BCtmp(i,17:20)=BCtmp(i,21:24)
              BCtmp(i,21:24)=tmpbc(1:4)
           endif
        enddo
        where (ibits(iBC,14,3) .eq. 6)
          tmp         = sqrt (BCtmp(:, 17)**2 + BCtmp(:, 18)**2
     &                                       + BCtmp(:, 19)**2)
          BCtmp(:, 17) = BCtmp(:, 17) / tmp
          BCtmp(:, 18) = BCtmp(:, 18) / tmp
          BCtmp(:, 19) = BCtmp(:, 19) / tmp
          BCtmp(:, 20) = BCtmp(:, 20) * tmp
c
          tmp         = sqrt (BCtmp(:, 21)**2 + BCtmp(:, 22)**2
     &                                       + BCtmp(:,23)**2)
          BCtmp(:, 21) = BCtmp(:, 21) / tmp
          BCtmp(:, 22) = BCtmp(:, 22) / tmp
          BCtmp(:, 23) = BCtmp(:, 23) / tmp
          BCtmp(:, 24) = BCtmp(:, 24) * tmp
c
          BCtmp(:, 17) = BCtmp(:, 23) * BCtmp(:, 17)
     &                 - BCtmp(:, 19) * BCtmp(:, 21)
          BCtmp(:, 18) = BCtmp(:, 23) * BCtmp(:, 18)
     &                 - BCtmp(:, 19) * BCtmp(:, 22)
          BCtmp(:, 20) = BCtmp(:, 23) * BCtmp(:, 20)
     &                 - BCtmp(:, 19) * BCtmp(:, 23)
          BC(:,elas3)  = BCtmp(:, 20) / BCtmp(:, 18) * Delt(1)
          BC(:,elas4)  = BCtmp(:, 17) / BCtmp(:, 18)
c
          BCtmp(:, 21) = BCtmp(:, 18) * BCtmp(:, 21)
     &                 - BCtmp(:, 22) * BCtmp(:, 17)
          BCtmp(:, 23) = BCtmp(:, 18) * BCtmp(:, 23)
          BCtmp(:, 24) = BCtmp(:, 18) * BCtmp(:, 24)
     &                 - BCtmp(:, 22) * BCtmp(:, 20)
          BC(:,elas5)  = BCtmp(:, 24) / BCtmp(:, 23) * Delt(1)
          BC(:,elas6)  = BCtmp(:, 21) / BCtmp(:, 23)
        endwhere
        endif
c
c.... if all velocities are specified
c

        if (nsd .eq. 3) then
        where (ibits(iBC,14,3) .eq. 7)
          BC(:,elas3) = BCtmp(:,20) * BCtmp(:,17) * Delt(1)
          BC(:,elas4) = BCtmp(:,20) * BCtmp(:,18) * Delt(1)
          BC(:,elas5) = BCtmp(:,20) * BCtmp(:,19) * Delt(1)
        endwhere
        endif
      endif ! end if iALE = 2
c----------------End Mesh Elastic Velocity---------------
c
c.... end
        return
        end


text:
! A parser for standard IO commands: printf, sprintf.

module stdinoutparser

  use fsystem
  use genoutput
  use storage
  use collection
  use typedsymbol

  use commandparserbase

  implicit none

  private

  public :: sioprs_qualifyString

contains

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

  !<subroutine>

  subroutine sioprs_qualifyString (sstring,Rvariables)

  !<description>
    ! Resolves a string. Replaces all qualifiers in the string by the
    ! values in the variables.
  !</description>

  !<inputoutput>
    ! The string to resolve. Will be replaced by the fully qualified string.
    ! The string should not be dequoted.
    character(len=*), intent(inout) :: sstring
  !</inputoutput>

  !<input>
    ! List of variable values to be inserted into the string.
    type(t_symbolValue), dimension(:), intent(in), optional :: Rvariables
  !</input>

!</subroutine>

    ! local variables
    character(len=len(sstring)) :: soutput
    character(len=SYS_STRLEN) :: sformat
    integer :: iout, iin, ilen, istrlen, inextvar
    logical :: bescape,berror

    ! Loop through the input characters
    bescape = .false.
    iout = 1
    iin = 1
    inextvar = 1
    soutput = ""
    istrlen = len_trim(sstring)
    do
      if (iin .gt. istrlen) exit

      if (bescape) then
        ! Take the next character as it is.
        soutput(iout:iout) = sstring(iin:iin)
        iout = iout + 1
      else
        ! What is the character?
        if (sstring(iin:iin) .eq. "\") then
          ! Escape the next character
          bescape = .true.

          ! Transfer the escape character, so later dequoting works correctly.
          soutput(iout:iout) = sstring(iin:iin)
          iout = iout + 1

        else if (sstring(iin:iin) .eq. "%") then

          ! Get the full format qualifier.
          call getFormat (sstring(iin:),sformat,ilen)

          ! Skip it in the input
          iin = iin + ilen - 1

          ! Plug in the next variable (correctly formatted) and continue
          ! in the output.
          if (present(Rvariables)) then
            if (inextvar .le. ubound(Rvariables,1)) then
              if ((Rvariables(inextvar)%ctype .ne. STYPE_INTEGER) .and. &
                  (Rvariables(inextvar)%ctype .ne. STYPE_DOUBLE) .and. &
                  (Rvariables(inextvar)%ctype .ne. STYPE_STRING)) then
                call output_line ("Ignoring invalid variable.");
                berror = .true.
              else
                call formatVariable (Rvariables(inextvar),soutput(iout:),sformat,ilen,berror)
                inextvar = inextvar + 1
              end if

              ! Ignore wrong qualifiers.
              if (.not. berror) then
                iout = iout + ilen
              end if
            end if
          end if
        else
          ! Transfer the character
          soutput(iout:iout) = sstring(iin:iin)
          iout = iout + 1
        end if

      end if

      bescape = .false.
      iin = iin + 1
    end do

    ! Return the string
    sstring = soutput

  end subroutine

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

  !<subroutine>

  subroutine getFormat (sstring,sformat,ilen)

  !<description>
    ! Local routine. Extracts a format qualifier from a string.
    ! The qualifier must have the format "%... " ending with a space or
    ! the string end.
  !</description>

  !<input>
    ! String to analyse
    character(len=*), intent(in) :: sstring
  !</input>

  !<output>
    ! Format string.
    character(len=*), intent(out) :: sformat

    ! Length of the format string
    integer, intent(out) :: ilen
  !</output>

!</subroutine>

    ! Position of the first format specifier
    ilen = scan(sstring,"dfs")
    if (ilen .eq. 0) ilen = len_trim(sstring)

    sformat = sstring(1:ilen)
  end subroutine

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

  !<subroutine>

  subroutine formatVariable (rvariable,sstring,sformat,ilen,berror)

  !<description>
    ! Local routine. Formats a variable according to the format string.
  !</description>

  !<input>
    ! Variable to format.
    type(t_symbolValue), intent(in) :: rvariable

    ! Format string.
    character(len=*), intent(in) :: sformat
  !</input>

  !<output>
    ! Output string
    character(len=*), intent(inout) :: sstring

    ! Length of the output string
    integer, intent(out) :: ilen

    ! Set to TRUE if the format is wrong.
    logical, intent(out) :: berror
  !</output>

!</subroutine>

    integer :: ilenformat,ipostype,ios
    character (len=SYS_NAMELEN) :: snewformat
    character (len=SYS_STRLEN) :: stemp

    berror = .false.
    ilenformat = len_trim(sformat)

    ! Find the type
    ipostype = scan(sformat,"ldfs")

    ! Verify that the format is ok.
    if (verify(sformat(2:ipostype-1),"0123456789.") .ne. 0) then
      ! Incorrect string format
      berror = .true.
      return
    end if

    if (ipostype .gt. 2) then
      ! Get the format string
      snewformat = sformat(2:ipostype-1)
    else
      snewformat = ""
    end if

    ! What type should we write?
    if (sformat(ilenformat:ilenformat) .eq. "s") then
      ! Just a string.
      if (snewformat .eq. "") then
        ! Without length qualifier
        call cmdprs_dequoteStd (rvariable%svalue,sstring,ilen)
      else
        ! With length qualifier
        if (index(snewformat,".") .eq. 0) then
          read (snewformat,*) ilen
          snewformat = "(A"//trim(snewformat)//")"
          call cmdprs_dequoteStd (rvariable%svalue,stemp)
          write (sstring,snewformat,IOSTAT=ios) stemp
        else
          ! Otherwise: Format invalid.
          ilen = 0
          berror = .true.
        end if
      end if
      return
    end if

    if (sformat(ilenformat:ilenformat) .eq. "d") then
      ! Integer value.
      if (snewformat .ne. "") then
        snewformat = "(I"//trim(snewformat)//")"
        write (sstring,snewformat,IOSTAT=ios) rvariable%ivalue
      else
        write (sstring,*) rvariable%ivalue
        sstring = adjustl(sstring)
      end if
      ilen = len_trim(sstring)
      return
    end if

    if (sformat(ilenformat:ilenformat) .eq. "f") then
      ! Float value.
      if (snewformat .ne. "") then
        snewformat = "(F"//trim(snewformat)//")"
        write (sstring,snewformat,IOSTAT=ios) rvariable%dvalue
      else
        write (sstring,*) rvariable%dvalue
        sstring = adjustl(sstring)
      end if
      ilen = len_trim(sstring)
      return
    end if

    if (sformat(ilenformat:ilenformat) .eq. "e") then
      ! Float value.
      if (snewformat .ne. "") then
        snewformat = "(E"//trim(snewformat)//")"
        write (sstring,snewformat,IOSTAT=ios) rvariable%dvalue
      else
        write (sstring,*) rvariable%dvalue
        sstring = adjustl(sstring)
      end if
      ilen = len_trim(sstring)
      return
    end if

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense