MODULE aerosol_thermodynamics
   IMPLICIT NONE

   ! Prototype module for calculating water content and vapour pressures for
   ! organic species condensing in a chamber study.

   ! Prototype module for efficient calculation of the partitioning of species between
   ! gas and aerosol. Built on a newly developed inorganic compound partitioning code
   ! developed by Dave Topping.
   !
   ! Written by Dave.Topping. Pure organic component properties predicted by Mark Barley
   ! based on VOCs predicted in MCM simulations performed by Mike Jenkin. Delivered by
   ! Gordon McFiggans as Deliverable D22 from WP1.4 in the EU FP6 EUCAARI Integrated Project.
   !
   ! Queries concerning the use of this code through Gordon McFiggans, g.mcfiggans@manchester.ac.uk,
   ! Ownership: D. Topping, Centre for Atmospheric Sciences, University of Manchester, 2007

CONTAINS

   SUBROUTINE inorganic_pdfite(rh,temp,ions,water_total,Press_HNO3,Press_HCl,Press_NH3,gamma_out,mols_out)
      IMPLICIT NONE

   !the following subroutine calculates the water content of a mixed inorganic/organic
   !particle as well as equilibrium vapour pressures above the solution (HNO3,HCl,NH3 and
   !representative organic compounds)
   !Inorganic ions treated (H-NH4-Na-SO4-NO3-Cl) - bisulphate ion concentration calculated
   !as described in a shortly to appear publication.
   !currently, one representative secondary organic compound is treated
   !- commented sections indicate how 16 components (to follow, after completion of D30 and
   !D34). The one representative SOA compound currently carries the abundance weighted
   !average properties from >3000 potentially condensable organic compounds predicted in
   !the base case simulation using the Master Chemical Mechanism (MCM). This module accounts for
   !no interaction between the inorganic & organic components
   !
   !--CALCULATION PROCEDURE----
   !
   !the following code determines inorganic solutes and 'equivalent' mole fractions from the
   !ionic inputs (ions), for use in calculating the water content (water_total) and activity
   !coefficients. after activity coefficients are calculated these are then used for calculating
   !vapour pressures using appropriate equilibrium constants (Press_HNO3,Press_NH3,Press_HCl).
   !the one representative organic is treated as ideal for contributions to the total water and for
   !calculating its equilibrium vapour pressure using raoults law.
   !note: the organic pure component vapour pressures for >3000 potentially condensable organic
   !compounds from the base case simulation of the Master Chemical Mechanism (MCM) are calculated
   !from boiling point predictions. absorptive partitioning of all compounds at their total predicted
   !atmospheric burden to a similar organic mass is used to calculate the relative abundance of organic
   !species. the properties of the single representative organic component in this simulation are
   !calculated as the condensed phase abundance-weighted average of all condensing components
   !
   !the following text describe the sections of code individually which are labelled 1,2,3,4,5
   !
   !1) - COMPOSITION DEFINITIONS
   !
   ! a)Inorganic ion pairing
   !
   !  in order to calculate the water content, which is also used in calculating vapour
   !  pressures, one needs to pair the anions and cations for use in the ZSR mixing rule
   !  The equation provided by Clegg et al (2001) is used for ion pairing.
   !  The solutes chosen comprise of 9 inorganic salts and acids which provide a pairing between
   !  each anion and cation:
   !  (NH4)2SO4,NH4NO3,NH4Cl,Na2SO4,NaNO3,NaCl,H2SO4,HNO3,HCl
   !  the organic compound is treated as a seperate solute
   !
   ! b)Inorganic equivalent fractions
   !  These values are calculated so that activity coefficients can be expressed by
   !  a linear additive rule, thus allowing more efficient calculations and future expansion
   !  (see more detailed description below)
   !
   !2) - WATER CALCULATION
   !
   ! a)The water content is calculated using the ZSR rule with solute concentrations calculated
   !  using 1a above
   !  Whilst the usual approximation of ZSR relies on binary data consisting of 5th or higher order
   !  polynomials, in this code 4 different RH regimes are used, each housing cubic equations for
   !  the water associated with each solute listed above.
   !  Binary water contents for inorganic components were calculated using AIM online (Clegg et al 1998)
   !  The water associated with the organic compound is calculated assuming ideality and that aw=RH
   !
   ! b)molality of each inorganic ion and organic solute (initial input) is calculated for use in
   !  vapour pressure calculation
   !
   !3) - BISULPHATE DISSOCIATION CALCULATION
   !
   ! The dissociation of the bisulphate ion is calculated explicitly
   !  A solution to the equilibrium euqation between the bisulphate ion, hydrogen ion and sulphate ion
   !  is found using tabulated equilibrium constants (referenced)
   !  It is necessary to calculate the activity coefficients of HHSO4 and H2SO4 in a noniterative
   !  manner. these are calulated using the same format as described in 4) below, where both activity
   !  coefficients were fit to the output from ADDEM (Topping et al 2005a,b) covering an extensive
   !  composition space, providing the activity coefficients and bisulphate ion dissociation as a
   !  function of equivalent mole fractions and relative humidity.
   !
   !4) ACTIVITY COEFFICIENTS -for vapour pressures of HNO3,HCl and NH3
   !
   ! This section evaluates activity coefficients and vapour pressures using the water content
   !  calculated above) for each inorganic condensing species
   !	a - HNO3
   !	b - NH3
   !	c - HCl
   !  the following procedure is used:
   !
   !  Zaveri et al (2005) found that one could express the variation of activity coefficients
   !  linearly in log space if equivalent mole fractions were used. So, by a taylor series expansion
   !  log(activity coefficient)=log(binary activity coefficient at a given RH)+
   !   (equivalent mole fraction compound A)*('interaction' parameter between A and condensing specie)+
   !   (equivalent mole fraction compound B)*('interaction' parameter between B and condensing specie)+
   !  here, the interaction parameters have been fit to ADDEM by searching the whole compositon space
   !  and fit usign the Levenberg-Marquardt nonlinear least squares algorithm.
   !
   !  They are given as a function of RH and vary with complexity ranging from linear to 5th order
   !  polynomial expressions, the binary activity coefficients were calculated using AIM online.
   !  note: for NH3, no binary activity coefficient was used and the data were fit to the ratio of
   !  the activity coefficients for the ammonium and hydrogen ions
   !  Once the activity coefficients are obtained the vapour pressure can be easily calculated
   !  using tabulated equilibrium constants (referenced)
   !  This procedure differs from that of Zaveri et al (2005) in that it is not assumed one can carry
   !  behaviour from binary mixtures in multicomponent systems. To this end we have fit the
   ! 'interaction' parameters explicitly to a general inorganic equilibrium model (ADDEM - Topping et
   !  al 2005a,b). Such parameters take into account bisulphate ion dissociation and water content.
   !  This also allows us to consider one regime for all composition space, rather than defining
   !  sulphate rich and sulphate poor regimes
   !
   !5) - The vapour pressure of the condensing organic species are calculated in this prototype by
   !  assuming ideality and using Raoults law. In the final module, fits to organic composition space
   !  within UNIFAC will be used to estimate component activities (including water associated with the
   !  organic compounds) and hence the organic vapour pressures.
   !
   !-----------------------------------------------------------------------------------

      !VARIABLE DECLARATIONS

      ! this array is for the output of the activity coefficients for calculating the
      ! non-ideal dissociation constants
      ! 1: gamma_HNO3
      ! 2: gamma_HCl
      ! 3: gamma_NH4+/gamma_H+ ("gamma_NH3")
      ! 4: (gamma_HHSO4**2)/gamma_H2SO4 (for H2SO4 dissociation)
      ! 5: (gamma_H2SO4**3)/(gamma_HHSO4**2) (for HSO4 dissociation)

              ! Juha: Added the following
              ! 6: gamma_NH4HSO2
              ! 7: gamma_HHSO4

      REAL, DIMENSION(7) :: gamma_out  ! Juha: Changed dimension 5->7

      REAL, DIMENSION(7) :: ions,ions_mol

      REAL, DIMENSION(7) :: mols_out  ! Juha: put out ion molalities

      REAL :: charge_sum, nitric_acid, hydrochloric_acid, sulphuric_acid,&
              ammonium_sulphate, ammonium_nitrate, ammonium_chloride, sodium_sulphate,&
              sodium_nitrate, sodium_chloride, solutes, hydrochloric_acid_eq_frac,&
              sulphuric_acid_eq_frac, ammonium_sulphate_eq_frac, ammonium_nitrate_eq_frac,&
              ammonium_chloride_eq_frac, sodium_sulphate_eq_frac, sodium_nitrate_eq_frac,&
              sodium_chloride_eq_frac, nitric_acid_eq_frac, RH, Temp

      real::hno3_w, hcl_w, h2so4_w, nh42so4_w, nh4no3_w, nh4cl_w,&
            na2so4_w, nano3_w, nacl_w, water_total, h_out, hso4_out, so4_out

      real::nh43hso42_w, nh4hso4_w, na3hso42_w, nahso4_w

      real::hno3_water, hcl_water, h2so4_water, nh42so4_water, nh4no3_water, nh4cl_water,&
            na2so4_water, nano3_water, nacl_water, nh43hso42_water, nh4hso4_water, na3hso42_water, nahso4_water

      real::nNa, nNH4, nSulf, Xt, fNa, fNH4

      integer::binary_case, full_complexity

      real::binary_hno3

      real::HCL_hno3, H2SO4_hno3, NH42SO4_hno3, NH4NO3_hno3, NH4Cl_hno3,&
            Na2SO4_hno3, NaNO3_hno3, NaCl_hno3

      real::Ln_HNO3_act, gamma_hno3, Press_HNO3, K_hno3

      real::binary_hcl

      real::HNO3_hcl, H2SO4_hcl, NH42SO4_hcl, NH4NO3_hcl, NH4Cl_hcl,&
            Na2SO4_hcl, NaNO3_hcl, NaCl_hcl

      real::Ln_HCl_act, gamma_hcl, Press_HCl

      real::HNO3_nh3, HCl_nh3, H2SO4_nh3, NH42SO4_nh3, NH4NO3_nh3, NH4Cl_nh3,&
            Na2SO4_nh3, NaNO3_nh3, NaCl_nh3

      real::Ln_NH3_act, gamma_nh3, Press_NH3

      real::K_hcl

      real::Kh, Knh4, Kw, molality_ratio_nh3

      !----variables used in calculating bisulphate dissociation----

      real::binary_hhso4

      real::HNO3_hhso4, HCL_hhso4, NH42SO4_hhso4, NH4NO3_hhso4, NH4Cl_hhso4,&
            Na2SO4_hhso4, NaNO3_hhso4, NaCl_hhso4

      real::Ln_hhso4_act, gamma_hhso4

      real::binary_h2so4

      real::HNO3_h2so4, HCl_h2so4, NH42SO4_h2so4, NH4NO3_h2so4, NH4Cl_h2so4,&
            Na2SO4_h2so4, NaNO3_h2so4, NaCl_h2so4

      real::Ln_h2so4_act, gamma_h2so4

      real::act_product, a, b, c, root1, root2

      real::h_real, hso4_real, so4_real

      real::Mol_frac_hno3,Mol_frac_hcl,Mol_frac_h2so4,Mol_frac_nh42so4,Mol_frac_nh4no3,&
            Mol_frac_nh4cl,Mol_frac_na2so4,Mol_frac_nano3,Mol_frac_nacl

      real::Mol_frac_nh43hso42, Mol_frac_nh4hso4, Mol_frac_na3hso42, Mol_frac_nahso4

      real:: henrys_temp_dep

      !---new nh3 variables---
      real::binary_nh4hso4,HNO3_nh4hso4,HCL_nh4hso4,H2SO4_nh4hso4,NH42SO4_nh4hso4,&
            NH4NO3_nh4hso4,NH4Cl_nh4hso4,Na2SO4_nh4hso4,NaNO3_nh4hso4,NaCl_nh4hso4
      real::Ln_NH4HSO4_act,gamma_nh4hso4

      !----------------------------------------------------------------------------------
      !
      !VALUE INITIALISATION
      henrys_temp_dep = (1/temp - 1/298d0)

      HCL_hno3=1.0d0;H2SO4_hno3=1.0d0;NH42SO4_hno3=1.0d0;NH4NO3_hno3=1.0d0;NH4Cl_hno3=1.0d0;
      Na2SO4_hno3=1.0d0;NaNO3_hno3=1.0d0;NaCl_hno3=1.0d0;
      HNO3_hcl=1.0d0;H2SO4_hcl=1.0d0;NH42SO4_hcl=1.0d0;NH4NO3_hcl=1.0d0;NH4Cl_hcl=1.0d0;
      Na2SO4_hcl=1.0d0;NaNO3_hcl=1.0d0;NaCl_hcl=1.0d0;
      HNO3_nh3=1.0d0;HCl_nh3=1.0d0;H2SO4_nh3=1.0d0;NH42SO4_nh3=1.0d0;NH4NO3_nh3=1.0d0;
      NH4Cl_nh3=1.0d0;Na2SO4_nh3=1.0d0;NaNO3_nh3=1.0d0;NaCl_nh3=1.0d0;
      HNO3_hhso4=1.0d0;HCL_hhso4=1.0d0;NH42SO4_hhso4=1.0d0;NH4NO3_hhso4=1.0d0;NH4Cl_hhso4=1.0d0;
      Na2SO4_hhso4=1.0d0;NaNO3_hhso4=1.0d0;NaCl_hhso4=1.0d0
      HNO3_h2so4=1.0d0;HCl_h2so4=1.0d0;NH42SO4_h2so4=1.0d0;NH4NO3_h2so4=1.0d0;NH4Cl_h2so4=1.0d0;
      Na2SO4_h2so4=1.0d0;NaNO3_h2so4=1.0d0;NaCl_h2so4=1.0d0;
      !--new nh3 variables--
      HNO3_nh4hso4=1.0d0;HCL_nh4hso4=1.0d0;H2SO4_nh4hso4=1.0d0;NH42SO4_nh4hso4=1.0d0;
      NH4NO3_nh4hso4=1.0d0;NH4Cl_nh4hso4=1.0d0;Na2SO4_nh4hso4=1.0d0;NaNO3_nh4hso4=1.0d0;NaCl_nh4hso4=1.0d0;

      ! Juha: added
      mols_out = 0.d0

      Press_HNO3=0.0d0
      Press_HCl=0.0d0
      Press_NH3=0.0d0  !Initialising vapour pressure over the multicomponent particle
      gamma_out = 1d0  ! i.e. don't alter the ideal mixing ratios if there's nothing there.

      !------------------------------------------------------------------------
      !
      !1)-COMPOSITION DEFINITIONS
      !
      ! a)Inorganic ion pairing
      !  pair cations and anions into solutes according to Clegg et al (2001)

      charge_sum=ions(1)+ions(2)+ions(3)+2.0d0*ions(4)+ions(5)+ions(6)+ions(7)
      nitric_acid=0.0d0;hydrochloric_acid=0.0d0;sulphuric_acid=0.0d0
      ammonium_sulphate=0.0d0;ammonium_nitrate=0.0d0;ammonium_chloride=0.0d0
      sodium_sulphate=0.0d0;sodium_nitrate=0.0d0;sodium_chloride=0.0d0
      nitric_acid=(2.0d0*ions(1)*ions(6)*((1.0d0/1.0d0)**0.5))/(charge_sum)
      hydrochloric_acid=(2.0d0*ions(1)*ions(7)*((1.0d0/1.0d0)**0.5))/(charge_sum)
      sulphuric_acid=(2.0d0*ions(1)*ions(4)*((2.0d0/2.0d0)**0.5))/(charge_sum)
      ammonium_sulphate=(2.0d0*ions(2)*ions(4)*((2.0d0/2.0d0)**0.5))/(charge_sum)
      ammonium_nitrate=(2.0d0*ions(2)*ions(6)*((1.0d0/1.0d0)**0.5))/(charge_sum)
      ammonium_chloride=(2.0d0*ions(2)*ions(7)*((1.0d0/1.0d0)**0.5))/(charge_sum)
      sodium_sulphate=(2.0d0*ions(3)*ions(4)*((2.0d0/2.0d0)**0.5))/(charge_sum)
      sodium_nitrate=(2.0d0*ions(3)*ions(6)*((1.0d0/1.0d0)**0.5))/(charge_sum)
      sodium_chloride=(2.0d0*ions(3)*ions(7)*((1.0d0/1.0d0)**0.5))/(charge_sum)

      ! b) - Inorganic equivalent fractions

      solutes=0.0d0
      solutes=3.0d0*sulphuric_acid+2.0d0*hydrochloric_acid+2.0d0*nitric_acid+3.0d0*ammonium_sulphate+&
           &2.0d0*ammonium_nitrate+2.0d0*ammonium_chloride+3.0d0*sodium_sulphate+2.0d0*sodium_nitrate+&
           &2.0d0*sodium_chloride

      nitric_acid_eq_frac=2.0d0*nitric_acid/(solutes)
      hydrochloric_acid_eq_frac=2.0d0*hydrochloric_acid/(solutes)
      sulphuric_acid_eq_frac=3.0d0*sulphuric_acid/(solutes)
      ammonium_sulphate_eq_frac=3.0d0*ammonium_sulphate/(solutes)
      ammonium_nitrate_eq_frac=2.0d0*ammonium_nitrate/(solutes)
      ammonium_chloride_eq_frac=2.0d0*ammonium_chloride/(solutes)
      sodium_sulphate_eq_frac=3.0d0*sodium_sulphate/(solutes)
      sodium_nitrate_eq_frac=2.0d0*sodium_nitrate/(solutes)
      sodium_chloride_eq_frac=2.0d0*sodium_chloride/(solutes)

      !--inorganic ion molalities

      ions_mol(:)=0.0d0
      ions_mol(1)=ions(1)/(water_total*18.01528d-3)  !H
      ions_mol(2)=ions(2)/(water_total*18.01528d-3)  !NH4
      ions_mol(3)=ions(3)/(water_total*18.01528d-3)  !Na
      ions_mol(4)=ions(4)/(water_total*18.01528d-3)  !SO4
      ions_mol(5)=ions(5)/(water_total*18.01528d-3)  !HSO4
      ions_mol(6)=ions(6)/(water_total*18.01528d-3)  !NO3
      ions_mol(7)=ions(7)/(water_total*18.01528d-3)  !Cl

      !***
      !At this point we may need to introduce a method for prescribing H+ when
      !there is no 'real' value for H+..i.e. in the sulphate poor domain
      !this will give a value for
      !-solve quadratic proposed by Zaveri et al 2005

      !-----------------------------------------------------------------------------------------
      !
      !3) BISULPHATE ION DISSOCIATION
      !
      !Note: the flags "binary_case" and "full_complexity" are not used in this prototype
      !They are used for simplification of the fit expressions when using limited composition regions
      !
      !--this section of code calculates the bisulphate ion concentration

      IF (ions(1) > 0.0 .AND. ions(4) > 0.0) THEN

         !----------HHSO4-------------
         binary_case = 1;

         IF (RH > 0.1 .AND. RH < 0.9) THEN
            binary_hhso4 = -4.9521*(RH**3)+9.2881*(RH**2)-10.777*RH+6.0534;
         ELSE IF (RH >= 0.9 .AND. RH < 0.955) THEN
            binary_hhso4 = -6.3777*RH+5.962;
         ELSE IF (RH >= 0.955 .AND. RH < 0.99) THEN
            binary_hhso4 = 2367.2*(RH**3)-6849.7*(RH**2)+6600.9*RH-2118.7;
         ELSE IF (RH >= 0.99 .AND. RH < 0.9999) THEN
            binary_hhso4 = 3E-07*(RH**5)-2E-05*(RH**4)+0.0004*(RH**3)-0.0035*(RH**2)+0.0123*RH-0.3025;
         END IF

         IF (nitric_acid > 0.0) THEN
            HNO3_hhso4 = -4.2204*(RH**4)+12.193*(RH**3)-12.481*(RH**2)+6.459*RH-1.9004;
         END IF

         IF (hydrochloric_acid > 0.0) THEN
            HCL_hhso4 = -54.845*(RH**7)+209.54*(RH**6)-336.59*(RH**5)+&
               294.21*(RH**4)-150.07*(RH**3)+43.767*(RH**2)-6.5495*RH+0.60048;
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            NH42SO4_hhso4 = 16.768*(RH**3)-28.75*(RH**2)+20.011*RH-8.3206;
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            NH4NO3_hhso4 = -17.184*(RH**4)+56.834*(RH**3)-65.765*(RH**2)+35.321*RH-9.252;
         END IF

         IF (ammonium_chloride > 0.0) THEN
            IF (RH < 0.2 .AND. RH >= 0.1) THEN
               NH4Cl_hhso4 = 3.2809*RH-2.0637;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NH4Cl_hhso4 = -1.2981*(RH**3)+4.7461*(RH**2)-2.3269*RH-1.1259;
            END IF
         END IF

         IF (sodium_sulphate > 0.0) THEN
            Na2SO4_hhso4 = 118.87*(RH**6)-358.63*(RH**5)+435.85*(RH**4)-272.88*(RH**3)+94.411*(RH**2)-18.21*RH+0.45935;
         END IF

         IF (sodium_nitrate > 0.0) THEN
            IF (RH < 0.2 .AND. RH >= 0.1) THEN
               NaNO3_hhso4 = 4.8456*RH-2.5773;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NaNO3_hhso4 = 0.5964*(RH**3)-0.38967*(RH**2)+1.7918*RH-1.9691;
            END IF
         END IF

         IF (sodium_chloride > 0.0) THEN
            IF (RH < 0.2) THEN
               NaCl_hhso4 = 0.51995*RH-1.3981;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NaCl_hhso4 = 1.6539*RH-1.6101;
            END IF
         END IF

         Ln_hhso4_act = binary_hhso4+nitric_acid_eq_frac*HNO3_hhso4+hydrochloric_acid_eq_frac*HCL_hhso4+&
                        ammonium_sulphate_eq_frac*NH42SO4_hhso4+ammonium_nitrate_eq_frac*NH4NO3_hhso4+&
                        ammonium_chloride_eq_frac*NH4Cl_hhso4+sodium_sulphate_eq_frac*Na2SO4_hhso4+&
                        sodium_nitrate_eq_frac*NaNO3_hhso4+sodium_chloride_eq_frac*NaCl_hhso4

         gamma_hhso4 = exp(Ln_hhso4_act)

         !----------H2SO4--------------

         IF (RH >= 0.1 .AND. RH < 0.9) THEN
            binary_h2so4 = 2.4493*(RH**2)-6.2326*RH+2.1763;
         ELSE IF (RH >= 0.9 .AND. RH < 0.98) THEN
            binary_h2so4 = 914.68*(RH**3)-2502.3*(RH**2)+2281.9*RH-695.11;
         ELSE IF (RH >= 0.98 .AND. RH < 0.9999) THEN
            binary_h2so4 = 3E-08*(RH**4)-5E-06*(RH**3)+0.0003*(RH**2)-0.0022*RH-1.1305;
         END IF

         IF (nitric_acid > 0.0) THEN
            HNO3_h2so4 = -16.382*(RH**5)+46.677*(RH**4)-54.149*(RH**3)+34.36*(RH**2)-12.54*RH+2.1368;
         END IF

         IF (hydrochloric_acid > 0.0) THEN
            HCl_h2so4 = -14.409*(RH**5)+42.804*(RH**4)-47.24*(RH**3)+24.668*(RH**2)-5.8015*RH+0.084627;
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            NH42SO4_h2so4 = 66.71*(RH**5)-187.5*(RH**4)+210.57*(RH**3)-121.04*(RH**2)+39.182*RH-8.0606;
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            NH4NO3_h2so4 = -22.532*(RH**4)+66.615*(RH**3)-74.647*(RH**2)+37.638*RH-6.9711;
         END IF

         IF (ammonium_chloride > 0.0) THEN
            IF (RH >= 0.1 .AND. RH < 0.2) THEN
               NH4Cl_h2so4 = -0.32089*RH+0.57738;
            ELSE IF (RH >= 0.2 .AND. RH < 0.9) THEN
               NH4Cl_h2so4 = 18.089*(RH**5)-51.083*(RH**4)+50.32*(RH**3)-17.012*(RH**2)-0.93435*RH+1.0548;
            ELSE IF (RH >= 0.9 .AND. RH < 0.99) THEN
               NH4Cl_h2so4 = -1.5749*RH+1.7002;
            END IF
         END IF

         IF (sodium_sulphate > 0.0) THEN
            Na2SO4_h2so4 = 29.843*(RH**4)-69.417*(RH**3)+61.507*(RH**2)-29.874*RH+7.7556;
         END IF

         IF (sodium_nitrate > 0.0) THEN
            NaNO3_h2so4 = -122.37*(RH**6)+427.43*(RH**5)-604.68*(RH**4)+443.08*(RH**3)-178.61*(RH**2)+37.242*RH-1.9564;
         END IF

         IF (sodium_chloride > 0.0) THEN
            NaCl_h2so4 = -40.288*(RH**5)+115.61*(RH**4)-129.99*(RH**3)+72.652*(RH**2)-22.124*RH+4.2676;
         END IF

         Ln_h2so4_act = binary_h2so4+nitric_acid_eq_frac*HNO3_h2so4+hydrochloric_acid_eq_frac*HCl_h2so4+&
                        ammonium_sulphate_eq_frac*NH42SO4_h2so4+ammonium_nitrate_eq_frac*NH4NO3_h2so4+&
                        ammonium_chloride_eq_frac*NH4Cl_h2so4+sodium_sulphate_eq_frac*Na2SO4_h2so4+&
                        sodium_nitrate_eq_frac*NaNO3_h2so4+sodium_chloride_eq_frac*NaCl_h2so4;

         gamma_h2so4 = exp(Ln_h2so4_act)

	                                           ! export activity coefficients
         IF(gamma_H2SO4 > 0) gamma_out(4) = (gamma_HHSO4**2)/gamma_H2SO4
         IF(gamma_HHSO4 > 0) gamma_out(5) = (gamma_H2SO4**3)/(gamma_HHSO4**2)

         !--solve the quadratic

         act_product = (gamma_h2so4**3)/(gamma_hhso4**2)

         a = 1;
         b = -1*(ions(4)+ions(1)+((water_total*18.0e-3)/(99*act_product)));
         c = ions(4)*ions(1);
         root1 = ( (-1*b) + ( ( sqrt( b**2-4*a*c ) ) ) ) /(2*a)
         root2 = ( (-1*b) - ( ( sqrt( b**2-4*a*c ) ) ) ) /(2*a)

         IF (root1 > ions(1) .OR. root1 < 0.0) THEN
            root1 = 0.0
         END IF

         IF (root2 > ions(1) .OR. root2 < 0.0) THEN
            root2 = 0.0
         END IF

         !--calculate your new hydrogen ion, bisulphate ion and sulphate ion concentration

         hso4_REAL = 0
         h_REAL = ions(1)
         so4_REAL = ions(4)
         IF (root1 == 0.0) THEN
            hso4_REAL = root2
         ELSE IF (root2 == 0.0) THEN
            hso4_REAL = root1
         END IF
         h_REAL = ions(1)-hso4_REAL
         so4_REAL = ions(4)-hso4_REAL

         !--recalculate ion molalities
         ions_mol(1) = h_REAL/(water_total*18.01528e-3) !sulphuric_acid*2+hydrochloric_acid+nitric_acid  !H
         ions_mol(4) = so4_REAL/(water_total*18.01528e-3)  !sulphuric_acid+ammonium_sulphate+sodium_sulphate  !SO4
         ions_mol(5) = hso4_REAL/(water_total*18.01528e-3)  !0.0  !HSO4

         h_out =h_REAL
         hso4_out = hso4_REAL
         so4_out = so4_REAL

      ELSE IF (ions(1) == 0.0 .OR. ions(4) == 0.0) THEN

         h_out = ions(1)
         hso4_out = 0.0
         so4_out = ions(4)

      END IF

      !-----------------------------------------------------------------------------------------
      !
      !4) ACTIVITY COEFFICIENTS -for vapour pressures of HNO3,HCl and NH3
      !
      !Note: the flags "binary_CASE" and "full_complexity" are not used in this prototype
      !They are used for simplification of the fit expressions when using limited composition regions
      !
      !a) - ACTIVITY COEFF/VAPOUR PRESSURE - HNO3
      !
      !-----------binary hno3 act coeff---------------
      IF (ions(1) > 0 .AND. ions(6) > 0) THEN

         binary_case = 1;

         IF (RH > 0.1 .AND. RH < 0.98) THEN
                    !10 to 98! RH full run
            IF (binary_case == 1) THEN
               binary_hno3 = 1.8514*(RH**3)-4.6991*(RH**2)+1.5514*(RH)+0.90236
                    !40 to 98 full
            ELSE IF (binary_case == 2) THEN
               binary_hno3 = -1.1751*(RH**2)-0.53794*(RH)+1.2808
            END IF
         ELSE IF (RH >= 0.98 .AND. RH < 0.9999) THEN
                    !98 to 100!
            binary_hno3 = 1244.69635941351*(RH**3)-2613.93941099991*(RH**2)+1525.0684974546*(RH)-155.946764059316;
         END IF

         !------Contributions from other solutes---------

         full_complexity = 1

         IF (hydrochloric_acid > 0.0) THEN
            !---HCL---
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               HCL_hno3 = 16.051*(RH**4)-44.357*(RH**3)+45.141*(RH**2)-21.638*(RH)+ 4.8182
            !40! +
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               HCL_hno3 = -1.5833*(RH)+1.5569
	                                                    !y = 1.0392*x**{2} - 3.071*x + 2.0431
            END IF
         END IF

         IF (sulphuric_acid > 0.0) THEN
            !---H2SO4----
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               H2SO4_hno3 = -3.0849*(RH**3)+5.9609*(RH**2)-4.468*(RH)+1.5658;
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               !40! +
               H2SO4_hno3 = -0.93473*RH+0.9363
	                                                    !y = - 0.33243*x**{2} - 0.45881*x + 0.78079
            END IF
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            !---NH42SO4----
            NH42SO4_hno3 = 16.821*(RH**3)-28.391*(RH**2)+18.133*(RH)-6.7356
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            !---NH4NO3-----
            NH4NO3_hno3 = 11.01*(RH**3)-21.578*(RH**2)+14.808*RH-4.2593
         END IF

         IF (ammonium_chloride > 0.0) THEN
            !---NH4Cl------
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NH4Cl_hno3 = -1.176*(RH**3)+5.0828*(RH**2)-3.8792*(RH)-0.05518
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               !40! +
               NH4Cl_hno3 = 2.6219*(RH**2)-2.2609*RH-0.38436
            END IF
         END IF

         IF (sodium_sulphate > 0.0) THEN
            !---Na2SO4-----
            Na2SO4_hno3 = 35.504*(RH**4)-80.101*(RH**3)+67.326*(RH**2)-28.461*(RH)+5.6016
         END IF

         IF (sodium_nitrate > 0.0) THEN
            !---NaNO3------
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NaNO3_hno3 = 23.659*(RH**5)-66.917*(RH**4)+74.686*(RH**3)-40.795*(RH**2)+10.831*(RH)-1.4701
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               !40!+
               NaNO3_hno3 = 14.749*(RH**4)-35.237*(RH**3)+31.196*(RH**2)-12.076*(RH)+1.3605
            END IF
         END IF

         IF (sodium_chloride > 0.0) THEN
    	                                                    !---NaCl------
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NaCl_hno3 = 13.682*(RH**4)-35.122*(RH**3)+33.397*(RH**2)-14.586*(RH)+2.6276
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               !40!
               NaCl_hno3 = 1.1882*(RH**3)-1.1037*(RH**2)-0.7642*(RH)+0.6671
            END IF
         END IF

         Ln_HNO3_act = binary_hno3+hydrochloric_acid_eq_frac*HCL_hno3+sulphuric_acid_eq_frac*H2SO4_hno3+&
                       ammonium_sulphate_eq_frac*NH42SO4_hno3+ammonium_nitrate_eq_frac*NH4NO3_hno3+&
                       ammonium_chloride_eq_frac*NH4Cl_hno3+sodium_sulphate_eq_frac*Na2SO4_hno3+&
                       sodium_nitrate_eq_frac*NaNO3_hno3+sodium_chloride_eq_frac*NaCl_hno3

         gamma_hno3 = exp(Ln_HNO3_act)
         gamma_out(1) = gamma_hno3

         !---PARTIAL PRESSURE CALCULATION
         !K_hno3=2.51*(10**6);
         !K_hno3=2.628145923d6  !calculated using AIM online (Clegg et al 1998)
         K_hno3 = 2.6e6*exp(8700*henrys_temp_dep) !after Chameides (1984) (and NIST database)
         Press_HNO3 = (ions_mol(1)*ions_mol(6)*(gamma_hno3**2))/K_hno3

      END IF

      !b) - ACTIVITY COEFF/VAPOUR PRESSURE - NH3

      IF (ions(2) > 0 .AND. ions(1) > 0) THEN
         !follow the two solute approach of zaveri

         !A)-------------------------NH4HSO4-----------------------------
         binary_nh4hso4 = 56.907*(RH**6)-155.32*(RH**5)+142.94*(RH**4)-32.298*(RH**3)-27.936*(RH**2)+19.502*(RH)-4.2618
         !---------binary---------

         IF (nitric_acid > 0.0) THEN
            !---HNO3---
            HNO3_nh4hso4 = 104.8369*(RH**8)-288.8923*(RH**7)+129.3445*(RH**6)+373.0471*(RH**5)-571.0385*(RH**4)+ &
                           326.3528*(RH**3)-74.169*(RH**2)-2.4999*(RH)+3.17
         END IF

         IF (hydrochloric_acid > 0.0) THEN
            !---HCl---
            HCL_nh4hso4 = -7.9133*(RH**8)+126.6648*(RH**7)-460.7425*(RH**6)+731.606*(RH**5)-582.7467*(RH**4)+ &
                           216.7197*(RH**3)-11.3934*(RH**2)-17.7728*(RH)+5.75
         END IF

         IF (sulphuric_acid > 0.0) THEN
            !---H2SO4---
            H2SO4_nh4hso4 = 195.981*(RH**8)-779.2067*(RH**7)+1226.3647*(RH**6)-964.0261*(RH**5)+391.7911*(RH**4)- &
                            84.1409*(RH**3)+20.0602*(RH**2)-10.2663*(RH)+3.5817
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            !---NH42SO4---
            NH42SO4_nh4hso4 = 617.777*(RH**8)-2547.427*(RH**7)+4361.6009*(RH**6)-4003.162*(RH**5)+ &
                              2117.8281*(RH**4)-640.0678*(RH**3)+98.0902*(RH**2)-2.2615*(RH)-2.3811
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            !---NH4NO3---
            NH4NO3_nh4hso4 = -104.4504*(RH**8)+539.5921*(RH**7)-1157.0498*(RH**6)+1322.4507*(RH**5)- &
                              852.2475*(RH**4)+298.3734*(RH**3)-47.0309*(RH**2)+1.297*(RH)-0.8029
         END IF

         IF (ammonium_chloride > 0.0) THEN
            !---NH4Cl---
            NH4Cl_nh4hso4 = 258.1792*(RH**8)-1019.3777*(RH**7)+1592.8918*(RH**6)-1221.0726*(RH**5)+ &
                            442.2548*(RH**4)-43.6278*(RH**3)-7.5282*(RH**2)-3.8459*(RH)+2.2728
         END IF

         IF (sodium_sulphate > 0.0) THEN
            !---Na2SO4---
            Na2SO4_nh4hso4 = 225.4238*(RH**8)-732.4113*(RH**7)+843.7291*(RH**6)-322.7328*(RH**5)- &
                             88.6252*(RH**4)+72.4434*(RH**3)+22.9252*(RH**2)-25.3954*(RH)+4.6971
         END IF

         IF (sodium_nitrate > 0.0) THEN
            !---NaNO3---
            NaNO3_nh4hso4 = 96.1348*(RH**8)-341.6738*(RH**7)+406.5314*(RH**6)-98.5777*(RH**5)- &
                            172.8286*(RH**4)+149.3151*(RH**3)-38.9998*(RH**2)-0.2251*(RH)+0.4953
         END IF

         IF (sodium_chloride > 0.0) THEN
            !---NaCl---
            NaCl_nh4hso4 = 91.7856*(RH**8)-316.6773*(RH**7)+358.2703*(RH**6)-68.9142*(RH**5)- &
                           156.5031*(RH**4)+116.9592*(RH**3)-22.5271*(RH**2)-3.7716*(RH)+1.56
         END IF

         Ln_NH4HSO4_act = binary_nh4hso4+nitric_acid_eq_frac*HNO3_nh4hso4+hydrochloric_acid_eq_frac*HCL_nh4hso4+ &
                          sulphuric_acid_eq_frac*H2SO4_nh4hso4+ammonium_sulphate_eq_frac*NH42SO4_nh4hso4+&
                          ammonium_nitrate_eq_frac*NH4NO3_nh4hso4+ammonium_chloride_eq_frac*NH4Cl_nh4hso4+&
                          sodium_sulphate_eq_frac*Na2SO4_nh4hso4+sodium_nitrate_eq_frac*NaNO3_nh4hso4+&
                          sodium_chloride_eq_frac*NaCl_nh4hso4

         gamma_nh4hso4 = exp(Ln_NH4HSO4_act)
         ! Juha: Put this out as well
         gamma_out(6) = gamma_nh4hso4
                 ! --

         !  -------------------------------------------------------------
         !B)-------------------------HHSO4-------------------------------
         binary_hhso4 = -4.9521*(RH**3)+9.2881*(RH**2)-10.777*RH+6.0534;
         !---------binary---------

         IF (nitric_acid > 0.0) THEN
            !---HNO3---
            HNO3_hhso4 = -4.2204*(RH**4)+12.193*(RH**3)-12.481*(RH**2)+6.459*RH-1.9004;
         END IF

         IF (hydrochloric_acid > 0.0) THEN
            !---HCl---
            HCL_hhso4 = -54.845*(RH**7)+209.54*(RH**6)-336.59*(RH**5)+&
               294.21*(RH**4)-150.07*(RH**3)+43.767*(RH**2)-6.5495*RH+0.60048;
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            !---NH42SO4---
            NH42SO4_hhso4 = 16.768*(RH**3)-28.75*(RH**2)+20.011*RH-8.3206;
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            !---NH4NO3---
            NH4NO3_hhso4 = -17.184*(RH**4)+56.834*(RH**3)-65.765*(RH**2)+35.321*RH-9.252;
         END IF

         IF (ammonium_chloride > 0.0) THEN
            !---NH4Cl---
            IF (RH < 0.2 .AND. RH >= 0.1) THEN
               NH4Cl_hhso4 = 3.2809*RH-2.0637;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NH4Cl_hhso4 = -1.2981*(RH**3)+4.7461*(RH**2)-2.3269*RH-1.1259;
            END IF
         END IF

         IF (sodium_sulphate > 0.0) THEN
            !---Na2SO4---
            Na2SO4_hhso4 = 118.87*(RH**6)-358.63*(RH**5)+435.85*(RH**4)-272.88*(RH**3)+94.411*(RH**2)-18.21*RH+0.45935;
         END IF

         IF (sodium_nitrate > 0.0) THEN
            !---NaNO3---
            IF (RH < 0.2 .AND. RH >= 0.1) THEN
               NaNO3_hhso4 = 4.8456*RH-2.5773;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NaNO3_hhso4 = 0.5964*(RH**3)-0.38967*(RH**2)+1.7918*RH-1.9691;
            END IF
         END IF

         IF (sodium_chloride > 0.0) THEN
            !---NaCl---
            IF (RH < 0.2) THEN
               NaCl_hhso4 = 0.51995*RH-1.3981;
            ELSE IF (RH >= 0.2 .AND. RH < 0.99) THEN
               NaCl_hhso4 = 1.6539*RH-1.6101;
            END IF
         END IF

         Ln_HHSO4_act = binary_hhso4+nitric_acid_eq_frac*HNO3_hhso4+hydrochloric_acid_eq_frac*HCL_hhso4+&
                        ammonium_sulphate_eq_frac*NH42SO4_hhso4+ammonium_nitrate_eq_frac*NH4NO3_hhso4+&
                        ammonium_chloride_eq_frac*NH4Cl_hhso4+sodium_sulphate_eq_frac*Na2SO4_hhso4+&
                        sodium_nitrate_eq_frac*NaNO3_hhso4+sodium_chloride_eq_frac*NaCl_hhso4

         gamma_hhso4 = exp(Ln_HHSO4_act)
         ! Juha: Put this out as well
         gamma_out(7) = gamma_hhso4
                 ! --
         !  -------------------------------------------------------------

         gamma_nh3 = (gamma_nh4hso4**2)/(gamma_hhso4**2)  !exp(Ln_NH3_act)
         gamma_out(3) = gamma_nh3

         !this actually represents the ratio of the ammonium to hydrogen ion activity coefficients
         !(see Zaveri paper) - multiply this by the ratio of the ammonium to hydrogen ion molality
         !and the ratio of appropriate equilibrium constants

	                                           !equilibrium constants
         !Kh=57.64  !Zaveri et al (2005)
         Kh = 5.8e1*exp(4085*henrys_temp_dep) ! after Chameides (1984) (and NIST database)
         !Knh4=1.81e-5  !Zaveri et al (2005)
         Knh4 = 1.7e-5*exp(-4325*henrys_temp_dep) ! after Chameides (1984)
         !Kw=1.01e-14   !Zaveri et al (2005)
         Kw = 1e-14*exp(-6716*henrys_temp_dep) ! after Chameides (1984)

         molality_ratio_nh3 = ions_mol(2)/ions_mol(1)
         Press_NH3 = molality_ratio_nh3*gamma_nh3*(Kw/(Kh*Knh4))

      END IF

   	                          !c) - ACTIVITY COEFF/VAPOUR PRESSURE - HCl

      IF (ions(1) > 0 .AND. ions(7) > 0) THEN
         !-----------binary hcl act coeff---------------
         binary_case = 1;
         IF (RH > 0.1 .AND. RH < 0.98) THEN
                    !10 to 98! RH full run
            IF (binary_case == 1) THEN
               !binary = 1.7526*(RH**2) - 6.9942*(RH) + 5.0046;
               binary_hcl = -5.0179*(RH**3)+9.8816*(RH**2)-10.789*(RH)+5.4737
            !40 to 98 full
            ELSE IF (binary_case == 2) THEN
               binary_hcl = -4.6221*(RH)+4.2633
            END IF
         ELSE IF (RH >= 0.98 .AND. RH < 0.9999) THEN
            !98 to 100!
            binary_hcl = 775.6111008626*(RH**3)-2146.01320888771*(RH**2)+1969.01979670259*(RH)-598.878230033926
         END IF

         IF (nitric_acid > 0.0) THEN
            !---HNO3---
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               HNO3_hcl = 9.6256*(RH**4)-26.507*(RH**3)+27.622*(RH**2)-12.958*(RH)+2.2193
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               HNO3_hcl = 1.3242*(RH**2)-1.8827*(RH)+0.55706
            END IF
         END IF

         IF (sulphuric_acid > 0.0) THEN
            !--h2so4--
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               H2SO4_hcl = 1.4406*(RH**3)-2.7132*(RH**2)+1.014*RH+0.25226
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               H2SO4_hcl = 0.30993*(RH**2)-0.99171*RH+0.66913
            END IF
         END IF

         IF (ammonium_sulphate > 0.0) THEN
            !--NH42SO4---
            NH42SO4_hcl = 22.071*(RH**3)-40.678*(RH**2)+27.893*RH-9.4338
         END IF

         IF (ammonium_nitrate > 0.0) THEN
            !--nh4no3---
            NH4NO3_hcl = 19.935*(RH**3)-42.335*(RH**2)+31.275*(RH)-8.8675
         END IF

         IF (ammonium_chloride > 0.0) THEN
            !--nh4cl---
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NH4Cl_hcl = 2.8048*(RH**3)-4.3182*(RH**2)+3.1971*(RH)-1.6824
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               NH4Cl_hcl = 1.2304*(RH**2)-0.18262*(RH)-1.0643
            END IF
        END IF

         IF (sodium_sulphate > 0.0) THEN
            !--Na2SO4---
            Na2SO4_hcl = 36.104*(RH**4)-78.658*(RH**3)+63.441*(RH**2)-26.727*RH+5.7007
         END IF

         IF (sodium_nitrate > 0.0) THEN
            !--NaNO3---
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NaNO3_hcl = 54.471*(RH**5)-159.42*(RH**4)+180.25*(RH**3)-98.176*(RH**2)+25.309*RH-2.4275
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               NaNO3_hcl = 21.632*(RH**4)-53.088*(RH**3)+47.285*(RH**2)-18.519*RH+2.6846
            END IF
         END IF

         IF (sodium_chloride > 0.0) THEN
            !--NaCl---
            IF (full_complexity == 1 .OR. RH <= 0.4) THEN
               NaCl_hcl = 5.4138*(RH**4)-12.079*(RH**3)+9.627*(RH**2)-3.3164*RH+0.35224
            ELSE IF (full_complexity == 0 .AND. RH > 0.4) THEN
               NaCl_hcl = 2.432*(RH**3)-4.3453*(RH**2)+2.3834*RH-0.4762
            END IF
         END IF

         Ln_HCl_act = binary_hcl+nitric_acid_eq_frac*HNO3_hcl+sulphuric_acid_eq_frac*H2SO4_hcl+&
                      ammonium_sulphate_eq_frac*NH42SO4_hcl+ammonium_nitrate_eq_frac*NH4NO3_hcl+&
                      ammonium_chloride_eq_frac*NH4Cl_hcl+sodium_sulphate_eq_frac*Na2SO4_hcl+&
                      sodium_nitrate_eq_frac*NaNO3_hcl+sodium_chloride_eq_frac*NaCl_hcl

         gamma_hcl = exp(Ln_HCl_act)
         gamma_out(2) = gamma_hcl

	                                           !equilibrium constant
         !K_hcl=1.97d6 !Zaveri et al (2005)
         K_hcl = 2e6*exp(9000.*henrys_temp_dep) ! after Wagman et al (1982) (and NIST database)
         Press_HCl = (ions_mol(1)*ions_mol(7)*(gamma_hcl**2))/K_hcl

      END IF

      ! Juha: for ion malility output
      mols_out = ions_mol

   !REFERENCES
   !Clegg et al (1998) A Thermodynamic Model of the System H+-NH4+-Na+-SO42- -NO3--Cl--H2O at 298.15 K, J. Phys. Chem., 102A, 2155-2171.
   !Clegg et al (2001) Thermodynamic modelling of aqueous aerosols containing electrolytes and dissolved organic compounds. Journal of Aerosol Science 2001;32(6):713-738.
   !Topping et al (2005a) A curved multi-component aerosol hygroscopicity model framework: Part 1 - Inorganic compounds. Atmospheric Chemistry and Physics 2005;5:1205-1222.
   !Topping et al (2005b) A curved multi-component aerosol hygroscopicity model framework: Part 2 - Including organic compounds. Atmospheric Chemistry and Physics 2005;5:1223-1242.
   !Zaveri et al (2005). A new method for multicomponent activity coefficients of electrolytes in aqueous atmospheric aerosols, JGR, 110, 2201, 2005.
   END SUBROUTINE inorganic_pdfite

END MODULE aerosol_thermodynamics


text:
	SUBROUTINE IN_PWND  ( winsym, wintyp, skyflg, condtn, iret )
C************************************************************************
C* IN_PWND								*
C*									*
C* This subroutine decodes the input for a wind barb, a wind arrow, or  *
C* a directional arrow.  The input must be in the form                  *
C*                      : size : width : type : headsiz                 *
C* The arrow/barb size is a multiple of the base size.  Type 1 plots a  *
C* circle or an arrowhead for calm winds.  Type 2 does not plot         *
C* anything for calm winds.  The arrowhead size is a multiple of the    *
C* base arrowhead size.						        *
C* If the user has entered a condition, it must precede the first :     *
C*									*
C* IN_PWND  ( WINSYM, WINTYP, SKYFLG, CONDTN, IRET )			*
C*									*
C* Input parameters:							*
C*	WINSYM		CHAR*		Wind symbol input		*
C*	WINTYP		CHAR*1		Wind type			*
C*					  B = wind barb 		*
C*					  A = wind arrow		*
C*					  D = directional arrow		*
C*	SKYFLG		LOGICAL		Flag set for sky cover plot	*
C*									*
C* Output parameters:							*
C*	CONDTN		CHAR*		Condition			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/NMC	11/91	From IN_WIND				*
C* M. desJardins/NMC	 3/92	Set head size to arrow size if not	*
C*				explicitly set				*
C* D. Kidwell/NCEP	 5/98   Added GSDARR, cleaned up prologue       *
C************************************************************************
	CHARACTER*(*)	winsym, wintyp, condtn
	LOGICAL		skyflg
C*
	REAL		rwind (4)
	CHARACTER	symbol*24
C------------------------------------------------------------------------
	iret  = 0
C
C*	Check for colon.
C
	ipos = INDEX ( winsym, ':' )
	IF  ( ipos .eq. 0 )  THEN
	    condtn = winsym
	    symbol = ' '
	  ELSE
	    condtn = winsym ( 1:ipos-1 )
	    symbol = winsym ( ipos+1:  )
	END IF
C
C*	Decode wind attributes.
C
	CALL ST_RLST  ( symbol, ':', 0., 4, rwind, n, ier )
	size   = rwind (1)
	iwidth = NINT ( rwind (2) )
	itype  = NINT ( rwind (3) )
        sizehd = rwind (4)
	IF  ( sizehd .eq. 0. )  sizehd = size
C
C*	Set calm circle size if required.
C
	ipart = itype / 10
	ipart = MOD ( ipart, 10 )
	IF  ( ipart .eq. 0 )  THEN
	    IF  ( skyflg )  THEN
		itype = itype + 30
	      ELSE
		itype = itype + 10
	    END IF
	END IF
C
C*	Set size, width, type, and/or arrow head size.
C
	IF  ( wintyp .eq. 'B' )  THEN
	    CALL GSBARB  ( size, iwidth, itype, ier )
	  ELSE IF ( wintyp .eq. 'A' )  THEN
	    CALL GSARRW  ( size, sizehd, iwidth, itype, ier )
	  ELSE IF ( wintyp .eq. 'D' )  THEN
	    CALL GSDARR  ( size, sizehd, iwidth, itype, ier )
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! { dg-do compile }
! PR 18883: Fake result variables of non-constant length, in module
module foo
contains
    function s_to_c(chars)
        character, pointer :: chars(:)
        character(len=len(chars)) :: s_to_c
        s_to_c = 'a'
    end function s_to_c
end module foo

program huj

    use foo

    implicit none
    character, pointer :: c(:)
    character(3) :: s

    allocate(c(5))
    c = (/"a", "b", "c" /)
    s = s_to_c(c)

end program huj

! { dg-final { cleanup-modules "foo" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!
! Copyright (C) 2004-2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
!
MODULE atom
  !
  ! ... The variables needed to describe the atoms and related quantities
  !
  USE radial_grids, ONLY : radial_grid_type
  !
  SAVE
  !
  type(radial_grid_type), allocatable, target :: &
       rgrid(:)                ! the information on atomic radial grids.
                               ! NB: some of the subsequent data are therefore redundant
                               ! and will be eliminated in due course asap
  INTEGER, ALLOCATABLE :: &
       msh(:)                  ! the point at rcut
  !
END MODULE atom

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
  do jdw=1,DimDw
     mdw  = Hs(2)%map(jdw)
     ibdw  = bdecomp(mdw,Ns)
     !
     !
     !> H_imp: Off-diagonal elements, i.e. non-local part.
     !remark: iorb=jorb cant have simultaneously n=0 and n=1 (Jcondition)
     do ilat=1,Nlat
         do jlat=1,Nlat
             do iorb=1,Norb
                do jorb=1,Norb
                   is = imp_state_index(ilat,iorb)
                   js = imp_state_index(jlat,jorb)
                   Jcondition = (impHloc(ilat,jlat,Nspin,Nspin,iorb,jorb)/=0d0).AND.(ibdw(js)==1).AND.(ibdw(is)==0)
                   if (Jcondition) then
                      call c(js,mdw,k1,sg1)
                      call cdg(is,k1,k2,sg2)
                      idw = binary_search(Hs(2)%map,k2)
                      htmp = impHloc(ilat,jlat,Nspin,Nspin,iorb,jorb)*sg1*sg2
                      !
                      call sp_insert_element(spH0dws(1),htmp,idw,jdw)
                      !
                   endif
                enddo
             enddo
         enddo
      enddo
   !enddo
     !> H_Bath: inter-orbital bath hopping contribution.
     !if(bath_type=="replica") then
     !   do kp=1,Nbath
     !      do iorb=1,Norb
      !        do jorb=1,Norb
      !           !
       !          ialfa = getBathStride(iorb,kp)
       !          ibeta = getBathStride(jorb,kp)
     !            Jcondition = &
      !                (vca_bath%h(Nspin,Nspin,iorb,jorb,kp)/=zero) &
      !                .AND. (Ndw(ibeta)==1) .AND. (Ndw(ialfa)==0)
                 !
                 !if (Jcondition)then
                 !   call c(ibeta,mdw,k1,sg1)
                 !   call cdg(ialfa,k1,k2,sg2)
                 !   jdw = binary_search(Hs(2)%map,k2)
                 !   htmp = vca_bath%h(Nspin,Nspin,iorb,jorb,kp)*sg1*sg2
                   ! !
                 !   call sp_insert_element(spH0dws(1),htmp,idw,jdw)
                    !
                 !endif
              !enddo
           !enddo
        !enddo
        !
     !endif
     !
     !
     !>H_hyb: hopping terms for a given spin (imp <--> bath)
     if(Nbath>0)then
       do ilat=1,Nlat
         do iorb=1,Norb
            do kp=1,Nbath
               ialfa=getBathStride(ilat,iorb,kp)
               is = imp_state_index(ilat,iorb) !imp site
               if( (diag_hybr(ilat,Nspin,iorb,kp)/=0d0) &
                    .AND. (ibdw(is)==1) .AND. (ibdw(ialfa)==0) )then
                  call c(is,mdw,k1,sg1)
                  call cdg(ialfa,k1,k2,sg2)
                  idw=binary_search(Hs(2)%map,k2)
                  htmp=diag_hybr(ilat,Nspin,iorb,kp)*sg1*sg2
                  !
                  call sp_insert_element(spH0dws(1),htmp,idw,jdw)
                  !
               endif
               if( (diag_hybr(ilat,Nspin,iorb,kp)/=0d0) &
                    .AND. (ibdw(is)==0) .AND. (ibdw(ialfa)==1) )then
                  call c(ialfa,mdw,k1,sg1)
                  call cdg(is,k1,k2,sg2)
                  idw=binary_search(Hs(2)%map,k2)
                  htmp=diag_hybr(ilat,Nspin,iorb,kp)*sg1*sg2
                  !
                  call sp_insert_element(spH0dws(1),htmp,idw,jdw)
                  !
               endif
            enddo
         enddo
       enddo
     endif
     !
   enddo

max_stars_repo_name: QcmPlab/VCA

max_issues_repo_name: QcmPlab/VCA

max_forks_repo_name: QcmPlab/VCA


text:
      SUBROUTINE WASP4
C
C CHANGE RECORD
C **  SUBROUTINE WASPOUT WRITES OUTPUT FILES PROVIDING ADVECTIVE AND
C **  DIFFUSIVE TRANSPORT FIELDS FOR THE WASP4  WATER QUALITY MODEL
C
      USE GLOBAL

      INTEGER,SAVE,ALLOCATABLE,DIMENSION(:)::LDTMP
      INTEGER,SAVE,ALLOCATABLE,DIMENSION(:)::LUTMP
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::QTMP

      IF(.NOT.ALLOCATED(LDTMP))THEN
		ALLOCATE(LDTMP(KCM*LCM))
		ALLOCATE(LUTMP(KCM*LCM))
		ALLOCATE(QTMP(KCM*LCM))
	    LDTMP=0
	    LUTMP=0
	    QTMP=0.0
	ENDIF
C
C **  WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C **  THE VALUE OF X IN THE F10.X FORMATS MAY NEED TO BE CHANGED
C **  FROM PROBLEM TO PROBLEM.  A PRELIMINARY RUN USING E10.3
C **  CAN BE USED TO SPEED THE ADJUSTMENT
C **  READ CONTROL DATA FOR WRITING TO WASP COMPATIABLE FILES
C
      SVPT=1.
      IF(NTSMMT.LT.NTSPTC)SVPT=0.
      IF(JSWASP.EQ.1)THEN
        OPEN(1,FILE='EFDC.WSP',STATUS='UNKNOWN')
        READ(1,1)
        READ(1,1)
        READ(1,*) IVOPT,IBEDV,SCALV,CONVV,VMULT,VEXP,DMULT,DEXP
        READ(1,1)
        READ(1,1)
        READ(1,*) NRFLD,SCALR,CONVR
        READ(1,1)
        READ(1,1)
        READ(1,*) IQOPT,NFIELD,SCALQ,CONVQ
        READ(1,1)
        READ(1,1)
        READ(1,*) DEPSED
        CLOSE(1)
      ENDIF
    1 FORMAT (80X)
C
C **  WRITE HORIZONTAL POSITION AND LAYER FILE WASPP.OUT
C **  WRITE INITIAL VOLUME FILE WASPC.OUT
C **  FILE WASPC.OUT IS CONSISTENT WITH DATA GROUP C SPECIFICATIONS
C **  ON PAGE 172 OF THE WASP4 MANUAL PB88-185095, JAN 1988
C **  FILE WASPP.OUT DEFINES THE LAYER (1 IS SURFACE WATER LAYER, WITH
C **  LAYER NUMBERING INCREASING WITH DEPTH IN WATER COLUMN) AND
C **  HORIZONTAL POSITIONS IN LON,LAT OR UTME, UTMN OF THE WATER
C **  QUALITY (LONG TERM TRANSPORT) CELLS OR SEGEMENTS
C
      IF(JSWASP.EQ.1)THEN
        OPEN(90,FILE='WASPP.OUT',STATUS='UNKNOWN')
        OPEN(93,FILE='WASPC.OUT',STATUS='UNKNOWN')
        CLOSE(90,STATUS='DELETE')
        CLOSE(93,STATUS='DELETE')
        OPEN(90,FILE='WASPP.OUT',STATUS='UNKNOWN')
        OPEN(93,FILE='WASPC.OUT',STATUS='UNKNOWN')
C
C       IVOPT=2
C       IBEDV=0
C
        WRITE(93,1031)IVOPT,IBEDV
C
C       SCALV=1.
C       CONVV=1.
C
        WRITE(93,1032)SCALV,CONVV
C
C       VMULT=0.
C       VEXP=0.
C       DMULT=0.
C       DEXP=0.
C
        LCLTM2=LCLT-2
        LWASP=0
        IF(KC.GT.1)THEN
          LTYPE=1
          KWASP=1
          DO LT=2,LALT
            LWASP=LWASP+1
            LBELOW=LWASP+LCLTM2
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            VOLUME=DXYP(L)*HLPF(L)*DZC(KC)
            WRITE(90,1001)LWASP,KWASP,I,J,DLON(L),DLAT(L)
            WRITE(93,1033)LWASP,LBELOW,LTYPE,VOLUME,VMULT,VEXP,
     &          DMULT,DEXP
          ENDDO
          LTYPE=2
          DO K=KS,2,-1
            KWASP=KC-K+1
            DO LT=2,LALT
              LWASP=LWASP+1
              LBELOW=LWASP+LCLTM2
              I=ILLT(LT)
              J=JLLT(LT)
              L=LIJ(I,J)
              VOLUME=DXYP(L)*HLPF(L)*DZC(K)
              WRITE(90,1001)LWASP,KWASP,I,J,DLON(L),DLAT(L)
              WRITE(93,1033)LWASP,LBELOW,LTYPE,VOLUME,VMULT,VEXP,
     &            DMULT,DEXP
            ENDDO
          ENDDO
        ENDIF
        LTYPE=2
        IF(KC.EQ.1) LTYPE=1
        KWASP=KC
        DO LT=2,LALT
          LWASP=LWASP+1
C
C        LBELOW=0
C
          LBELOW=LWASP+LCLTM2
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          VOLUME=DXYP(L)*HLPF(L)*DZC(1)
          WRITE(90,1001)LWASP,KWASP,I,J,DLON(L),DLAT(L)
          WRITE(93,1033)LWASP,LBELOW,LTYPE,VOLUME,VMULT,VEXP,
     &        DMULT,DEXP
        ENDDO
        LTYPE=3
        KWASP=KC+1
        DO LT=2,LALT
          LWASP=LWASP+1
          LBELOW=0
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          VOLUME=DXYP(L)*DEPSED
          WRITE(90,1001)LWASP,KWASP,I,J,DLON(L),DLAT(L)
          WRITE(93,1033)LWASP,LBELOW,LTYPE,VOLUME,VMULT,VEXP,
     &        DMULT,DEXP
        ENDDO
        CLOSE(90)
        CLOSE(93)
      ENDIF
 1001 FORMAT(4I5,2F10.4)
 1031 FORMAT(2I5)
 1032 FORMAT(2F10.4)
 1033 FORMAT(3I10,5E10.3)
C
C **  WRITE DIFFUSIVE AND DISPERSIVE TRANSPORT FILE WASPB.OUT
C **  FILE WASPB.OUT IS CONSISTENT WITH DATA GROUP B SPECIFICATIONS
C **  ON PAGE 170 OF THE WASP4 MANUAL PB88-185095, JAN 1988
C
      IF(JSWASP.EQ.1)THEN
        OPEN(91,FILE='WASPB.OUT',STATUS='UNKNOWN')
        CLOSE(91,STATUS='DELETE')
        OPEN(91,FILE='WASPB.OUT',STATUS='UNKNOWN')
C
C       NRFLD=1
C
        WRITE(91,1011)NRFLD
        NTEX=NTS/NTSMMT
C
C       SCALR=1.
C       CONVR=1.
C
        WRITE(91,1012)NTEX,SCALR,CONVR
        CLOSE(91)
      ENDIF
      OPEN(91,FILE='WASPB.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      LCLTM2=LCLT-2
      NORSH=0
      NORSV=0
      DO LT=2,LALT
        I=ILLT(LT)
        J=JLLT(LT)
        L=LIJ(I,J)
        NORSH=NORSH+INT(SUB(L))+INT(SVB(L))
        NORSV=NORSV+INT(SPB(L))
      ENDDO
      NORS=KC*NORSH+KS*NORSV
      WRITE(91,1013)NORS
      UNITY=1.
      DO K=KC,1,-1
        KMUL=KC-K
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SUB(L).EQ.1.)THEN
            LWASP=LT-1+KMUL*LCLTM2
            LWASPW=LWASP-1
            LW=LWEST(L)
            ADDLW=DYU(L)*AHULPF(L,K)*DZC(K)*0.5*(HLPF(L)
     &          +HLPF(LW))*DXIU(L)
            WRITE(91,1014)ADDLW,UNITY,LWASPW,LWASP
          ENDIF
        ENDDO
      ENDDO
      UNITY=1.
      DO K=KC,1,-1
        KMUL=KC-K
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SVB(L).EQ.1.)THEN
            LWASP=LT-1+KMUL*LCLTM2
            LSLT=LSCLT(LT)
            LWASPS=LSLT-1+KMUL*LCLTM2
            LS=LSC(L)
            ADDLS=DXV(L)*AHVLPF(L,K)*DZC(K)*0.5*(HLPF(L)
     &          +HLPF(LS))*DYIV(L)
            WRITE(91,1014)ADDLS,UNITY,LWASPS,LWASP
          ENDIF
        ENDDO
      ENDDO
      IF(KC.GT.1)THEN
        UNITY=1.
        DO K=KS,1,-1
          KMUL1=KS-K
          KMUL2=KMUL1+1
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            IF(SPB(L).EQ.1.)THEN
              LWASP=LT-1+KMUL1*LCLTM2
              LBELOW=LT-1+KMUL2*LCLTM2
              ADDL=DXYP(L)*ABLPF(L,K)*DZIG(K)
              WRITE(91,1014)ADDL,UNITY,LWASP,LBELOW
            ENDIF
          ENDDO
        ENDDO
      ENDIF
      NBRK=6
      WRITE(91,1015)NBRK
      IF(ISDYNSTP.EQ.0)THEN
        TSTOP=DT*FLOAT(N)+TCON*TBEGIN
        TSTART=TSTOP-DT*FLOAT(NTSMMT)
      ELSE
        TSTOP=TENDRNSEC
        TSTART=TSTOP-DT*FLOAT(NTSMMT)
      ENDIF
      TSTOP=TSTOP/86400.
      TSTART=TSTART/86400.
      TSMALL=1.E-5
      D1=0.
      T1=0.-2*TSMALL
      D2=0.
      T2=TSTART-TSMALL
      D3=1.
      T3=TSTART+TSMALL
      D4=1.
      T4=TSTOP-TSMALL
      D5=0.
      T5=TSTOP+TSMALL
      D6=0.
      T6=2*TSMALL+(DT*FLOAT(NTS)+TBEGIN*TCON)/86400.
      WRITE(91,1016)D1,T1,D2,T2,D3,T3,D4,T4
      WRITE(91,1016)D5,T5,D6,T6
      CLOSE(91)
 1011 FORMAT(I5)
 1012 FORMAT(I5,2F10.4)
 1013 FORMAT(I5)
 1014 FORMAT(2E10.3,2I5)
 1015 FORMAT(I5)
 1016 FORMAT(4(2F10.5))
 1017 FORMAT(16I5)
C
C **  WRITE ADVECTIVE TRANSPORT FILE WASPD.OUT
C **  FILE WASPD.OUT IS CONSISTENT WITH DATA GROUP D.1 SPECIFICATIONS
C **  ON PAGE 174 OF THE WASP4 MANUAL PB88-185095, JAN 1988
C
      IF(JSWASP.EQ.1)THEN
        OPEN(92,FILE='WASPD.OUT',STATUS='UNKNOWN')
        CLOSE(92,STATUS='DELETE')
        OPEN(92,FILE='WASPD.OUT',STATUS='UNKNOWN')
C
C       IQOPT=1
C       NFIELD=1
C
        WRITE(92,1021)IQOPT,NFIELD
        NINQ=NTS/NTSMMT
C
C       SCALQ=1
C       CONVQ=1
C
        WRITE(92,1022)NINQ,SCALQ,CONVQ
        CLOSE(92)
      ENDIF
      OPEN(92,FILE='WASPD.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      LCLTM2=LCLT-2
      NOQSH=0
      NOQSV=0
      DO LT=2,LALT
        I=ILLT(LT)
        J=JLLT(LT)
        L=LIJ(I,J)
        NOQSH=NOQSH+INT(SUB(L))+INT(SVB(L))
        NOQSV=NOQSV+INT(SWB(L))
      ENDDO
      NOQS=KC*NOQSH+KS*NOQSV
      WRITE(92,1023)NOQS
      LL=0
      DO K=KC,1,-1
        KMUL=KC-K
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SUB(L).EQ.1.)THEN
            LL=LL+1
            LDTMP(LL)=LT-1+KMUL*LCLTM2
            LUTMP(LL)=LDTMP(LL)-1
            QTMP(LL)=DYU(L)*(UHLPF(L,K)+SVPT*UVPT(L,K))*DZC(K)
          ENDIF
        ENDDO
      ENDDO
      DO L=1,LL,4
        WRITE(92,1024)QTMP(L),  LUTMP(L),  LDTMP(L),
     &      QTMP(LEAST(L)),LUTMP(LEAST(L)),LDTMP(LEAST(L)),
     &      QTMP(L+2),LUTMP(L+2),LDTMP(L+2),
     &      QTMP(L+3),LUTMP(L+3),LDTMP(L+3)
      ENDDO
      LL=0
      DO K=KC,1,-1
        KMUL=KC-K
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SVB(L).EQ.1.)THEN
            LL=LL+1
            LSLT=LSCLT(LT)
            LDTMP(LL)=LT-1+KMUL*LCLTM2
            LUTMP(LL)=LSLT-1+KMUL*LCLTM2
            QTMP(LL)=DXV(L)*(VHLPF(L,K)+SVPT*VVPT(L,K))*DZC(K)
          ENDIF
        ENDDO
      ENDDO
      DO L=1,LL,4
        WRITE(92,1024)QTMP(L),  LUTMP(L),  LDTMP(L),
     &      QTMP(LEAST(L)),LUTMP(LEAST(L)),LDTMP(LEAST(L)),
     &      QTMP(L+2),LUTMP(L+2),LDTMP(L+2),
     &      QTMP(L+3),LUTMP(L+3),LDTMP(L+3)
      ENDDO
      IF(KC.GT.1)THEN
        LL=0
        DO K=KS,1,-1
          KMUL1=KS-K
          KMUL2=KMUL1+1
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            IF(SWB(L).EQ.1.)THEN
              LL=LL+1
              LUTMP(LL)=LT-1+KMUL1*LCLTM2
              LDTMP(LL)=LT-1+KMUL2*LCLTM2
              QTMP(LL)=-DXYP(L)*(WLPF(L,K)+SVPT*WVPT(L,K))
            ENDIF
          ENDDO
        ENDDO
      ENDIF
      IF(KC.GT.1)THEN
        DO L=1,LL,4
          WRITE(92,1024) QTMP(L),  LUTMP(L),  LDTMP(L),
     &        QTMP(LEAST(L)),LUTMP(LEAST(L)),LDTMP(LEAST(L)),
     &        QTMP(L+2),LUTMP(L+2),LDTMP(L+2),
     &        QTMP(L+3),LUTMP(L+3),LDTMP(L+3)
        ENDDO
      ENDIF
      NBRKQ=6
      WRITE(92,1025)NBRKQ
      WRITE(92,1026)D1,T1,D2,T2,D3,T3,D4,T4
      WRITE(92,1026)D5,T5,D6,T6
      CLOSE(92)
 1021 FORMAT(2I5)
 1022 FORMAT(I5,2F10.4)
 1023 FORMAT(I5)
 1024 FORMAT(4(E10.3,2I5))
 1025 FORMAT(I5)
 1026 FORMAT(4(2F10.5))
C
C **  WRITE TO DYNHYD.HYD EMULATION FILES WASPDH.OUT AND WASPDHU.OUT
C
      IF(JSWASP.EQ.1)THEN
        OPEN(90,FILE='WASPDHD.OUT',STATUS='UNKNOWN')
        OPEN(94,FILE='WASPDH.OUT',STATUS='UNKNOWN')
        OPEN(95,FILE='WASPDHU.OUT',STATUS='UNKNOWN',
     &      FORM='UNFORMATTED')
        CLOSE(90,STATUS='DELETE')
        CLOSE(94,STATUS='DELETE')
        CLOSE(95,STATUS='DELETE')
        OPEN(90,FILE='WASPDHD.OUT',STATUS='UNKNOWN')
        OPEN(94,FILE='WASPDH.OUT',STATUS='UNKNOWN')
        OPEN(95,FILE='WASPDHU.OUT',STATUS='UNKNOWN',
     &      FORM='UNFORMATTED')
        KCLC=KC*LCLT
        LCLTM2=LCLT-2
        DO KL=1,KCLC
          NCHNC(KL)=0
        ENDDO
        DO M=1,10
          DO KL=1,KCLC
            LCHNC(KL,M)=0
          ENDDO
        ENDDO
        NJUN=KC*(LCLT-2)
        NCHNH=0
        NCHNV=0
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          NCHNH=NCHNH+INT(SUB(L))+INT(SVB(L))
          NCHNV=NCHNV+INT(SWB(L))
        ENDDO
        NCHN=KC*NCHNH+KS*NCHNV
        ISTMP=0
        NODYN=NFLTMT
        TZERO=TBEGIN*TCON/86400.
        WRITE(90,901)NJUN,NCHN
        WRITE(94,941)NJUN,NCHN,DT,ISTMP,NTS,ISTMP,NODYN,TZERO
        WRITE(95)NJUN,NCHN,DT,ISTMP,NTS,ISTMP,NODYN,TZERO
C
C **  CHANNEL DATA
C
        RMNDUM=0.
        LCHN=0
        DO K=KC,1,-1
          KMUL=KC-K
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            IF(SUB(L).EQ.1.)THEN
              LDTM=LT-1+KMUL*LCLTM2
              LUTM=LDTM-1
              RLENTH=DXU(L)
              WIDTH=DYU(L)
              LCHN=LCHN+1
              NCHNC(LDTM)=NCHNC(LDTM)+1
              NCHNC(LUTM)=NCHNC(LUTM)+1
              LCHNC(LDTM,NCHNC(LDTM))=LCHN
              LCHNC(LUTM,NCHNC(LUTM))=LCHN
              WRITE(90,902)LCHN,RLENTH,WIDTH,RMNDUM,LUTM,LDTM
              WRITE(94,942)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
              WRITE(95)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
            ENDIF
          ENDDO
        ENDDO
        DO K=KC,1,-1
          KMUL=KC-K
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            IF(SVB(L).EQ.1.)THEN
              LSLT=LSCLT(LT)
              LDTM=LT-1+KMUL*LCLTM2
              LUTM=LSLT-1+KMUL*LCLTM2
              RLENTH=DYV(L)
              WIDTH=DXV(L)
              LCHN=LCHN+1
              NCHNC(LDTM)=NCHNC(LDTM)+1
              NCHNC(LUTM)=NCHNC(LUTM)+1
              LCHNC(LDTM,NCHNC(LDTM))=LCHN
              LCHNC(LUTM,NCHNC(LUTM))=LCHN
              WRITE(90,902)LCHN,RLENTH,WIDTH,RMNDUM,LUTM,LDTM
              WRITE(94,942)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
              WRITE(95)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
            ENDIF
          ENDDO
        ENDDO
        IF(KC.GT.1)THEN
          DO K=KS,1,-1
            KMUL1=KS-K
            KMUL2=KMUL1+1
            DO LT=2,LALT
              I=ILLT(LT)
              J=JLLT(LT)
              L=LIJ(I,J)
              IF(SWB(L).EQ.1.)THEN
                LUTM=LT-1+KMUL1*LCLTM2
                LDTM=LT-1+KMUL2*LCLTM2
                RLENTH=HLPF(L)*DZG(K)
                WIDTH=SQRT(DXYP(L))
                LCHN=LCHN+1
                NCHNC(LDTM)=NCHNC(LDTM)+1
                NCHNC(LUTM)=NCHNC(LUTM)+1
                LCHNC(LDTM,NCHNC(LDTM))=LCHN
                LCHNC(LUTM,NCHNC(LUTM))=LCHN
                WRITE(90,902)LCHN,RLENTH,WIDTH,RMNDUM,LUTM,LDTM
                WRITE(94,942)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
                WRITE(95)RLENTH,WIDTH,RMNDUM,LUTM,LDTM
              ENDIF
            ENDDO
          ENDDO
        ENDIF
C
C **  JUNCTION DATA
C
        DO K=KC,1,-1
          KMUL=KC-K
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            LCELL=LT-1+KMUL*LCLTM2
            WRITE(90,904)LCELL,DXYP(L),(LCHNC(LCELL,M),M=1,10)
            WRITE(94,944)DXYP(L),(LCHNC(LCELL,M),M=1,10)
            WRITE(95)DXYP(L),(LCHNC(LCELL,M),M=1,10)
          ENDDO
        ENDDO
        CLOSE(90)
        CLOSE(94)
        CLOSE(95)
      ENDIF
C
C **  WRITE TIME STEP, VOLUME AND FLOW DATA
C
      OPEN(94,FILE='WASPDH.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      OPEN(95,FILE='WASPDHU.OUT',POSITION='APPEND',STATUS='UNKNOWN',
     &    FORM='UNFORMATTED')
      LCLTM2=LCLT-2
      IZERO=0
      RZERO=0
      NSTEP=N-NTSMMT
      WRITE(94,945)NSTEP
      DO K=KC,1,-1
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          LN=LNC(L)
          VOLUM=DXYP(L)*HLPF(L)*DZC(K)
          QIN=QSUMELPF(L)*DZC(K)
          FLOWXI=DYU(L)*(UHLPF(L,K)+SVPT*UVPT(L,K))*DZC(K)
          FLOWYI=DXV(L)*(VHLPF(L,K)+SVPT*VVPT(L,K))*DZC(K)
          FLOWZI=DXYP(L)*(WLPF(L,K-1)+SVPT*WVPT(L,K-1))
          FLOWXO=DYU(LEAST(L))*(UHLPF(LEAST(L),K)+SVPT*UVPT(LEAST(L),K))*DZC(K)
          FLOWYO=DXV(LN)*(VHLPF(LN,K)+SVPT*VVPT(LN,K))*DZC(K)
          FLOWZO=DXYP(L)*(WLPF(L,K)+SVPT*WVPT(L,K))
          QQSUM=QIN+FLOWXI+FLOWYI+FLOWZI-FLOWXO-FLOWYO-FLOWZO
          DEPTH=HLPF(L)*DZC(K)
          VELX=0.5*(UHLPF(L,K)+SVPT*UVPT(L,K)
     &        +UHLPF(LEAST(L),K)+SVPT*UVPT(LEAST(L),K))/HLPF(L)
          VELY=0.5*(VHLPF(L,K)+SVPT*VVPT(L,K)
     &        +VHLPF(LN,K)+SVPT*VVPT(LN,K))/HLPF(L)
          VELZ=0.5*(WLPF(L,K-1)+SVPT*WVPT(L,K-1)
     &        +WLPF(L,K)+SVPT*WVPT(L,K))
          VELMAG=SQRT(VELX*VELX+VELY*VELY+VELZ*VELZ)
          WRITE(94,946)VOLUM,QIN,QSUM,DEPTH,VELMAG
          WRITE(95)VOLUM,QIN,QQSUM,DEPTH,VELMAG
        ENDDO
      ENDDO
      DO K=KC,1,-1
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SUB(L).EQ.1.)THEN
            FLOWX=DYU(L)*(UHLPF(L,K)+SVPT*UVPT(L,K))*DZC(K)
            WRITE(94,946)FLOWX
            WRITE(95)FLOWX
          ENDIF
        ENDDO
      ENDDO
      DO K=KC,1,-1
        DO LT=2,LALT
          I=ILLT(LT)
          J=JLLT(LT)
          L=LIJ(I,J)
          IF(SVB(L).EQ.1.)THEN
            FLOWY=DXV(L)*(VHLPF(L,K)+SVPT*VVPT(L,K))*DZC(K)
            WRITE(94,946)FLOWY
            WRITE(95)FLOWY
          ENDIF
        ENDDO
      ENDDO
      IF(KC.GT.1)THEN
        DO K=KS,1,-1
          DO LT=2,LALT
            I=ILLT(LT)
            J=JLLT(LT)
            L=LIJ(I,J)
            IF(SWB(L).EQ.1.)THEN
              FLOWZ=-DXYP(L)*(WLPF(L,K)+SVPT*WVPT(L,K))
              WRITE(94,946)FLOWZ
              WRITE(95)FLOWZ
            ENDIF
          ENDDO
        ENDDO
      ENDIF
      CLOSE(94)
      CLOSE(95)
  901 FORMAT(2I5,E12.4,4I5,E12.4)
  902 FORMAT(I5,2X,3E12.4,2I5)
  903 FORMAT(3E12.4,2I5)
  904 FORMAT(I5,2X,E12.4,10I5)
  905 FORMAT(I5)
  906 FORMAT(5E12.4)
  941 FORMAT(2I5,E12.4,4I5,E12.4)
  942 FORMAT(3E12.4,2I5)
  943 FORMAT(3E12.4,2I5)
  944 FORMAT(E12.4,10I5)
  945 FORMAT(I5)
  946 FORMAT(5E12.4)
      JSWASP=0
      RETURN
      END


text:
*DECK CTANH
      COMPLEX FUNCTION CTANH (Z)
C***BEGIN PROLOGUE  CTANH
C***PURPOSE  Compute the complex hyperbolic tangent.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4C
C***TYPE      COMPLEX (CTANH-C)
C***KEYWORDS  ELEMENTARY FUNCTIONS, FNLIB, HYPERBOLIC TANGENT
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C CTANH(Z) calculates the complex hyperbolic tangent of complex
C argument Z.  Z is in units of radians.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CTAN
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  CTANH
      COMPLEX Z, CI, CTAN
      SAVE CI
      DATA CI /(0.,1.)/
C***FIRST EXECUTABLE STATEMENT  CTANH
      CTANH = -CI*CTAN(CI*Z)
C
      RETURN
      END


text:
C (C) Copyright 1996-2016 ECMWF.
C
C This software is licensed under the terms of the Apache Licence Version 2.0
C which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation nor
C does it submit to any jurisdiction.
C

         PROGRAM METOPS_MSL_AND_PRECIP

      PARAMETER (NLEV=8)
      DIMENSION  RLEV (NLEV)
      DATA       RLEV /0., 1., 4., 8., 20., 50., 100., 200./
      CHARACTER*25 CTAB
      DIMENSION  CTAB  (NLEV)
      DATA       CTAB  /'NONE',
     +                  'GREEN',
     +                  'BLUE',
     +                  'HSL(270,1,0.75)',
     +                  'HSL(300,1,0.75)',
     +                  'HSL(330,1,0.75)',
     +                  'HSL(0,1,0.75)',
     +                  'HSL(0,1,0.75)'/

C     Open MAGICS and set the output device

      CALL POPEN
      CALL PARSE_COMMAND_LINE ('metops_msl_and_precip')

C     Area specification (SOUTH, WEST, NORTH, EAST )

C      CALL PSETC ('SUBPAGE_MAP_PROJECTION', 'POLAR_STEREOGRAPHIC')
      CALL PSETR ('SUBPAGE_LOWER_LEFT_LATITUDE',    35.0)
      CALL PSETR ('SUBPAGE_LOWER_LEFT_LONGITUDE',  -33.0)
      CALL PSETR ('SUBPAGE_UPPER_RIGHT_LATITUDE',   51.0)
      CALL PSETR ('SUBPAGE_UPPER_RIGHT_LONGITUDE',  72.0)
      CALL PSETR ('SUBPAGE_UPPER_RIGHT_LONGITUDE',  0.0)

C     First, load and plot the temperature data, shaded

C     Pass the data to MAGICS

      CALL PSETC ('GRIB_INPUT_TYPE',      'FILE')
      CALL PSETC ('GRIB_INPUT_FILE_NAME', 'data/total_precip.grib')
      CALL PGRIB

C     Set up the coastline attributes

      CALL PSETC ('MAP_COASTLINE_COLOUR', 'BLACK')
      CALL PSETC ('MAP_GRID_COLOUR',      'BLACK')

C     Define the contour

      CALL PSETC  ('LEGEND',                       'ON')
      CALL PSETC  ('CONTOUR',                      'ON')
      CALL PSETC  ('CONTOUR_HIGHLIGHT',            'OFF')
      CALL PSETC  ('CONTOUR_LEVEL_SELECTION_TYPE', 'LEVEL_LIST')
      CALL PSET1R ('CONTOUR_LEVEL_LIST',            RLEV, 2)
      CALL PSETC  ('CONTOUR_SHADE_COLOUR_METHOD',  'LIST')
      CALL PSET1C ('CONTOUR_SHADE_COLOUR_LIST',     CTAB, 2)
      CALL PSETC  ('CONTOUR_SHADE',                'ON')
      CALL PSETC  ('CONTOUR_SHADE_TECHNIQUE',      'POLYGON_SHADING')
      CALL PSETC  ('CONTOUR_SHADE_METHOD',         'AREA_FILL')
      CALL PSETC  ('CONTOUR_HILO',                 'OFF')
      CALL PSETC  ('CONTOUR_LABEL',                'ON')
      CALL PCONT

C     Plot the title text and the coastlines

      CALL PSETR  ('TEXT_REFERENCE_CHARACTER_HEIGHT',  0.35)
      CALL PTEXT
      CALL PCOAST

C     Now load the MSL data

      CALL PSETC ('GRIB_INPUT_FILE_NAME', 'data/msl.grib')
c     CALL PGRIB

C     Redefine the contouring parameters for the MSL

      CALL PSETC  ('CONTOUR',                      'ON')
      CALL PSETR  ('CONTOUR_INTERVAL',              5.0)
      CALL PSETC  ('CONTOUR_LEVEL_SELECTION_TYPE', 'INTERVAL')
      CALL PSETC  ('CONTOUR_LINE_COLOUR',          'BLACK')
      CALL PSETI  ('CONTOUR_LINE_THICKNESS',        2)
      CALL PSETC  ('CONTOUR_HIGHLIGHT_COLOUR',     'BLACK')
      CALL PSETI  ('CONTOUR_HIGHLIGHT_THICKNESS',   4)
      CALL PSETC  ('CONTOUR_LABEL_COLOUR',         'BLACK')
      CALL PSETC  ('CONTOUR_SHADE',                'OFF')
      CALL PSETC  ('CONTOUR_HILO',                 'ON')
      CALL PSETR  ('CONTOUR_HILO_SUPPRESS_RADIUS',  40.0)
      CALL PSETR  ('CONTOUR_HILO_HEIGHT',           0.25)
      CALL PSETC  ('CONTOUR_HI_COLOUR',            'BLACK')
      CALL PSETC  ('CONTOUR_LO_COLOUR',            'BLACK')
      CALL PSETC  ('CONTOUR_LABEL',                'ON')
c      CALL PCONT

      CALL PCLOSE

      STOP
      END

C --------------------------------------------------------------------
C     PARSE_COMMAND_LINE
C     Checks the command-line for any arguments.
C     Arguments can come in pairs. Currently supported arguments are:
C     PROJECTION
C     DEVICE
C     e.g. Run the program with:
C      PROJECTION POLAR_STEREOGRAPHIC
C      PROJECTION CYLINDRICAL   DEVICE SVG
C --------------------------------------------------------------------

      SUBROUTINE PARSE_COMMAND_LINE (OUTROOTNAME)

      CHARACTER*32 ARG
      CHARACTER*64 ID_TEXT
      CHARACTER*32 PROJECTION
      CHARACTER*32 DEVICE
      CHARACTER*48 EXENAME
      CHARACTER*64 OUTNAME
      CHARACTER*(*) OUTROOTNAME
      INTEGER NUM_ARGS
      INTEGER DEVICE_SET

      DEVICE_SET = 0

      ID_TEXT = ''

      NUM_ARGS = IARGC()

      I = 1

20    IF (I.LE.NUM_ARGS) THEN
          CALL GETARG ( I, ARG )

C         Set the projection?

          IF (ARG.EQ.'PROJECTION') THEN
              I = I + 1
              CALL GETARG ( I, PROJECTION )
              CALL PSETC ('SUBPAGE_MAP_PROJECTION', PROJECTION)

C        Set the device?

          ELSEIF (ARG.EQ.'DEVICE') THEN
              I = I + 1
              CALL GETARG ( I, DEVICE )

C             Set the output filename

              IF     (DEVICE.EQ.'PS')  THEN
                OUTNAME = OUTROOTNAME //   '.ps'
                CALL PSETC ('DEVICE',       DEVICE)
                CALL PSETC ('PS_DEVICE',   'ps_a4')
                CALL PSETC ('PS_FILE_NAME', OUTNAME)
                CALL PSETC ('PS_PDF',      'BOTH')
              ELSEIF (DEVICE.EQ.'GIF') THEN
                CALL PSETC ('DEVICE',       'GD')
                OUTNAME = OUTROOTNAME //   '.gif'
                CALL PSETC ('GD_FORMAT',   'ANIMATION')
                CALL PSETC ('GD_FILE_NAME', OUTNAME)
                CALL PSETI ('GD_DELAY',     150)
              ELSEIF (DEVICE.EQ.'PNG') THEN
                CALL PSETC ('DEVICE',       'GD')
                OUTNAME = OUTROOTNAME //   '.png'
                CALL PSETC ('GD_FILE_NAME', OUTNAME)
              ELSEIF (DEVICE.EQ.'SVG') THEN
                CALL PSETC ('DEVICE',       DEVICE)
                OUTNAME = OUTROOTNAME //   '.svg'
                CALL PSETC ('SVG_FILE_NAME', OUTNAME)
              ELSE
                WRITE(*, '(A)') 'BAD DEVICE: ', DEVICE
              ENDIF

            DEVICE_SET = 1

C        Split the PostScript pages into separate files?

          ELSEIF (ARG.EQ.'PS_SPLIT') THEN
                CALL PSETC ('PS_SPLIT',     'ON')

C        Run using linear contouring?

          ELSEIF (ARG.EQ.'LINEAR') THEN
                CALL PSETC ('CONTOUR_METHOD', 'LINEAR')
                CALL PSETC ('PAGE_ID_LINE_USER_TEXT_PLOT', 'ON')
                CALL PSETC ('PAGE_ID_LINE_USER_TEXT', 'LINEAR')
          ENDIF

          I = I + 1
          GOTO 20
      ENDIF

C     If no device has been set, then use PostScript by default

      IF (DEVICE_SET.EQ.0) THEN
        OUTNAME = OUTROOTNAME // '.ps'
        CALL PSETC ('PS_DEVICE',    'ps_a4')
        CALL PSETC ('PS_FILE_NAME', OUTNAME)
      ENDIF

      END

max_stars_repo_name: dtip/magics

max_issues_repo_name: dtip/magics

max_forks_repo_name: dtip/magics


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     This module is used for I/O.
!
!--------------------------------------------------------------------

      MODULE CHNLMOD
      USE ISO_FORTRAN_ENV
      USE UTILMOD
      IMPLICIT NONE

!     Channels cases: standard output, error output, warning output,
!     debugging output
      INTEGER, PARAMETER :: CHNL_O = 601
      INTEGER, PARAMETER :: CHNL_E = 602
      INTEGER, PARAMETER :: CHNL_W = 603
      INTEGER, PARAMETER :: CHNL_D = 604

!     Whether to use color in printing outputs
      LOGICAL :: pClr = .TRUE.

!     A general counter for file ID
      INTEGER :: gFID = 314

!     Appended path to all files that are going to be saved
      CHARACTER(LEN=stdL) :: appPath = ""

!     Channel type, used in I/O
      TYPE chnlType
!        Whether it is open to the screen
         LOGICAL :: oTS = .FALSE.
!        Whether it is open to the file
         LOGICAL :: oTF = .FALSE.
!        Channel identifier
         INTEGER id
!        File ID
         INTEGER fId
!        File address
         CHARACTER(LEN=stdL) :: fName = "histor"
!        Channel tag
         CHARACTER(LEN=stdL) :: tag = ""
      CONTAINS
!        Creates a new channel
         PROCEDURE :: new => newChnl
!        Closes the channel
         PROCEDURE :: close => closeChnl
!        To send a string to channel
         PROCEDURE chnlAssign
         GENERIC :: ASSIGNMENT(=) => chnlAssign
      END TYPE chnlType

!     Only to group four channels, in case I rather have them as one
!     variable
      TYPE ioType
!        Standard output
         TYPE(chnlType) :: o
!        Error
         TYPE(chnlType) :: e
!        Warning
         TYPE(chnlType) :: w
!        Debugging
         TYPE(chnlType) :: d
!        Status file
         TYPE(chnlType) :: s
      CONTAINS
!        Opens all as standard channels
         PROCEDURE :: new => newIO
!        Closes the channel
         PROCEDURE :: close => closeIO
      END TYPE ioType

      INTERFACE ASSIGNMENT(=)
         MODULE PROCEDURE chnlAssignChnl, ioAssignIO
      END INTERFACE ASSIGNMENT(=)

      CONTAINS

!####################################################################

      SUBROUTINE newChnl(chnl,id,fName,tag,oTS,oTF)
      IMPLICIT NONE
      CLASS(chnlType), INTENT(OUT) :: chnl
      INTEGER, INTENT(IN) :: id
      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: fName, tag
      LOGICAL, INTENT(IN), OPTIONAL :: oTS, oTF

      LOGICAL flag

      IF (ALL(id.NE.(/CHNL_O,CHNL_E,CHNL_W,CHNL_D/)))
     2   STOP "Non-existent channel"

      chnl%id = id
      IF (PRESENT(fName)) chnl%fName = ADJUSTL(fName)
      IF (PRESENT(tag)) THEN
         chnl%tag = ADJUSTL(tag)
         IF (chnl%tag .NE. "") chnl%fName =
     2      TRIM(chnl%fName)//"-"//chnl%tag
      END IF

!     By defult outputs are open for error and warning
      IF (PRESENT(oTS)) THEN
         chnl%oTS = oTS
      ELSE
         IF (id.EQ.CHNL_E .OR. id.EQ.CHNL_W) chnl%oTS = .TRUE.
      END IF
      IF (PRESENT(oTF)) THEN
         chnl%oTF = oTF
      ELSE
         IF (id.EQ.CHNL_E .OR. id.EQ.CHNL_W) chnl%oTF = .TRUE.
      END IF
      chnl%fName = TRIM(chnl%fName)//".dat"

      gFID     = gFID + 1
      chnl%fId = gFID
      DO
         INQUIRE(UNIT=chnl%fId, OPENED=flag)
         IF (.NOT.flag) EXIT
         chnl%fId = chnl%fId + 1
      END DO

      RETURN
      END SUBROUTINE newChnl
!--------------------------------------------------------------------
      SUBROUTINE closeChnl(chnl)
      IMPLICIT NONE
      CLASS(chnlType), INTENT(INOUT) :: chnl

      LOGICAL flag

      INQUIRE(UNIT=chnl%fId, OPENED=flag)
      IF (flag) CLOSE(chnl%fId)

      RETURN
      END SUBROUTINE closeChnl
!--------------------------------------------------------------------
      SUBROUTINE chnlAssignChnl(s,r)
      IMPLICIT NONE
      CLASS(chnlType), INTENT(OUT) :: s
      TYPE(chnlType), INTENT(IN) :: r

      s%oTS   = r%oTS
      s%oTF   = r%oTF
      s%id    = r%id
      s%fId   = r%fId
      s%fName = r%fName
      s%tag   = r%tag

      RETURN
      END SUBROUTINE chnlAssignChnl
!--------------------------------------------------------------------
      SUBROUTINE ioAssignIO(s,r)
      IMPLICIT NONE
      CLASS(ioType), INTENT(OUT) :: s
      TYPE(ioType), INTENT(IN) :: r

      s%o = r%o
      s%e = r%e
      s%w = r%w
      s%d = r%d
      s%s = r%s

      RETURN
      END SUBROUTINE ioAssignIO

!####################################################################
!     Sending strings to correct channel
      SUBROUTINE chnlAssign(chnl,sInp)
      IMPLICIT NONE
      CLASS(chnlType), INTENT(INOUT) :: chnl
      CHARACTER(LEN=*), INTENT(IN) :: sInp

      SELECT CASE(chnl%Id)
      CASE (CHNL_O)
         CALL CHNLOUTPUT(chnl,sInp)
      CASE (CHNL_E)
         CALL CHNLERROR(chnl,sInp)
      CASE (CHNL_W)
         CALL CHNLWARNING(chnl,sInp)
      CASE (CHNL_D)
         CALL CHNLDEBUGGING(chnl,sInp)
      END SELECT

      RETURN
      END SUBROUTINE chnlAssign
!--------------------------------------------------------------------
!     This routine is used for keeping track of what is printed on the
!     screen and history file
      SUBROUTINE CHNLOUTPUT(chnl,isTmp)
      IMPLICIT NONE
      TYPE(chnlType), INTENT(INOUT) :: chnl
      CHARACTER(LEN=*), INTENT(IN) :: isTmp

      LOGICAL flag
      CHARACTER(LEN=LEN(isTmp)) sTmp, fName

      IF (.NOT.chnl%oTS .AND. .NOT.chnl%oTF) RETURN

!     Searching and removing all the color codes
      IF (chnl%oTF .OR. (chnl%oTS.AND.(.NOT.pClr))) sTmp = RMCLR(isTmp)

      IF (chnl%oTF) THEN
         INQUIRE(UNIT=chnl%fId, OPENED=flag)
         IF (.NOT.flag) THEN
            IF (appPath .NE. "") CALL SYSTEM("mkdir -p "//TRIM(appPath))
            fName = TRIM(appPath)//TRIM(chnl%fName)
            INQUIRE(FILE=fName, OPENED=flag)
            IF (.NOT.flag) THEN
               OPEN(chnl%fId, FILE=fName, POSITION="APPEND")
            ELSE
               INQUIRE(FILE=fName, NUMBER=chnl%fId)
            END IF
         END IF
         WRITE(chnl%fId,"(A)") TRIM(sTmp)
         CALL FLUSH(chnl%fId)
      END IF

      IF (chnl%oTS) THEN
         IF (chnl%tag .NE. "") THEN
            IF (pClr) THEN
               WRITE(*,"(A)") " "//CLR(TRIM(chnl%tag)//" >>",4)//" "//
     2            TRIM(isTmp)
            ELSE
               WRITE(*,"(A)") " "//TRIM(chnl%tag)//" >> "//TRIM(sTmp)
            END IF
         ELSE
            IF (pClr) THEN
               WRITE(*,"(A)") TRIM(isTmp)
            ELSE
               WRITE(*,"(A)") TRIM(sTmp)
            END IF
         END IF
         CALL FLUSH(OUTPUT_UNIT)
      END IF

      RETURN
      END SUBROUTINE CHNLOUTPUT
!--------------------------------------------------------------------
!     This routine is used for ERROR handling
      SUBROUTINE CHNLERROR(chnl,isTmp)
      IMPLICIT NONE
      TYPE(chnlType), INTENT(INOUT) :: chnl
      CHARACTER(LEN=*), INTENT(IN) :: isTmp

      CHARACTER(LEN=LEN(isTmp)) sTmp

      IF (chnl%oTS .OR. chnl%oTF) THEN
         sTmp = RMCLR(isTmp)

         CALL CHNLOUTPUT(chnl,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"//
     2      "!!!!!!!!")
         CALL CHNLOUTPUT(chnl,"ERROR occurred, see below for more"//
     2      " explanation")
         CALL CHNLOUTPUT(chnl,CLR("ERROR: "//TRIM(ADJUSTL(sTmp))))
         CALL CHNLOUTPUT(chnl,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"//
     2      "!!!!!!!!")
      END IF
      STOP "All processors are forced to stop by a fatal error"

      END SUBROUTINE CHNLERROR
!--------------------------------------------------------------------
!     This routine is used for ERROR handling
      SUBROUTINE CHNLWARNING(chnl,sTmp)
      IMPLICIT NONE
      TYPE(chnlType), INTENT(INOUT) :: chnl
      CHARACTER(LEN=*), INTENT(IN) :: sTmp

      IF (chnl%oTS .OR. chnl%oTF) THEN
         IF (pClr) THEN
            CALL CHNLOUTPUT(chnl,CLR(" WARNING:",4)//" "//
     2         ADJUSTL(sTmp))
         ELSE
            CALL CHNLOUTPUT(chnl," WARNING: "//ADJUSTL(sTmp))
         END IF
      END IF

      RETURN
      END SUBROUTINE CHNLWARNING
!--------------------------------------------------------------------
!     This routine is used for ERROR handling
      SUBROUTINE CHNLDEBUGGING(chnl,sTmp)
      IMPLICIT NONE
      TYPE(chnlType), INTENT(INOUT) :: chnl
      CHARACTER(LEN=*), INTENT(IN) :: sTmp

      IF (chnl%oTS .OR. chnl%oTF) THEN
         CALL CHNLOUTPUT(chnl,CLR(" DEBUG-INFO:",3)//" "//ADJUSTL(sTmp))
      END IF

      RETURN
      END SUBROUTINE CHNLDEBUGGING

!####################################################################

      SUBROUTINE newIO(io, iflag)
      IMPLICIT NONE
      CLASS(ioType), INTENT(OUT) :: io
      LOGICAL, INTENT(IN), OPTIONAL :: iflag

      LOGICAL flag

      IF (PRESENT(iflag)) flag = iflag

      CALL io%o%new(CHNL_O,oTS=flag,oTF=flag)
      CALL io%e%new(CHNL_E)
      CALL io%w%new(CHNL_W)
      CALL io%d%new(CHNL_D)
      CALL io%s%new(CHNL_O,oTF=flag)

      RETURN
      END SUBROUTINE newIO
!--------------------------------------------------------------------
      SUBROUTINE closeIO(io)
      IMPLICIT NONE
      CLASS(ioType), INTENT(INOUT) :: io

      CALL io%o%close()
      CALL io%e%close()
      CALL io%w%close()
      CALL io%d%close()
      CALL io%s%close()

      RETURN
      END SUBROUTINE closeIO

      END MODULE CHNLMOD


text:
Bioletti runs NorthSouth between Hutchison Drive and La Rue.  There is a gate in the middle, so cars arent able to drive from one end to the other.  This is because of the proximity to the Tercero dorms  University officials decided it would be too risky to have everyone driving down Bioletti without slowing down near the dorms.  Everyone who can get past the gate essentially knows about the large number of Freshmen that cross the street, so having the gate there makes it much safer.

It borders the Academic Surge Buildings, Tercero, Kemper Hall, Meyer Hall and Parking Lot 44.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search