MODULE Interp_table
!
!
     USE def_variables
     USE def_constants
     USE properties
     USE interp_functions
     USE non_linear_solvers, ONLY:New_Rap1D
     USE var_const, ONLY : f_out!,xloc, yloc, mess1,mess2

!
     IMPLICIT NONE
!
     PRIVATE
     PUBLIC  CO2BLLT_EQUI !,CO2BLLT_META, CO2BCLT_EQUI, CO2BCLT_META
!
!
     CONTAINS
!
!===============================================================================
!
!SUBROUTINE CO2BLLT_EQUI(p_out, T_out, c_out, x_out, a_out, xxx, u_in,v_in,flag)
SUBROUTINE CO2BLLT_EQUI(p_out, T_out, c_out, x_out, vp_out, xxx, u_in,v_in,flag)
!
!===============================================================================
!
!                       CO2 BILINEAR LOOK-UP TABLES
!
! Input: u_in (the specific internal energy)
!        v_in (the specific volume v_in)
! Output: p_out (pressure)
!         T_out (temerature)
!         c_out (speed of sound)
!         x_out (thermodynamic quality for two-phase regions)
!
! NOTE: TP region tabulated and mix with NewtonRaph to optimise the errors except
!       a small region close to the critical point 27/11/2017.
!       The boundary between LH and HT is optimised. 25/11/2017.
!    f_out = 0(init),1(LL),2(LH),3(R),4(HT),6(LP),7(TPH),8(TPM),9(TPL),10(S),11(S)
!===============================================================================
REAL(pr),INTENT(OUT)  :: p_out, T_out, x_out, c_out, vp_out
INTEGER, INTENT(OUT)  :: flag
REAL(pr),INTENT(IN)   :: xxx !p_guess,res
!
!Local Variables
INTEGER :: i, j, flag_TP, j_sat, i_R, i_L, flag_loca,Niter,exitflag
REAL(pr) :: v_min, v_max, v_sat, v_sat_log, delta,&
&            qual,press,temp,sound,T_guess,out2, a_out
REAL(pr) :: duL_dp,  duV_dp,  dvL_dp,  dvV_dp
REAL(pr) :: vL, uL, vV, uV, pp, ratio, du_dp_x, dv_dp_x, entrov, entrol
REAL(pr) :: x_check, v_check
REAL(pr) :: gamma_pg, e_pg, v_in, u_in,speed
REAL(pr) :: resnorm,out3,T_gs, p_gs, c_gs, u_ori,v_ori,dummy
REAL(pr) :: soundv, soundl
!
!=======================================================================================
!
! 1 step: locating the couple (v,e) in the domain physique
!        - HT: if the input value of internal energy is higher than the value of the
!          maximum point on the saturation curve, the point lies in the HT
!          Region.
!
!          If the input value of internal energy is lower than the value which
!          corresponds to the maximum saturation point, it is necessary to
!          compare the input value of the specific volume with the specific
!          volume of the maximum point of the curve. If the input volume is lower,
!          the input point is in the left part with respect to the maximum, so the
!          loop enters the LL_LH section
!
!        - LL_LH:it is necessary to ensure that the point is inside the defined
!          and tabulated range of validity. The loop stops working if the input
!          point is outside (internal energy lower than the one referred to the
!          triple point, pressure higher than the maximum allowed).
!          If it inside, a comparison between its value of specific internal energy
!          and the critical one, leads to the choice between LL Region and LH Region.

!        - if the point is instead in the right part with respect to the maximum,
!          a similar reasoning is required. It is checked that the point
!          belongs to the tabulated range (its pressure cannot be lower than the
!          Triple Point pressure). Then, through a comparison between values of
!          internal energy and specific volume with respect to the ones at the
!          saturation, R Region or TP Region are chosen respectively.
!
!
!2 step:   choose suitable subroutine to evaluate properties.
!
!==========================================================================================
!flag_TP = 1
flag_loca = 0
x_out   = 1_pr
a_out   = 1_pr
!res     = 0_pr
f_out   = 0
u_ori   = u_in
v_ori   = v_in
!check input
IF ( (u_in /= u_in) .OR. (v_in /= v_in) .OR. (v_in < 0.0) ) THEN
!
!   write(*,*) 'INPUT ERRORS FOR CO2 TABLE from ', mess2,' CALLED by ',mess1
!   write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
    write(*,*) 'INPUT ERRORS FOR CO2 TABLE ','u_in=',u_in,' v_in=',v_in, ' x= ',xxx
   STOP
ENDIF

IF (u_in .GT. e_umax) THEN
!
          IF (u_in .GT. u_end) THEN
!             write(*,*) 'Too high eint for CO2_TABLE from ',mess2,' CALLED by',mess1
!             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
             write(*,*) 'Too high eint for CO2_TABLE ','u_in=',u_in,' v_in=',v_in,' x= ',xxx
             STOP
          ENDIF
!
! To evaluate the interval on the vertical axis, BC of HT
          delta = y_mesh_HT(2) - y_mesh_HT(1)
          i     = INT((u_in    - y_mesh_HT(1))/delta) + 1 !location indice
!
          v_min = 0_pr
          v_max = 0_pr
          DO j = 1, ord_spline + 1
             v_min = v_min + spline_left_HT (ord_spline+2-j,i) * u_in**(j-1)
             v_max = v_max + spline_right_HT(ord_spline+2-j,i) * u_in**(j-1)
          ENDDO
!
          IF (v_in .LT. v_min) THEN
!             write(*,*) 'Pressure higher than 50 MPa in HT, from ',mess2,'CALLED by',mess1
!             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
              write(*,*) 'Pressure higher than 50 MPa in HT ','u_in=',u_in,'v_in=',v_in, ' x= ',xxx
             STOP
          ELSEIF (v_in .GT. v_max) THEN
! In Low pressure (LP)
             flag_loca = 6
          ELSE
! In HT
             flag_loca = 4
          END IF
ELSEIF ((u_in .LE. e_umax) .AND. (v_in .LE. v_umax)) THEN
! Left part:   LL or LH or TP
       IF (u_in .LE. e_cr) THEN
! In LL
           IF (u_in .LT. e_tri_L) THEN
!             write(*,*) 'Too low eint in LL from ',mess2,' CALLED by',mess1
!             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
              write(*,*) 'Too low eint in LL ','u_in=',u_in,' v_in=',v_in, ' x= ',xxx
              STOP
           ENDIF
!
! To evaluate the interval on the vertical axis, BC of LL
          delta  = y_mesh_LL(2) - y_mesh_LL(1)
          i      = INT((u_in    - y_mesh_LL(1))/delta) + 1
!
          v_min = 0_pr;          v_sat = 0_pr
          DO j = 1, ord_spline + 1
             v_min = v_min + spline_pmax_LL(ord_spline+2-j,i) * u_in**(j-1)
             v_sat = v_sat + spline_Lsat_LL(ord_spline+2-j,i) * u_in**(j-1)
!             v_sat_log = v_sat_log + spline_Lsat_LL(ord_spline+2-j,i) *u_in**(j-1)
!             v_sat     = 10_pr ** v_sat_log
          ENDDO
!
          IF (v_in .LT. v_min) THEN
!             write(*,*) 'Pressure higher than 50 MPa in LL from ',mess2,' CALLED by',mess1
!             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
              write(*,*) 'Pressure higher than 50 MPa ','u_in=',u_in,'v_in=',v_in, ' x= ',xxx
             STOP
          ELSEIF (v_in .GT. v_sat) THEN
! Region Two-phase (TPL)
             x_check = (u_in - e_tri_L)/(e_tri_R - e_tri_L)
             v_check = x_check*v_tri_R + (1_pr-x_check)*v_tri_L
             IF (v_in .GT. v_check) THEN
! In Solid, we put it back to TP
                v_in = v_check - 1e-5_pr
                f_out = 10
             ENDIF
             flag_loca = 5
          ELSE
! Region Left Low (LL)
             flag_loca = 1
          ENDIF
! End in LL and begin in LH
       ELSE
!
! To evaluate the interval on the vertical axis, BC of LH
          delta = y_mesh_LH(2)- y_mesh_LH(1)
          i     = INT((u_in - y_mesh_LH(1))/delta) + 1
!
          v_min = 0_pr;          v_sat = 0_pr
          v_sat_log = 0_pr
          DO j = 1, ord_spline + 1
             v_min     = v_min     + spline_pmax_LH(ord_spline+2-j,i) *u_in**(j-1)
             v_sat     = v_sat     + spline_Lsat_LH(ord_spline+2-j,i) *u_in**(j-1)
!             v_sat_log = v_sat_log + spline_Lsat_LH(ord_spline+2-j,i) *u_in**(j-1)
!             v_sat     = 10_pr ** v_sat_log
          ENDDO
!
          IF (v_in .LT. v_min) THEN
!             write(*,*) 'Pressure higher than 50 MPa in LH from ',mess2,' CALLED by',mess1
!             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
             write(*,*) 'Pressure higher than 50 MPa ', 'u_in=',u_in,'v_in=',v_in, ' x= ',xxx
             STOP
          ELSEIF (v_in .GT. v_sat) THEN
! Region Two-phase (TPH + TPM)
             x_check = (u_in - e_tri_L)/(e_tri_R - e_tri_L)
             v_check = x_check*v_tri_R + (1_pr-x_check)*v_tri_L
             IF (v_in .GT. v_check) THEN
! In Solid, we put it back to TP
                v_in = v_check - 1e-5_pr
                f_out = 10
             ENDIF
                flag_loca = 5
          ELSE
! Region Left High (LH)
                flag_loca = 2
          ENDIF
       ENDIF
! Right part: R or TP or LP or solid
ELSE
!
       IF ((u_in .LT. e_tri_R) .AND. (u_in .GT. e_tri_L)) THEN
          x_check = (u_in - e_tri_L)/(e_tri_R - e_tri_L)
          v_check = x_check*v_tri_R + (1_pr-x_check)*v_tri_L
          IF (v_in .GT. v_check) THEN
! In Solid, we put it back to TP
             v_in = v_check - 1e-5_pr
             flag_loca = 5
             f_out = 10
          ENDIF
! In TP
          flag_loca = 5
       ELSEIF (u_in .LT. e_tri_L) THEN
          x_check = (v_in - v_tri_L)/(v_tri_R - v_tri_L)
          u_in    = x_check*e_tri_R + (1_pr-x_check)*e_tri_L
          flag_loca = 5

       ELSE
! To evaluate the interval on the vertical axis, BC of R
          delta = y_mesh_R(2) - y_mesh_R(1)
          i     = INT((u_in   - y_mesh_R(1))/delta) + 1
!
          v_sat = 0_pr;       v_max = 0_pr
          DO j = 1, ord_spline + 1
            v_sat = v_sat + spline_Vsat(ord_spline+2-j,i) * u_in**(j-1)
            v_max = v_max + spline_pmin(ord_spline+2-j,i) * u_in**(j-1)
          ENDDO
!
          IF (v_in .GT. v_max) THEN
! In LP
            flag_loca=6
          ELSEIF (v_in .LT. v_sat) THEN
! In TP
            flag_loca = 5
          ELSE
! In Region Right (R)
            flag_loca = 3
          ENDIF
       ENDIF
!
ENDIF
!
!print*, "location", flag_loca
!
flag = flag_loca
IF (f_out==10) THEN
   write(*,*) 'solid phase !',' u_ori=',u_ori,' v_ori=',v_ori, ' x= ',xxx

!,mess2,' CALLED by',mess1
!   write(*,*) 'u_ori=',u_ori,' v_ori=',v_ori, ' at node, x=',xloc,' y=',yloc
!   v_in = v_tri_R
!   u_in = e_tri_R
!   T_out= T_tri
!   p_out= p_tri
!   c_out= c_tri_R
!   x_out= 1.0
!   a_out= 1.0
   IF (v_ori > v_tri_R) THEN
      v_in  = v_ori
      u_in  = e_tri_R
      flag_loca = 6
      flag = 6
   ENDIF
!   STOP
ENDIF
f_out= flag_loca
!
SELECT CASE (flag_loca)
!
CASE( 0 )
write(*,*)  'Locating the points in the physical domaion failed '
STOP
!
!###### LL #####
CASE( 1 )
     CALL Lin_int_Left_Low(T_out, p_out, c_out, u_in, v_in)

! sound speed correction
!
     IF ( (u_in >-200.0e3_pr) .AND. (v_in >1.9e-3_pr) ) THEN
        CALL sound_speed(T_out,v_in,c_out)
     ENDIF
!
!     CALL entropy(T_out, v_in, entro)
     x_out   = 0_pr
     a_out   = 0_pr
!
!###### LH #####
CASE( 2 )
!
!   IF ( (u_in > e_F) .OR. ((u_in < e_G) .AND. (v_in>v_G)) ) THEN
    IF ( u_in > e_F) THEN
    CALL Lin_int_Left_High(T_gs, p_gs, c_gs, u_in, v_in)
!
    CALL  New_Rap1D(1, T_out, out2, resnorm, Niter,&
     &              exitflag, u_in,T_gs,v_in, out3)
    CALL  pressure(T_out,v_in,p_out)
    CALL  sound_speed(T_out,v_in,c_out)
!
   ELSE
     CALL Lin_int_Left_High(T_out, p_out, c_out, u_in, v_in)
   ENDIF
! speed sound correction
   IF ( (u_in < e_G) .AND. (v_in >v_G) ) THEN
        CALL sound_speed(T_out,v_in,c_out)
   ENDIF
!
   IF ( (p_out > p_cr) .AND. (T_out > T_cr) ) THEN
       x_out   = 0_pr
       a_out   = 0_pr
   ENDIF
!
!
!###### R ######
CASE( 3 )
     delta        = (x_mesh_max - x_mesh_min)/(MMM_R-1)
     x_mesh_R     = x_mesh_min + (/(i*delta, i=0,MMM_R-1)/)
     CALL Lin_int(T_out, p_out, c_out, u_in, v_in, NNN_R, MMM_R, x_mesh_R, y_mesh_R, &
&         spline_pmin, spline_Vsat, TTT_R, ppp_R, ccc_R)

!## HT ##
CASE( 4 )
     delta       = (x_mesh_max - x_mesh_min)/(MMM_HT-1)
     x_mesh_HT   =  x_mesh_min + (/(i*delta, i=0,MMM_HT-1)/)
     CALL Lin_int_Log10(T_out, p_out, c_out, u_in, v_in, NNN_HT, MMM_HT, x_mesh_HT, y_mesh_HT, &
&         spline_right_HT, spline_left_HT, TTT_HT, ppp_HT, ccc_HT)
!
!
   IF ( (p_out > p_cr) .AND. (T_out > T_cr) ) THEN
       x_out   = 0_pr
       a_out   = 0_pr
   ENDIF
!
!
!######## TP #########
CASE( 5 )
!----------------- Spline Tsat and Psat construction flag_TP=0-------------------------
!  IF (flag_TP==0) THEN
!-----------------     ITERATIF PROCESS NOT USE ANY MORE    -------------------------
! initial guess for the saturation pressure (p_guess)
!
!        CALL New_Rap1D(2, press, qual, res, Niter,&
!     &                exitflag, u_in, p_guess, v_in, temp)
!!
!        IF (res > 1e-5) THEN
!          print*, "res IN Interp_table", res, "u,v",u_in,v_in,&
!     &            'T=',temp,'p=',press,'x=',qual
!          STOP '** Interpolation in two-phase region failed in Interp_table case(5)'
!        ENDIF
!p_out = press
!x_out = qual
!T_out = temp
!!
!! SPEED OF SOUND CALCULATION for the HEM:
!!
!       delta = saturP(2) - saturP(1)
!       j_sat = INT((p_out  - saturP(1))/delta) + 1
!!
!       duL_dp = 0_pr;  duV_dp = 0_pr;  dvL_dp = 0_pr;  dvV_dp = 0_pr
!       DO i = 1, ord_spline
!          pp      = p_out**(ord_spline - i)
!          duL_dp  = duL_dp + (ord_spline+1-i) * uL_psat_spline(i,j_sat) *pp
!          duV_dp  = duV_dp + (ord_spline+1-i) * uV_psat_spline(i,j_sat) *pp
!          dvL_dp  = dvL_dp + (ord_spline+1-i) * vL_psat_spline(i,j_sat) *pp
!          dvV_dp  = dvV_dp + (ord_spline+1-i) * vV_psat_spline(i,j_sat) *pp
!       ENDDO
!!
!!
!       vL = 0_pr; uL = 0_pr;  vV = 0_pr;  uV = 0_pr
!       DO i = 1, ord_spline+1
!          pp = p_out**(i-1)
!          vL = vL + vL_psat_spline(ord_spline+2-i, j_sat) * pp
!          uL = uL + uL_psat_spline(ord_spline+2-i, j_sat) * pp
!          vV = vV + vV_psat_spline(ord_spline+2-i, j_sat) * pp
!          uV = uV + uV_psat_spline(ord_spline+2-i, j_sat) * pp
!       ENDDO
!!
!       ratio = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!       du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!       dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!!
!       c_out = SQRT((p_out  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!!
!       a_out = x_out*vV/v_in
!!
!       IF (((a_out .GT. 1_pr) .AND. (x_out .GT. 8e-1_pr)) .OR. (x_out .GT. a_out)) THEN
!             a_out = x_out
!       ENDIF
!---------------------------------------------------------------------------------------------
!!                     LOOK_UP TABLE TWOPHASE in TPL, TPM and TPH, flag_TP=1
!---------------------------------------------------------------------------------------------
!  ELSEIF (flag_TP==1) THEN
!@@@@@@@@@@ TPH @@@@@@@@@@
   IF (u_in .GE. e_tri_R) THEN
!
         f_out=7
!print*, 'TPH','f_out',f_out,u_in,v_in
         delta       = (x_mesh_max - x_mesh_min)/(MMM_TPH-1)
         x_mesh_TPH  =  x_mesh_min + (/(i*delta, i=0,MMM_TPH-1)/)
         CALL Lin_TP_Log10(T_out, p_out, c_out, x_out, u_in, v_in, NNN_TPH, MMM_TPH,&
&             x_mesh_TPH, y_mesh_TPH,spline_Vsat, spline_left_TPH,&
&             TTT_TPH, ppp_TPH, ccc_TPH, xxx_TPH)
!
!         CALL New_Rap1D(2, press, qual, res, Niter,&
!&                       exitflag, u_in, p_out, v_in, temp)
!print*, p_out
!         p_out = press
!         T_out = temp
!         x_out = qual
        IF (u_in .GE. y_mesh_TPH(NNN_TPH-1)) THEN
           p_out = P_umax
           T_out = T_umax
           c_out = c_umax
        ENDIF
!
        CALL satprop (3, p_out, dummy, vV, vL, uV, uL)
!        CALL satderiv(3, p_out, duL_dp, duV_dp, dvL_dp, dvV_dp)
!!
        x_out   = (v_in - vL) / (vV - vL)
        IF ( (0.998<x_out) .AND. (x_out<1.01) )THEN
         x_out = 1.0
        ENDIF
!        ratio   = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!        du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!        dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!!
!       c_out = SQRT((p_out  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!!
        a_out = x_out*vV/v_in
        vp_out= vV
!!
         IF (((a_out .GT. 1_pr) .AND. (x_out .GT. 8e-1_pr)) .OR. (x_out .GT. a_out)) THEN
            a_out = x_out
         ENDIF

!! speed of sound Wood
        CALL sound_speed(T_out,vV,soundv)
        CALL sound_speed(T_out,vL,soundl)

        c_out = v_in * 1_pr/(a_out*vV/(soundv*soundv) + (1_pr-a_out)*vL/(soundl*soundl))
!! Naka
!        c_out = p_out / (a_out*(a_out/vV) + (1.0-a_out)/vL)
        if (c_out <0.0) then
           write(*,*) 'c < 0 in TPH ', 'u_in=',u_in,'v_in=',v_in, ' x= ',xxx
           STOP
        endif
        c_out = sqrt(c_out)

!@@@@@@@@@@@@@ TPL @@@@@@@@@@@@@
   ELSEIF (u_in .LE. e_cr) THEN
!
         f_out=9
         delta       = (x_mesh_max - x_mesh_min)/(MMM_TPL-1)
         x_mesh_TPL  =  x_mesh_min + (/(i*delta, i=0,MMM_TPL-1)/)
         CALL Lin_TP_Log10(T_out, p_out, c_out, x_out, u_in, v_in, NNN_TPL, MMM_TPL,&
&             x_mesh_TPL,y_mesh_TPL,spline_right_TPL, spline_Lsat_LL,&
&             TTT_TPL, ppp_TPL, ccc_TPL, xxx_TPL)
!
!      IF( (u_in < -199e3_pr) .OR.(v_in > 2.34e-3_pr) )THEN
!         CALL New_Rap1D(2, press, qual, res, Niter,&
!&                       exitflag, u_in, p_out, v_in, temp)
!!
!         IF (res > 1_pr) THEN
!            print*, 'TPL interpolation guess prob', res
!!           print*, 'u,v',u_in,v_in, 'p_iter',press, 'p_tb', p_out
!        ENDIF
!        p_out = press
!!       T_out = temp
!!       x_out = qual
         CALL satprop (3, p_out, dummy, vV, vL, uV, uL)
!        CALL satderiv(3, p_out, duL_dp, duV_dp, dvL_dp, dvV_dp)
!!
         x_out   = abs(v_in - vL) / (vV - vL)
!
!       ratio = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!       du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!!      dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!
!       c_out = SQRT((p_out  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!!
         a_out = x_out*vV/v_in
         vp_out= vV
!!
         IF ( abs(x_out) < 1e-4 ) THEN
            x_out =0.0
            flag=1
         ENDIF
!
         IF (((a_out .GT. 1_pr) .AND. (x_out .GT. 8e-1_pr)) .OR. (x_out .GT. a_out)) THEN
            a_out = x_out
         ENDIF
!! speed of sound Wood
        CALL sound_speed(T_out,vV,soundv)
        CALL sound_speed(T_out,vL,soundl)

        c_out = v_in * 1_pr/(a_out*vV/(soundv*soundv) + (1_pr-a_out)*vL/(soundl*soundl))
! Naka
!        c_out = p_out / (a_out*(a_out/vV) + (1.0-a_out)/vL)
        if (c_out <0.0) then
           write(*,*) 'c < 0 in TPL ', 'u_in=',u_in,'v_in=',v_in, ' x= ',xxx
           STOP
        endif
        c_out = sqrt(c_out)
!      ENDIF
!@@@@@@@@@@@@@ TPM @@@@@@@@@@@@@@@@
   ELSE
!
         f_out=8
         delta       = (x_mesh_max - x_mesh_min)/(MMM_TPM-1)
         x_mesh_TPM  =  x_mesh_min + (/(i*delta, i=0,MMM_TPM-1)/)
         CALL Lin_TP_Log10(T_out, p_out, c_out, x_out, u_in, v_in, NNN_TPM, MMM_TPM,&
&             x_mesh_TPM,y_mesh_TPM,spline_right_TPM, spline_left_TPM,&
&             TTT_TPM, ppp_TPM, ccc_TPM, xxx_TPM)
!
!       IF ( (u_in > -183e3_pr) .OR. (v_in >2.3408e-3_pr) ) THEN
!         CALL New_Rap1D(2, press, qual, res, Niter,&
!     &                exitflag, u_in, p_out, v_in, temp)
!
!         IF (res > 1_pr) THEN
!            print*, 'TPM interpolation guess prob', res
!            print*, 'u,v',u_in,v_in, 'p_iter',press, 'p_tb', p_out
!         ENDIF
!!        p_out = press
!         T_out = temp
!         x_out = qual
!
         CALL satprop (3, p_out, dummy, vV, vL, uV, uL)
!        CALL satderiv(3, p_out, duL_dp, duV_dp, dvL_dp, dvV_dp)
!
         x_out   = (v_in - vL) / (vV - vL)
!
!       ratio = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!       du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!       dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!
!       c_out = SQRT((p_out  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!!
        a_out = x_out*vV/v_in
        vp_out= vV
!
        IF (((a_out .GT. 1_pr) .AND. (x_out .GT. 8e-1_pr)) .OR. (x_out .GT. a_out)) THEN
            a_out = x_out
        ENDIF
!! speed of sound Wood
        CALL sound_speed(T_out,vV,soundv)
        CALL sound_speed(T_out,vL,soundl)

        c_out = v_in * 1_pr/(a_out*vV/(soundv*soundv) + (1_pr-a_out)*vL/(soundl*soundl))
! Naka
!        c_out = p_out / (a_out*(a_out/vV) + (1.0-a_out)/vL)
        if (c_out <0.0) then
           write(*,*) 'c < 0 in TPM ', 'u_in=',u_in,'v_in=',v_in, ' x= ',xxx
           STOP
        endif
        c_out = sqrt(c_out)
!    ENDIF
   ENDIF
! ENDIF
!################## LP ###############  perfect gas OR  span-wagner (iterative)
CASE( 6 )
!
        IF (u_in .LT. -123.74e3_pr) THEN
! In Solid phase, but we push it back to LP
         u_in  = -123.74e3_pr + 1e3_pr
         f_out = 11
         write(*,*) 'Maybe in solid phase in LP'
        ENDIF
! Using directly the EoS of span-wagner for small pressure region
!
!        CALL New_Rap1D(1,temp,out2,res,Niter,&
!                        exitflag,u_in,T_guess,v_in,out2)
!        IF (res > 10e-10) THEN
!          print*, "res IN Interp_table", res, "iter", Niter
!          STOP '** conversion from internal energy to temperature in small pressure  region failed in Interp_table case(6)'
!        ENDIF
!
!        CALL pressure(temp, v_in, press)
!        CALL sound_speed(temp, v_in, sound)
!
!        T_out = temp
!        p_out = press
!        c_out = sound
!
!        IF ((p_out .GE. 0.5_pr) .OR. (p_out .LT. 0_pr)) THEN
!         STOP '** Problem in small pressure region in Interp_table case (6)'
!        ENDIF
!--------------  perfect gas formulation---------------------------------------------
         gamma_pg = 1.313_pr
         e_pg     = 236.0294e3_pr
         p_out    = (gamma_pg-1_pr)*(u_in+e_pg)/v_in
         c_out    = sqrt(gamma_pg*p_out*v_in)
         T_out    = p_out*v_in/R
         x_out    = 1_pr
         a_out    = 1_pr
         IF ((p_out .GE. 0.6e6_pr) .OR. (p_out .LT. 1e-2_pr) .OR.&
&            (c_out .LT. 1e-2_pr ) .OR. (T_out .LT. 1e-2_pr)) THEN
!
             write(*,*) 'Problem Perfect gas formulation in LP'
             write(*,*) 'u_in=',u_in,' v_in=',v_in, ' x= ',xxx
             STOP
         ENDIF
END SELECT
!
!
!######### CHECK #########
! Pressure
  IF ( (p_out /= p_out) .OR. (p_out < 0.0) ) THEN
     write(*,*) 'LOOK_UP TABLE pressure negative or inf',p_out,' f_out=',f_out,'flag=',flag
     write(*,*) 'u_in=',u_in,' v_in=',v_in, ' x= ',xxx
     STOP
  ENDIF
! Temperature
  IF ( (T_out /= T_out) .OR. (T_out < 0.0) ) THEN
     write(*,*) 'LOOK_UP TABLE Temperature negative or inf',T_out,'f_out=',f_out, 'flag=',flag
     write(*,*) 'u_in=',u_in,' v_in=',v_in, ' x= ',xxx
     STOP
  ENDIF
! Speed of sound
  IF ( (c_out /= c_out) .OR. (c_out < 0.0) ) THEN
     write(*,*) 'LOOK_UP TABLE sound speed negative or inf',c_out,'f_out=',f_out, 'flag=',flag
     write(*,*) 'u_in=',u_in,' v_in=',v_in, ' x= ',xxx
     write(*,*) 'output, T,p,c,x',T_out, p_out, c_out, x_out
     write(*,*) '------------------------------'
     STOP
  ENDIF
! Quality and mass fraction
  IF ( (x_out /= x_out) .OR. (x_out < 0.0) .OR. &
&      (a_out /= a_out) .OR. (a_out < 0.0) ) THEN
     write(*,*) 'LOOK_UP TABLE quanlity and mass fraction  negative or inf',x_out,a_out,' f_out=',f_out, 'flag=',flag
     write(*,*) 'u_in=',u_in,' v_in=',v_in, ' x= ',xxx
     STOP
  ENDIF
!
  IF ( (T_out<200.0) .OR. (T_out>400) ) THEN
!       write(*,*) 'LOOK_UP TABLE anomal',' f_out=',f_out, 'flag=',flag,'u_in=',u_in,' v_in=',v_in
!       write(*,*) 'u_in=',u_in,' v_in=',v_in, ' at node, x=',xloc,' y=',yloc
  ENDIF

!
END SUBROUTINE CO2BLLT_EQUI
!
!
END MODULE Interp_table


! @Copyright 2007 Kristjan Haule
module RealBubble
  !  ###############################
  !  # Computing real axis Bubble  #
  !  ###############################
  IMPLICIT NONE
  REAL*8, allocatable :: Ome(:)
  REAL*8, allocatable :: chi0r0(:,:,:)
  INTEGER :: norb, nOme
CONTAINS

  SUBROUTINE RealBubble__Init__()
    use LinLogMesh
    use Qlist
    use greenk
    IMPLICIT NONE

    nOme=nom-zero_ind+1
    allocate( Ome(nOme) )
    Ome(:) = om(zero_ind:)
    norb=cixdm
    ! This is the zero frequency value
    allocate( chi0r0(nQ,norb,norb) )
  END SUBROUTINE RealBubble__Init__

  INTEGER FUNCTION find_index(list, element)
    INTEGER, dimension(:) :: list
    INTEGER, intent(in) :: element
    !
    INTEGER :: i
    do i=1,size(list)
       if (list(i).EQ.element) EXIT
    enddo
    find_index =i
    return
  END FUNCTION find_index

  SUBROUTINE CmpChiQ0_real(fileout)
    ! This routine actaually computes Bubble on real axis using formula
    ! P''(Om)_{q,ab} = \sum_k \int_0^{Om}  G''(x)_{k,ab} G''(x-Om)_{k-q,ba} /pi
    ! here G'' = (G - G^+)/(2*i)  and is a complex function for off-diagonal components
    ! We compute P'(Om) by generalized Kramar's-Kronig relation, i.e.,
    !
    !    1/2[chi(om+i*delta)+chi(om-i*delta)] =-1/pi Principal \int_{-\infty}^{\infty} [chi(x+i*delta)-chi(x-i*delta)]/(2*i)/(om-x)
    !
    use LinLogMesh
    use Qlist
    use Klist
    use greenk
    IMPLICIT NONE
    CHARACTER*200, intent(in) :: fileout
    !
    CHARACTER*10 :: skii
    CHARACTER*210:: FNAME
    INTEGER,PARAMETER :: CHUNK= 4
    INTEGER,PARAMETER :: fh_out= 99
    INTEGER :: Q(3)
    INTEGER :: iQ, Qi, ik, ikq, iOm, im, iom1, iom2, izero, dOm, level, iomfin, iorb, jorb
    REAL*8     :: PI, t1c, t1w, t2c, t2w, time_Q_c, time_Q_w, rf, dsum, dsum2
    COMPLEX*16 :: csum, irhok, irhoq, cf
    COMPLEX*16 :: ImBubO(norb,norb)
    INTEGER, allocatable:: om_idx(:)
    COMPLEX*16, allocatable :: ImBub(:,:,:), Bub(:,:,:)
    REAL*8, allocatable :: dx(:), Imf(:), Ref(:), OmKK(:)
    COMPLEX*16, allocatable :: Imc(:), Rec(:)
    COMPLEX*16, PARAMETER  :: IMAG = (0.0D0,1.0D0)

    PI=ACOS(-1.0D0)

    allocate( ImBub(norb,norb,nOme), Bub(norb,norb,nOme) )
    allocate( om_idx(nom), dx(nom) )

    ! Arrays for Kramars-Kronig
    allocate( Imf(2*nOme-1), Ref(nOme) )
    allocate( Imc(2*nOme-1), Rec(nOme) )
    allocate( OmKK(2*nOme-1) )
    ! Mesh for Kramars-Kronig
    do iOm=1,nOme
       OmKK(nOme+iOm-1) = Ome(iOm)
       OmKK(nOme-iOm+1) = -Ome(iOm)
    enddo

    time_Q_c=0
    time_Q_w=0
    do iQ=1,nQ
       Q = Qp(:3,iQ)
       Qi = k_index(Q)
       !print *, 'Qi=', Qi, Q

       call cputim(t1c)
       call walltim(t1w)

       ImBub(:,:,1) = 0

       do iOm=2,nOme
          level = (iOm-2)/Nd+1                      !# which linear mesh should we use?
          om_idx(:) = idx(level,:)                  !# index for the linear mesh on this level
          izero= nidx(level)/2+1                    !# zero_ind on this level

          iomfin = find_index(om_idx(:nidx(level)), zero_ind+iOm-1) ! finds index such that om(om_idx(iomfin))==Ome(iOm)
          dOm = iomfin - izero   !# om-Om in integer notation is im-dOm

          ! dx for trapezoid integration
          dx(izero)=0.5*(om(om_idx(izero+1))-om(om_idx(izero)))
          do im=izero+1,iomfin-1
             dx(im) = 0.5*(om(om_idx(im+1))-om(om_idx(im-1)))
          enddo
          dx(iomfin) = 0.5*(om(om_idx(iomfin))-om(om_idx(iomfin-1)))

          ImBubO=0.0
          !$OMP  PARALLEL DO SHARED(gk,nkp,norb,k_m_q,Qi,om_idx,dx)&
          !$OMP& PRIVATE(ik,iorb,jorb,im,ikq,iom1,iom2,csum,irhok,irhoq)&
          !$OMP& SCHEDULE(STATIC,CHUNK)&
          !$OMP& REDUCTION(+:ImBubO)
          do ik=1,nkp
             ikq  = k_m_q(ik,Qi)
             do iorb=1,norb
                do jorb=1,norb
                   csum=0.0
                   do im=izero,iomfin
                      iom1 = om_idx(im)
                      iom2 = om_idx(im-dOm)
                      irhok=(gk(iom1,iorb,jorb,ik) -conjg(gk(iom1,jorb,iorb,ik)))/2.0
                      irhoq=(gk(iom2,jorb,iorb,ikq)-conjg(gk(iom2,iorb,jorb,ikq)))/2.0
                      csum = csum - irhok*irhoq*dx(im) !  ImG_k * ImG_{k+q}
                   enddo
                   ImBubO(iorb,jorb) = ImBubO(iorb,jorb) + csum/(nkp*PI)
                enddo
             enddo
          enddo
          !$OMP END PARALLEL DO
          ImBub(:,:,iOm)=ImBubO(:,:)
       enddo

       do iorb=1,norb
          do jorb=1,norb
             dsum=0.
             do iOm=1,nOme
                dsum = dsum + abs(dimag(ImBub(iorb,jorb,iOm)))
             enddo
             ! Check if we need complex Kramars-Kronig
             if (dsum.lt.1e-4) then
                ! creating a function which is defined for both positive and negative frequencies
                Imf(nOme)=0.0
                do iOm=2,nOme
                   Imf(nOme+iOm-1) = dreal(ImBub(iorb,jorb,iOm))
                   Imf(nOme-iOm+1) = -dreal(ImBub(iorb,jorb,iOm))
                enddo
                ! Perform Kramars-Kronig for the real part
                !$OMP  PARALLEL DO SHARED(Imf,nOme,OmKK,Ref) PRIVATE(rf,iOm)
                do iOm=2,nOme
                   CALL kramarskronig(rf, Imf, OmKK, nOme+iOm-1, 2*nOme-1)
                   Ref(iOm) = rf
                enddo
                !$OMP END PARALLEL DO
                !
                ! zero frequency value done here, because it can not be obtained by subroutine
                do iOm=1,nOme-1
                   Imf(iOm) = Imf(iOm)/OmKK(iOm)
                enddo
                do iOm=nOme+1,2*nOme-1
                   Imf(iOm) = Imf(iOm)/OmKK(iOm)
                enddo
                Imf(nOme)=0.5*(Imf(nOme-1)+Imf(nOme+1))
                CALL integrate_trapz(rf, Imf, OmKK, 2*nOme-1)
                Ref(1) = rf/PI
                ! Saving the complex function result into Bub
                do iOm=1,nOme
                   Bub(iorb,jorb,iOm) = Ref(iOm) + IMAG*ImBub(iorb,jorb,iOm)
                enddo
             else
                print *, 'Imaginary Kramars-Kronig', iorb, jorb
                Imc(nOme)=0.0
                do iOm=2,nOme
                   Imc(nOme+iOm-1) = ImBub(iorb,jorb,iOm)
                   Imc(nOme-iOm+1) = -ImBub(iorb,jorb,iOm)
                enddo
                ! Perform Kramars-Kronig for the real part
                !$OMP  PARALLEL DO SHARED(Imc,nOme,OmKK,Rec) PRIVATE(cf,iOm)
                do iOm=2,nOme
                   CALL complex_kramarskronig(cf, Imc, OmKK, nOme+iOm-1, 2*nOme-1)
                   Rec(iOm) = cf
                enddo
                !$OMP END PARALLEL DO
                !
                ! zero frequency value done here, because it can not be obtained by subroutine
                do iOm=1,nOme-1
                   Imc(iOm) = Imc(iOm)/OmKK(iOm)
                enddo
                do iOm=nOme+1,2*nOme-1
                   Imc(iOm) = Imc(iOm)/OmKK(iOm)
                enddo
                Imc(nOme)=0.5*(Imc(nOme-1)+Imc(nOme+1))
                CALL complex_integrate_trapz(cf, Imc, OmKK, 2*nOme-1)
                Rec(1) = cf/PI
                ! Saving the result
                do iOm=1,nOme
                   Bub(iorb,jorb,iOm) = Rec(iOm) + IMAG*ImBub(iorb,jorb,iOm)
                enddo
             endif
          enddo
       enddo

       ! Writing out the result
       WRITE(skii,fmt='(I5)') (iQ-1)
       FNAME=TRIM(fileout)//ADJUSTL(TRIM(skii))
       open(fh_out, FILE=FNAME, STATUS='unknown')
       do iOm=1,nOme
          WRITE(fh_out,'(F12.6)',advance='no') Ome(iOm)
          do iorb=1,norb
             do jorb=1,norb
                WRITE(fh_out,'(F14.6,1x,F14.6,3x)',advance='no') dreal(Bub(iorb,jorb,iOm)),dimag(Bub(iorb,jorb,iOm))
             enddo
          enddo
          WRITE(fh_out,*)
       enddo
       close(fh_out)
       ! Timings
       call cputim(t2c)
       call walltim(t2w)
       time_Q_c=time_Q_c+t2c-t1c
       time_Q_w=time_Q_w+t2w-t1w
       WRITE(*,'(A,I4,A,f10.4,A,f10.4)') 'Q=', iQ, ' time[s]=', t2c-t1c, ' walltime[s]', t2w-t1w
    enddo

    deallocate( OmKK )
    deallocate( Imf, Ref )
    deallocate( Imc, Rec )

    deallocate( om_idx, dx )
    deallocate( ImBub, Bub )
    WRITE(*,'(A,f15.4,A,f15.4)') 'Total-time [in min] to calculate Bubble=', time_Q_c/60., ' Total-walltime', time_Q_w/60.

  end SUBROUTINE CmpChiQ0_real

  SUBROUTINE RealBubble__Destruct__
    DEALLOCATE( Ome , chi0r0 )
  END SUBROUTINE RealBubble__Destruct__
end module RealBubble

program RealAxisBubble
  use LinLogMesh
  use Qlist
  use Klist
  use greenk
  use RealBubble
  IMPLICIT NONE
  CHARACTER*200 :: filemesh, filegkr, fbubbleReal, fileQlist, fileKlist, fileout
  INTEGER :: nargs, j, iom
  INTEGER,PARAMETER :: fhi = 995
  CHARACTER*100, allocatable :: argv(:) ! Command-line arguments

  nargs = iargc()
  if (nargs.LT.1) then
     print *, 'Missing the input file to proceed. Create input file with the following lines'
     print *, 'Ba122.klist    # filename with k-list'
     print *, 'Qlist.dat      # filename with Qlist'
     print *, 'rmesh.dat      # real axis mesh'
     print *, 'G_k1r_         # file with real axis k-dependent Grens function'
     print *, 'G_local1r_     # file with real axis local Grens function'
     print *, 'chi0_real.     # name of the output Bubble on real axis'
     call exit(1)
  endif
  ALLOCATE (argv(nargs))                                 ! wouldn't like to parse
  WRITE(*,'(A,I2)') 'nargs=', nargs
  DO j=1,nargs
     CALL getarg(j, argv(j))
     WRITE(*,'(A,A)') 'argi=', TRIM(argv(j))
  ENDDO

  open(fhi, FILE=argv(1),status='old', ERR=900, form='formatted')
  READ(fhi,*)   fileKlist   ! k-list filename, Ba122.klist
  READ(fhi,*)   fileQlist   ! Q-list
  READ(fhi,*)   filemesh    ! real axis mesh
  READ(fhi,*)   filegkr     ! Green's function G_k on real axis
  READ(fhi,*)
  READ(fhi,*)   fileout     ! output filename, for example  'fchiQ.'
  close(fhi)

  ! Reading real axis mesh
  CALL LinLogMesh__Init__(filemesh)
  ! Generates real axis mesh
  CALL LinLogMesh_Generate()

  ! Reading Q list
  CALL Qlist__Init__(fileQlist)
  CALL Klist__Init__(fileKlist)
  CALL greenk__Init__(filegkr)
  if ( sum(abs(om-oml)).gt.1e-5 ) then
     print *, 'Mesh in greens function and real-mesh are not compatible. Diff=', sum(abs(om-oml))
     do iom=1,noml
        print *, iom, om(iom), oml(iom)
     enddo
  endif
  CALL RealBubble__Init__()

  CALL CmpChiQ0_real(fileout)

  CALL RealBubble__Destruct__()
  CALL greenk__Destruct__()
  CALL LinLogMesh__Destruct__()
  CALL Qlist__Destruct__()
  CALL Klist__Destruct__()

  deallocate( argv )
  STOP
900 print *, 'ERROR opening ',TRIM(argv(1)),' file'
end program RealAxisBubble


text: Sacramento Sacred Harp is a Singing Groups group that meets about every two months to sing old fashioned hymns in wiki:wikipedia:Sacred_Harpstyle four part harmony. For more information, contact Susan Jones or David Wright at mailto(davesue@calweb.com) or give them a call at (916) 7398906.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*DECK DXPNRM
      SUBROUTINE DXPNRM (NU1, NU2, MU1, MU2, PQA, IPQA, IERROR)
C***BEGIN PROLOGUE  DXPNRM
C***SUBSIDIARY
C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
C            This subroutine transforms an array of Legendre functions
C            of the first kind of negative order stored in array PQA
C            into normalized Legendre polynomials stored in array PQA.
C            The original array is destroyed.
C***LIBRARY   SLATEC
C***CATEGORY  C3A2, C9
C***TYPE      DOUBLE PRECISION (XPNRM-S, DXPNRM-D)
C***KEYWORDS  LEGENDRE FUNCTIONS
C***AUTHOR  Smith, John M., (NBS and George Mason University)
C***ROUTINES CALLED  DXADJ
C***REVISION HISTORY  (YYMMDD)
C   820728  DATE WRITTEN
C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
C   901019  Revisions to prologue.  (DWL and WRB)
C   901106  Changed all specific intrinsics to generic.  (WRB)
C           Corrected order of sections in prologue and added TYPE
C           section.  (WRB)
C   920127  Revised PURPOSE section of prologue.  (DWL)
C***END PROLOGUE  DXPNRM
      DOUBLE PRECISION C1,DMU,NU,NU1,NU2,PQA,PROD
      DIMENSION PQA(*),IPQA(*)
C***FIRST EXECUTABLE STATEMENT  DXPNRM
      IERROR=0
      L=(MU2-MU1)+(NU2-NU1+1.5D0)
      MU=MU1
      DMU=MU1
      NU=NU1
C
C         IF MU .GT.NU, NORM P =0.
C
      J=1
  500 IF(DMU.LE.NU) GO TO 505
      PQA(J)=0.D0
      IPQA(J)=0
      J=J+1
      IF(J.GT.L) RETURN
C
C        INCREMENT EITHER MU OR NU AS APPROPRIATE.
C
      IF(MU2.GT.MU1) DMU=DMU+1.D0
      IF(NU2-NU1.GT..5D0) NU=NU+1.D0
      GO TO 500
C
C         TRANSFORM P(-MU,NU,X) INTO NORMALIZED P(MU,NU,X) USING
C              NORM P(MU,NU,X)=
C                 SQRT((NU+.5)*FACTORIAL(NU+MU)/FACTORIAL(NU-MU))
C                              *P(-MU,NU,X)
C
  505 PROD=1.D0
      IPROD=0
      K=2*MU
      IF(K.LE.0) GO TO 520
      DO 510 I=1,K
      PROD=PROD*SQRT(NU+DMU+1.D0-I)
  510 CALL DXADJ(PROD,IPROD,IERROR)
      IF (IERROR.NE.0) RETURN
  520 DO 540 I=J,L
      C1=PROD*SQRT(NU+.5D0)
      PQA(I)=PQA(I)*C1
      IPQA(I)=IPQA(I)+IPROD
      CALL DXADJ(PQA(I),IPQA(I),IERROR)
      IF (IERROR.NE.0) RETURN
      IF(NU2-NU1.GT..5D0) GO TO 530
      IF(DMU.GE.NU) GO TO 525
      PROD=SQRT(NU+DMU+1.D0)*PROD
      IF(NU.GT.DMU) PROD=PROD*SQRT(NU-DMU)
      CALL DXADJ(PROD,IPROD,IERROR)
      IF (IERROR.NE.0) RETURN
      MU=MU+1
      DMU=DMU+1.D0
      GO TO 540
  525 PROD=0.D0
      IPROD=0
      MU=MU+1
      DMU=DMU+1.D0
      GO TO 540
  530 PROD=SQRT(NU+DMU+1.D0)*PROD
      IF(NU.NE.DMU-1.D0) PROD=PROD/SQRT(NU-DMU+1.D0)
      CALL DXADJ(PROD,IPROD,IERROR)
      IF (IERROR.NE.0) RETURN
      NU=NU+1.D0
  540 CONTINUE
      RETURN
      END


SUBROUTINE PrintForces2(forcea,fsph,fsph2,fnsp,fvdrho,forb,fextra)
  !     prints partial forces
  use structure, only : mult, natm
  USE com,   only : nat
  IMPLICIT NONE
  !IMPLICIT REAL*8 (A-H,O-Z)
  logical, intent(in) :: forcea(0:3,nat)
  real*8, intent(in)  :: fsph(3,natm),fsph2(3,natm),fnsp(3,natm),fvdrho(3,nat), forb(3,natm), fextra(3,nat)
  ! local variables
  real*8 :: fpsi(3,nat), fequ(3), ftot(3,nat)
  real*8 :: fequ_nrm, fnsp_nrm, fsph_nrm, fsph2_nrm, forb_nrm, fvdrho_nrm, fpsi_nrm, ftot_nrm, fextra_nrm
  INTEGER :: jatom, ia, ia1, ia2, ik

  fpsi(:,:)=0.d0
  DO jatom=1,nat
     ia2 = sum(mult(1:jatom))
     ia1 = ia2-mult(jatom)+1
     DO ia=ia1,ia2
        fequ(:)=0.d0
        DO ik=1,3
           if (forcea(ik,jatom)) fequ(ik)=(fsph(ik,ia)+fsph2(ik,ia)+forb(ik,ia)+fnsp(ik,ia))/mult(jatom)
        ENDDO
        fpsi(:,jatom) = fpsi(:,jatom) + fequ(:) ! averaged force in local coordiante system
        fequ_nrm = sqrt(dot_product(fequ,fequ))
        fnsp_nrm = sqrt(dot_product(fnsp(:,ia),fnsp(:,ia)))     ! Force2
        fsph_nrm = sqrt(dot_product(fsph(:,ia),fsph(:,ia)))     ! Force1
        fsph2_nrm= sqrt(dot_product(fsph2(:,ia),fsph2(:,ia)))   ! Force3
        forb_nrm = sqrt(dot_product(forb(:,ia),forb(:,ia)))     ! ForceU
        write(6,77) jatom,ia-ia1+1,'+ SPH', fsph_nrm,  (fsph(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ SP2', fsph2_nrm, (fsph2(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ ORB', forb_nrm,  (forb(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ NSP', fnsp_nrm,  (fnsp(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'> EQU', fequ_nrm,  (fequ(ik),ik=1,3)
     ENDDO
     !     calculate partial forces
     ftot(:,jatom) = fvdrho(:,jatom) + fpsi(:,jatom) + fextra(:,jatom)
     fvdrho_nrm    = sqrt(dot_product(fvdrho(:,jatom),fvdrho(:,jatom)))  ! Force4
     fpsi_nrm      = sqrt(dot_product(fpsi(:,jatom),fpsi(:,jatom)))
     fextra_nrm    = sqrt(dot_product(fextra(:,jatom),fextra(:,jatom)))  ! Extra force due to discountinuity of the basis functions
     ftot_nrm      = sqrt(dot_product(ftot(:,jatom),ftot(:,jatom)))
     write(6,77) jatom,1,'= PSI', fpsi_nrm,   (fpsi(ik,jatom),ik=1,3)
     write(6,77) jatom,1,'+ VDR', fvdrho_nrm, (fvdrho(ik,jatom),ik=1,3)
     write(6,77) jatom,1,'+ EXT', fextra_nrm, (fextra(ik,jatom),ik=1,3)
     write(6,77) jatom,1,'> TOT', ftot_nrm,   (ftot(ik,jatom),ik=1,3)
     write(21,78)
     write(21,79) jatom,jatom, ftot_nrm*1000, (ftot(ik,jatom)*1000,ik=1,3)
     write(6,*)

77   format(2i3,a7,4e15.7)
78   FORMAT (7x,'VALENCE-FORCE IN mRy/a.u. = |F|',3x,'Fx',13x,'Fy',13x,'Fz')
79   FORMAT (':FVA',i3.3,':',1x,i3,'.ATOM',4f17.9)
  ENDDO
  RETURN
END SUBROUTINE PrintForces2

SUBROUTINE PrintForces(forcea,fsph,fsph2,fnsp,fvdrho,forb)
  !     prints partial forces
  use structure, only : mult, natm
  USE com,   only : nat
  IMPLICIT NONE
  !IMPLICIT REAL*8 (A-H,O-Z)
  logical, intent(in) :: forcea(0:3,nat)
  real*8, intent(in)  :: fsph(3,natm),fsph2(3,natm),fnsp(3,natm),fvdrho(3,nat), forb(3,natm)
  ! local variables
  real*8 :: fpsi(3,nat), fequ(3), ftot(3,nat)
  real*8 :: fequ_nrm, fnsp_nrm, fsph_nrm, fsph2_nrm, forb_nrm, fvdrho_nrm, fpsi_nrm, ftot_nrm
  INTEGER :: jatom, ia, ia1, ia2, ik

  fpsi(:,:)=0.d0
  DO jatom=1,nat
     ia2 = sum(mult(1:jatom))
     ia1 = ia2-mult(jatom)+1
     DO ia=ia1,ia2
        fequ(:)=0.d0
        DO ik=1,3
           if (forcea(ik,jatom)) fequ(ik)=(fsph(ik,ia)+fsph2(ik,ia)+forb(ik,ia)+fnsp(ik,ia))/mult(jatom)
        ENDDO
        fpsi(:,jatom) = fpsi(:,jatom) + fequ(:) ! averaged force in local coordiante system
        fequ_nrm = sqrt(dot_product(fequ,fequ))
        fnsp_nrm = sqrt(dot_product(fnsp(:,ia),fnsp(:,ia)))     ! Force2
        fsph_nrm = sqrt(dot_product(fsph(:,ia),fsph(:,ia)))     ! Force1
        fsph2_nrm= sqrt(dot_product(fsph2(:,ia),fsph2(:,ia)))   ! Force3
        forb_nrm = sqrt(dot_product(forb(:,ia),forb(:,ia)))     ! ForceU
        write(6,77) jatom,ia-ia1+1,'+ SPH', fsph_nrm,  (fsph(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ SP2', fsph2_nrm, (fsph2(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ ORB', forb_nrm,  (forb(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'+ NSP', fnsp_nrm,  (fnsp(ik,ia),ik=1,3)
        write(6,77) jatom,ia-ia1+1,'> EQU', fequ_nrm,  (fequ(ik),ik=1,3)
     ENDDO
     !     calculate partial forces
     ftot(:,jatom) = fvdrho(:,jatom) + fpsi(:,jatom)
     fvdrho_nrm    = sqrt(dot_product(fvdrho(:,jatom),fvdrho(:,jatom)))  ! Force4
     fpsi_nrm      = sqrt(dot_product(fpsi(:,jatom),fpsi(:,jatom)))
     ftot_nrm      = sqrt(dot_product(ftot(:,jatom),ftot(:,jatom)))
     write(6,77) jatom,1,'= PSI', fpsi_nrm,   (fpsi(ik,jatom),ik=1,3)
     write(6,77) jatom,1,'+ VDR', fvdrho_nrm, (fvdrho(ik,jatom),ik=1,3)
     write(6,77) jatom,1,'> TOT', ftot_nrm,   (ftot(ik,jatom),ik=1,3)
     write(21,78)
     write(21,79) jatom,jatom, ftot_nrm*1000, (ftot(ik,jatom)*1000,ik=1,3)
     write(6,*)

77   format(2i3,a7,4e15.7)
78   FORMAT (7x,'VALENCE-FORCE IN mRy/a.u. = |F|',3x,'Fx',13x,'Fy',13x,'Fz')
79   FORMAT (':FVA',i3.3,':',1x,i3,'.ATOM',4f17.9)
  ENDDO
  RETURN
END SUBROUTINE PrintForces


      subroutine header_sd_channel

      use sd_channel_module
      use basin_module
      use hydrograph_module

      implicit none

!!!  SWAT-DEG CHANNEL
      if (sp_ob%chandeg > 0) then
        if (pco%sd_chan%d == "y") then
          open (2500,file="channel_sd_day.txt",recl = 1500)
          write (2500,*) bsn%name, prog
          write (2500,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
          write (2500,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
          !write (2500,*) ch_wbod_hdr, sd_hyd_hdr
          !write (2500,*) ch_wbod_hdr_units, sd_hyd_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_day.txt"
          if (pco%csvout == "y") then
            open (2504,file="channel_sd_day.csv",recl = 1500)
            write (2504,*) bsn%name, prog
            write (2504,'(*(G0.3,:,","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
            write (2504,'(*(G0.3,:,","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
            !write (2504,'(*(G0.3,:,","))') ch_wbod_hdr, sd_hyd_hdr
            !write (2504,'(*(G0.3,:,","))') ch_wbod_hdr_units, sd_hyd_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_day.csv"
          end if
        endif
      endif

        if (sp_ob%chandeg > 0) then
          if (pco%sd_chan%m == "y") then
          open (2501,file="channel_sd_mon.txt",recl = 1500)
          write (2501,*) bsn%name, prog
          write (2501,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
          write (2501,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
          !write (2501,*) ch_wbod_hdr, sd_hyd_hdr
          !write (2501,*) ch_wbod_hdr_units, sd_hyd_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_mon.txt"
          if (pco%csvout == "y") then
            open (2505,file="channel_sd_mon.csv",recl = 1500)
            write (2505,*) bsn%name, prog
            write (2505,'(*(G0.3,:,","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
            write (2505,'(*(G0.3,:,","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
            !write (2505,'(*(G0.3,:,","))') ch_wbod_hdr, sd_hyd_hdr
            !write (2505,'(*(G0.3,:,","))') ch_wbod_hdr_units, sd_hyd_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_mon.csv"
          end if
          end if
         end if

      if (sp_ob%chandeg > 0) then
        if (pco%sd_chan%y == "y") then
          open (2502,file="channel_sd_yr.txt",recl = 1500)
          write (2502,*) bsn%name, prog
          write (2502,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
          write (2502,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
          !write (2502,*) ch_wbod_hdr, sd_hyd_hdr
          !write (2502,*) ch_wbod_hdr_units, sd_hyd_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_yr.txt"
          if (pco%csvout == "y") then
            open (2506,file="channel_sd_yr.csv",recl = 1500)
            write (2506,*) bsn%name, prog
            write (2506,'(*(G0.3,:,","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
            write (2506,'(*(G0.3,:,","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
            !write (2506,'(*(G0.3,:,","))') ch_wbod_hdr, sd_hyd_hdr
            !write (2506,'(*(G0.3,:,","))') ch_wbod_hdr_units, sd_hyd_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_yr.csv"
          end if
        endif
      endif

        if (sp_ob%chandeg > 0) then
          if (pco%sd_chan%a == "y") then
          open (2503,file="channel_sd_aa.txt",recl = 1500)
          write (2503,*) bsn%name, prog
          write (2503,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
          write (2503,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
          !write (2503,*) ch_wbod_hdr, sd_hyd_hdr
          !write (2503,*) ch_wbod_hdr_units, sd_hyd_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_aa.txt"
          if (pco%csvout == "y") then
            open (2507,file="channel_sd_aa.csv",recl = 1500)
            write (2507,*) bsn%name, prog
            write (2507,'(*(G0.3,:,","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr, wtmp_hdr
            write (2507,'(*(G0.3,:,","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1, wtmp_units
            !write (2507,'(*(G0.3,:,","))') ch_wbod_hdr, sd_hyd_hdr
            !write (2507,'(*(G0.3,:,","))') ch_wbod_hdr_units, sd_hyd_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL          channel_sd_aa.csv"
          end if
          end if
         end if

!!!!!!!! SD_CHANMORPH
      if (sp_ob%chandeg > 0) then
        if (pco%sd_chan%d == "y") then
          open (4800,file="channel_sdmorph_day.txt",recl = 1500)
          write (4800,*) bsn%name, prog
          write (4800,*) sdch_hdr !! swat deg channel morph
          write (4800,*) sdch_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_day.txt"
          if (pco%csvout == "y") then
            open (4804,file="channel_sdmorph_day.csv",recl = 1500)
            write (4804,*) bsn%name, prog
            write (4804,'(*(G0.3,:,","))') sdch_hdr
            write (4804,'(*(G0.3,:,","))') sdch_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_day.csv"
          end if
        endif
      endif

        if (sp_ob%chandeg > 0) then
          if (pco%sd_chan%m == "y") then
          open (4801,file="channel_sdmorph_mon.txt",recl = 1500)
          write (4801,*) bsn%name, prog
          write (4801,*) sdch_hdr   !! swat deg channel morph
          write (4801,*) sdch_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_mon.txt"
          if (pco%csvout == "y") then
            open (4805,file="channel_mon_sdmorph.csv",recl = 1500)
            write (4805,*) bsn%name, prog
            write (4805,'(*(G0.3,:,","))') sdch_hdr
            write (4805,'(*(G0.3,:,","))') sdch_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_mon.csv"
          end if
          end if
         end if

      if (sp_ob%chandeg > 0) then
        if (pco%sd_chan%y == "y") then
          open (4802,file="channel_sdmorph_yr.txt",recl = 1500)
          write (4802,*) bsn%name, prog
          write (4802,*) sdch_hdr !! swat deg channel morph
          write (4802,*) sdch_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_yr.txt"
          if (pco%csvout == "y") then
            open (4806,file="channel_sdmorph_yr.csv",recl = 1500)
            write (4806,*) bsn%name, prog
            write (4806,'(*(G0.3,:,","))') sdch_hdr !! swat deg channel morph csv
            write (4806,'(*(G0.3,:,","))') sdch_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_yr.csv"
          end if
        endif
      endif

        if (sp_ob%chandeg > 0) then
          if (pco%sd_chan%a == "y") then
          open (4803,file="channel_sdmorph_aa.txt",recl = 1500)
          write (4803,*) bsn%name, prog
          write (4803,*) sdch_hdr   !! swat deg channel morph
          write (4803,*) sdch_hdr_units
          write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_aa.txt"
          if (pco%csvout == "y") then
            open (4807,file="channel_sdmorph_aa.csv",recl = 1500)
            write (4807,*) bsn%name, prog
            write (4807,'(*(G0.3,:,","))') sdch_hdr
            write (4807,'(*(G0.3,:,","))') sdch_hdr_units
            write (9000,*) "SWAT-DEG_CHANNEL_MORPH    channel_sdmorph_aa.csv"
          end if
          end if
         end if
!!!!!!!! SD_CHANMORPH

      return
      end subroutine header_sd_channel


text:
program tst

character (8) N

character (8) g

character (7) T

character (7) M

character (5) M1

g='gabinete'

T='teclado'

M='monitor'

M1='mouse'

10 continue

write(*,*) 'Escreva o nome do item'

read(*,*,iostat=ierr)N

if(ierr/=0) then

write(*,*) 'Erro, tente de novo'

go to 10

end if

select case(N)

case('gabinete')

write(*,*) 'Custa 99.90 reais'

go to 10

case('teclado')

write(*,*) 'Custa 20.00 reais'

go to 10

case('monitor')

write(*,*) 'Custa 450.00 reais'

go to 10

case('mouse')

write(*,*) 'Custa 15.00 reais'

go to 10

case default

write(*,*) 'Erro, tente novamente'

go to 10

end select

stop

end program

max_stars_repo_name: danielsanfr/fortran-study

max_issues_repo_name: danielsanfr/fortran-study

max_forks_repo_name: danielsanfr/fortran-study


function test1(n)
  integer  :: n
  character(n) :: test1
  character(n) :: bar1
  test1 = ""
  return
entry bar1()
  bar1 = ""
end function test1

function test2()
  character(1) :: test2
  character(1) :: bar2
  test2 = ""
  return
entry bar2()
  bar2 = ""
end function test2

function test3() ! { dg-warning "Obsolescent feature" }
  character(*) :: test3
  character(*) :: bar3 ! { dg-warning "Obsolescent feature" }
  test3 = ""
  return
entry bar3()
  bar3 = ""
end function test3

function test4(n) ! { dg-warning "returning variables of different string lengths" }
  integer  :: n
  character(n) :: test4
  character(*) :: bar4 ! { dg-warning "Obsolescent feature" }
  test4 = ""
  return
entry bar4()
  bar4 = ""
end function test4

function test5() ! { dg-warning "returning variables of different string lengths" }
  character(1) :: test5
  character(2) :: bar5
  test5 = ""
  return
entry bar5()
  bar5 = ""
end function test5

function test6() ! { dg-warning "Obsolescent feature|returning variables of different string lengths" }
  character(*) :: test6
  character(2) :: bar6
  test6 = ""
  return
entry bar6()
  bar6 = ""
end function test6


text:
      subroutine swtche(id,nodes,a,b,idir,n)
c     implicit double precision (a-h,o-z)                                    dp
      dimension id(2,*),nodes(*),a(*),b(*),idir(*)
      do 10 i=1,n
      ndf=idir(i)
      if (ndf.ne.3000000) then
      b(ndf)=a(ndf)
      a(ndf)=0.
      else
      nodnum=nodes(i)
      ndf1=id(1,nodnum)
      ndf2=id(2,nodnum)
      b(ndf1)=a(ndf1)
      b(ndf2)=a(ndf2)
      a(ndf1)=0.0
      a(ndf2)=0.0
      endif
   10 continue
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
!##############################################################################
!# ****************************************************************************
!# <name> cc2dmedium_method2 </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module is a demonstation program how to solve a stationary
!# Navier-Stokes problem
!#
!#              $$- \nu Laplace(u) + u*grad(u) + \Nabla p = f $$
!#              $$ \Nable \cdot p = 0$$
!#
!# on a 2D domain for a 2D function $u=(u_1,u_2)$ and a pressure $p$.
!#
!# The routine splits up the tasks of reading the domain, creating
!# triangulations, discretisation, solving, postprocessing and creanup into
!# different subroutines. The communication between these subroutines
!# is done using an application-specific structure saving problem data
!# as well as a collection structure for the communication with callback
!# routines.
!#
!# For the nonlinearity, the nonlinear solver is invoked. The
!# defect that is setted up there is preconditioned by a linear Multigrid
!# solver with a simple-VANCA smoother/preconditioner for
!# 2D saddle point problems, Jacobi-Type. As coarse grid solver,
!# UMFPACK is used.
!# </purpose>
!##############################################################################

module cc2dmedium_method2

  use fsystem
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use coarsegridcorrection
  use spdiscprojection
  use nonlinearsolver
  use paramlist
  use statistics

  use collection
  use convection

  use cc2dmediumm2basic
  use cc2dmediumm2init
  use cc2dmediumm2boundary
  use cc2dmediumm2discretisation
  use cc2dmediumm2postprocessing
  use cc2dmediumm2stationary

  implicit none

contains

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

!<subroutine>

  subroutine cc2dmedium2_getLogFiles (slogfile,serrorfile)

!<description>
  ! Temporarily reads the output DAT file to get the names of the output
  ! files.
!</description>

!<output>
  ! Name of the message log file.
  character(LEN=*), intent(OUT) :: slogfile

  ! Name of the error log file.
  character(LEN=*), intent(OUT) :: serrorfile
!</output>

!</subroutine>

    type(t_parlist) :: rparlist
    character(LEN=SYS_STRLEN) :: sstring

    ! Init parameter list that accepts parameters for output files
    call parlst_init (rparlist)

    ! Read parameters that configure the output
    call parlst_readfromfile (rparlist, './data/output.dat')

    ! Now the real initialisation of the output including log file stuff!
    call parlst_getvalue_string (rparlist,'GENERALOUTPUT',&
                                'smsgLog',sstring,'')
    read(sstring,*) slogfile

    call parlst_getvalue_string (rparlist,'GENERALOUTPUT',&
                                'serrorLog',sstring,'')
    read(sstring,*) serrorfile

    ! That temporary parameter list is not needed anymore.
    call parlst_done (rparlist)

    end subroutine

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

!<subroutine>

  subroutine cc2dmedium2_getDAT (rparamList)

!<description>
  ! Reads in all DAT files into the parameter list rparlist
!</description>

!<inputoutput>
  ! The parameter list where the values of the DAT files should be stored.
  ! The structure must have been initialised, the parameters are just added
  ! to the list.
  type(t_parlist), intent(INOUT) :: rparamList
!</inputoutput>

!</subroutine>

    ! Each 'readfromfile' command adds the parameter of the specified file
    ! to the parameter list.
    call parlst_readfromfile (rparamList, './data/discretisation.dat')
    call parlst_readfromfile (rparamList, './data/linsol_cc2d.dat')
    call parlst_readfromfile (rparamList, './data/nonlinsol_cc2d.dat')
    call parlst_readfromfile (rparamList, './data/output.dat')
    call parlst_readfromfile (rparamList, './data/paramtriang.dat')
    call parlst_readfromfile (rparamList, './data/bdconditions.dat')
    call parlst_readfromfile (rparamList, './data/timediscr.dat')
    call parlst_readfromfile (rparamList, './data/postprocessing.dat')

  end subroutine

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

!<subroutine>

  subroutine cc2dmedium2

!<description>
  ! This is a 'separated' Navier-Stokes solver for solving a Navier-Stokes
  ! problem. The different tasks of the problem are separated into
  ! subroutines. The problem uses a problem-specific structure for the
  ! communication: All subroutines add their generated information to the
  ! structure, so that the other subroutines can work with them.
  ! (This is somehow a cleaner implementation than using a collection!).
  ! For the communication to callback routines of black-box subroutines
  ! (matrix-assembly), a collection is used.
  !
  ! The following tasks are performed by the subroutines:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to GMV file
  ! 8.) Release all variables, finish
!</description>

!</subroutine>

    ! A problem structure for our problem
    type(t_problem), pointer :: p_rproblem

    ! A structure for the solution vector and the RHS vector of the problem.
    type(t_vectorBlock) :: rvector,rrhs

    ! A structure for the postprocessing.
    type(t_c2d2postprocessing) :: rpostprocessing

    ! Timer objects for stopping time
    type(t_timer) :: rtimerTotal
    type(t_timer) :: rtimerGridGeneration
    type(t_timer) :: rtimerMatrixGeneration
    type(t_timer) :: rtimerSolver

    integer :: i

    ! Ok, let's start.

    ! Initialise the timers by zero:
    call stat_clearTimer(rtimerTotal)
    call stat_clearTimer(rtimerGridGeneration)
    call stat_clearTimer(rtimerMatrixGeneration)
    call stat_clearTimer(rtimerSolver)

    ! Start the timer
    call stat_startTimer(rtimerTotal)

    ! Allocate memory fo rthe problem structure -- it's rather large!
    allocate (p_rproblem)

    ! Initialise the collection
    call collct_init (p_rproblem%rcollection)
    do i=1,NNLEV
      call collct_addlevel_all (p_rproblem%rcollection)
    end do

    ! Initialise the parameter list object. This creates an empty parameter list.
    call parlst_init (p_rproblem%rparamList)

    ! Add the parameter list to the collection so that the parameters
    ! from the DAT/INI files are available everywhere where we have the
    ! collection.
    call collct_setvalue_parlst(p_rproblem%rcollection,'INI',&
                                p_rproblem%rparamList,.true.)

    ! Read parameters from the INI/DAT files into the parameter list.
    call cc2dmedium2_getDAT (p_rproblem%rparamList)

    ! Ok, parameters are read in.
    ! Get the output levels during the initialisation phase and during the program.
    call c2d2_initOutput (p_rproblem)

    ! Print the configuration to the terminal
    if (p_rproblem%MSHOW_Initialisation .ge. 2) then
      call output_line ('Parameters:')
      call parlst_info (p_rproblem%rparamList)
    end if

    ! Evaluate these parameters and initialise global data in the problem
    ! structure for global access.
    call c2d2_initParameters (p_rproblem)

    ! So now the different steps - one after the other.
    !
    ! Initialisation
    !
    ! Parametrisation & Triangulation
    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising parametrisation / triangulation...')
    end if

    call stat_startTimer(rtimerGridGeneration)

    call c2d2_initParamTriang (p_rproblem)

    call stat_stopTimer(rtimerGridGeneration)
    call output_lbrk ()
    call output_line ("Time for mesh generation: "//&
      trim(sys_sdL(rtimerGridGeneration%delapsedReal,10)))

    ! Print mesh information
    if (p_rproblem%MSHOW_Initialisation .ge. 2) then
      call output_lbrk ()
      call output_line ('Mesh statistics:')
      call output_lbrk ()
      do i=p_rproblem%NLMIN,p_rproblem%NLMAX
        call tria_infoStatistics (p_rproblem%RlevelInfo(i)%rtriangulation,&
            i .eq. p_rproblem%NLMIN,i)
      end do
    end if

    ! Discretisation
    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising discretisation...')
    end if
    call c2d2_initDiscretisation (p_rproblem)

    if (p_rproblem%MSHOW_Initialisation .ge. 2) then
      call output_lbrk ()
      call output_line ('Discretisation statistics:')
      do i=p_rproblem%NLMIN,p_rproblem%NLMAX
        call output_lbrk ()
        call output_line ('Level '//sys_siL(i,5))
        call dof_infoDiscrBlock (p_rproblem%RlevelInfo(i)%p_rdiscretisation,.false.)
      end do
    end if

    ! And all the other stuff...
    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising postprocessing...')
    end if
    call c2d2_initPostprocessing (p_rproblem,rpostprocessing)

    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising matrices/vectors...')
    end if

    call stat_startTimer(rtimerMatrixGeneration)

    call c2d2_allocMatVec (p_rproblem,rvector,rrhs)

    call stat_stopTimer(rtimerMatrixGeneration)
    call output_lbrk ()
    call output_line ("Time for matrix initialisation: "//&
      trim(sys_sdL(rtimerGridGeneration%delapsedReal,10)))

    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising analytic boundary conditions...')
    end if
    call c2d2_initAnalyticBC (p_rproblem)

    ! On all levels, generate the static matrices used as templates
    ! for the system matrix (Laplace, B, Mass,...)
    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Generating basic matrices...')
    end if
    call c2d2_generateBasicMatrices (p_rproblem)

    ! Create the solution vector -- zero or read from file.
    if (p_rproblem%MSHOW_Initialisation .ge. 1) then
      call output_separator (OU_SEP_MINUS)
      call output_line('Initialising initial solution vector...')
    end if
    call c2d2_initInitialSolution (p_rproblem,rvector)

    ! Now choose the algorithm. Stationary or time-dependent simulation?
    if (p_rproblem%itimedependence .eq. 0) then

      ! Stationary simulation
      !
      ! Generate the RHS vector.
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Generating RHS vector...')
      end if
      call c2d2_generateBasicRHS (p_rproblem,rrhs)

      ! Generate discrete boundary conditions
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Generating discrete boundary conditions...')
      end if
      call c2d2_initDiscreteBC (p_rproblem,rvector,rrhs)

      ! Implementation of boundary conditions
      call c2d2_implementBC (p_rproblem,rvector,rrhs,.true.,.true.)

      ! Solve the problem
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Invoking stationary solver...')
        call output_separator (OU_SEP_MINUS)
      end if

      call stat_startTimer(rtimerSolver)

      call c2d2_solve (p_rproblem,rvector,rrhs)

      call stat_stopTimer(rtimerSolver)

      ! Postprocessing
      call c2d2_postprocessingStationary (p_rproblem,rvector,rpostprocessing)

    else

      ! Time dependent simulation with explicit time stepping.
      !
      ! Generate the RHS vector for the first time step.
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Generating RHS vector...')
      end if
      call c2d2_generateBasicRHS (p_rproblem,rrhs)

      ! Initialise the boundary conditions, but
      ! don't implement any boundary conditions as the nonstationary solver
      ! doesn't like this.
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Generating discrete boundary conditions of first time step...')
      end if
      call c2d2_initDiscreteBC (p_rproblem,rvector,rrhs)

      ! Call the nonstationary solver to solve the problem.
      if (p_rproblem%MSHOW_Initialisation .ge. 1) then
        call output_separator (OU_SEP_MINUS)
        call output_line('Invoking nonstationary solver...')
        call output_separator (OU_SEP_MINUS)
      end if

      call stat_startTimer(rtimerSolver)

      call c2d2_solveNonstationary (p_rproblem,rvector,rrhs,rpostprocessing)

      call stat_stopTimer(rtimerSolver)

    end if

    ! (Probably) write final solution vector
    call c2d2_writeSolution (p_rproblem,rvector)

    ! Cleanup
    call c2d2_doneMatVec (p_rproblem,rvector,rrhs)
    call c2d2_doneBC (p_rproblem)
    call c2d2_doneDiscretisation (p_rproblem)
    call c2d2_donepostprocessing (rpostprocessing)
    call c2d2_doneParamTriang (p_rproblem)

    ! Release parameters from the DAT/INI files from the problem structure.
    call c2d2_doneParameters (p_rproblem)

    ! Release the parameter list
    call collct_deleteValue (p_rproblem%rcollection,'INI')
    call parlst_done (p_rproblem%rparamList)

    ! Print some statistical data about the collection - anything forgotten?
    call output_lbrk ()
    call output_line ('Remaining collection statistics:')
    call output_line ('--------------------------------')
    call output_lbrk ()
    call collct_printStatistics (p_rproblem%rcollection)

    ! Finally release the collection and the problem structure.
    call collct_done (p_rproblem%rcollection)

    deallocate(p_rproblem)

    ! Stop the timer
    call stat_stopTimer(rtimerTotal)

    ! Print the time for the total computation
    call output_lbrk ()
    call output_line ("Time for initial mesh generation:       "//&
        trim(sys_sdL(rtimerGridGeneration%delapsedReal,10)))
    call output_line ("Time for initial matrix initialisation: "//&
      trim(sys_sdL(rtimerGridGeneration%delapsedReal,10)))
    call output_line ("Total Time for solver:                  "//&
      trim(sys_sdL(rtimerSolver%delapsedReal,10)))
    call output_line ("Total time:                             "//&
        trim(sys_sdL(rtimerTotal%delapsedReal,10)))

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
!! https://github.com/flang-compiler/flang/issues/240
!! Fixed as of 3/19/2018, but see flang-20180319{a,b}.f90
!!
!! $ flang -c flang-bug-20170923c.f90
!! clang-5.0: error: unable to execute command: Segmentation fault (core dumped)
!! clang-5.0: error: Fortran frontend to LLVM command failed due to signal (use -v to see invocation)
!! clang version 5.0.1 (https://github.com/flang-compiler/clang.git b11539abc46cbd19189c5719d1e30539de3a93b9) (https://github.com/llvm-mirror/llvm.git 81029f142231bde8e119becda112a2173f1459c9)
!! Target: x86_64-unknown-linux-gnu
!! Thread model: posix
!! InstalledDir: /opt/flang/5.0/bin
!! clang-5.0: note: diagnostic msg: PLEASE submit a bug report to  and include the crash backtrace, preprocessed source, and associated run script.
!! clang-5.0: note: diagnostic msg: Error generating preprocessed source(s) - no preprocessable inputs.

module map_type
  type :: item
    type(item), pointer :: next => null()
  end type
  type :: map
    type(item), pointer :: first => null()
  end type
  type :: parameter_list
    type(map) :: params = map() ! this default initializaton causes the ICE
  end type
end module

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
      SUBROUTINE MB04TU( N, X, INCX, Y, INCY, C, S )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To perform the Givens transformation, defined by C (cos) and S
C     (sin), and interchange the vectors involved, i.e.
C
C        |X(i)|    | 0   1 |   | C   S |   |X(i)|
C        |    | := |       | x |       | x |    |, i = 1,...N.
C        |Y(i)|    | 1   0 |   |-S   C |   |Y(i)|
C
C     REMARK. This routine is a modification of DROT from BLAS.
C             This routine is called only by the SLICOT routines MB04TX
C             and MB04VX.
C
C     NUMERICAL ASPECTS
C
C     The algorithm is backward stable.
C
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Apr. 1997.
C     Supersedes Release 2.0 routine MB04FU by Th.G.J. Beelen,
C     Philips Glass Eindhoven, Holland.
C
C     REVISIONS
C
C     January 26, 1998.
C
C     KEYWORDS
C
C     Othogonal transformation.
C
C     ******************************************************************
C
C     .. Scalar Arguments ..
      INTEGER           INCX, INCY, N
      DOUBLE PRECISION  C, S
C     .. Array Arguments ..
      DOUBLE PRECISION  X(*), Y(*)
C     .. Local Scalars ..
      DOUBLE PRECISION  DTEMP
      INTEGER           I, IX, IY
C     .. Executable Statements ..
C
      IF ( N.LE.0 ) RETURN
      IF ( ( INCX.NE.1 ) .OR. ( INCY.NE.1 ) ) THEN
C
C        Code for unequal increments or equal increments not equal to 1.
C
         IX = 1
         IY = 1
         IF ( INCX.LT.0 ) IX = (-N+1)*INCX + 1
         IF ( INCY.LT.0 ) IY = (-N+1)*INCY + 1
C
         DO 20 I = 1, N
            DTEMP = C*Y(IY) - S*X(IX)
            Y(IY) = C*X(IX) + S*Y(IY)
            X(IX) = DTEMP
            IX = IX + INCX
            IY = IY + INCY
   20    CONTINUE
C
      ELSE
C
C        Code for both increments equal to 1.
C
         DO 40 I = 1, N
            DTEMP = C*Y(I) - S*X(I)
            Y(I) = C*X(I) + S*Y(I)
            X(I) = DTEMP
   40    CONTINUE
C
      END IF
C
      RETURN
C *** Last line of MB04TU ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
Function fdbel(s)
  Parameter (srt0=2.012)
  If (s<=srt0**2) Then
    fdbel = 0.
  Else
    fdbel = 2500.*exp(-(s-7.93)**2/0.003) + 300.*exp(-(s-7.93)**2/0.1) + 10.
  End If
  Return
End Function fdbel

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
	REAL*4 FUNCTION IPPOWER
     .  (XALPHA,XP0,XP1,XN,XM,R,CASECTRL,ALT,UORF)
	IMPLICIT NONE
C
	REAL*4 XALPHA,XP0,XP1,XN,XM
	REAL*4 R
	INTEGER*4 CASECTRL, ALT, UORF
C
C	POWER CALCULATIONS FOR INDEPENDENT PROPORTIONS
C	NORMAL APPROXIMATION ASSUMED
C
C   Power formula from:  Schlesselman:  Case-control Studies: Design,
C   Conduct, Analysis.  New York: Oxford U. Press 1982: 144-152.
C	See also:  FRIEDMAN, FURBERG AND DEMETS:
C	FUNDAMENTALS OF CLINICAL TRIALS. BOSTON: WRIGHT PSG. 1982:75.
C
C INPUT:
C	ALPHA	TYPE I ERROR probability (TWO TAILED)
C	P0	Probability of exposure in control group when CASECTL=Y
C		Probability of event in control group when CASECTL=N
C	P1	Probability of exposure in case group when CASECTL=Y
C		Probability of event in case group when CASECTL=N
C	N	NUMBER OF CASE PATIENTS
C	M	Ratio of control to case
C	R	Odds ratio or relative risk.
C	CASECTL 1 for case-control studies
C		2 for prospective studies
C	ALT	=1 if input is two proportions
C		=2 if input is other (odds ration or relative risk)
C	UORF =1 for uncorrected chi-square test
C        =2 for Fisher's exact test
C OUTPUT:
C	POWER	The probability of correctly rejecting the null
C		hypothesis of equal event (exposure) probabilities
C		with type I error probability ALPHA given P0, P1, N,
C		and M defined above.
C
C  Globals
C
	REAL ALPHA,P1,P0,N,M
	COMMON/IPFISH/ALPHA,P1,P0,N,M
C
C    Declare functions.
C
	EXTERNAL FISHSIZB
	REAL FISHSIZB
C
C  Locals
C
	REAL ZALPHA,ZCRVALUE,PBAR
	REAL Q0,Q1,QBAR,S1,S2
	DOUBLE PRECISION PHI
	INTEGER ERR
	REAL BETA,EPS,E1,Y1,Y2
C
C Code:
C
	ALPHA =	XALPHA
	P0 = XP0
	P1 = XP1
	N = XN
	M = XM

	IF (ALT.EQ.2) THEN
	    IF (CASECTRL.EQ.1) THEN
	    	P1=P0*R/(1.+P0*(R-1.))
	    ELSE
	    	P1=P0*R
	    END IF
	END IF
C
C    Calculate a single value of POWER.
C
	ZALPHA=ZCRVALUE(ALPHA/2.)
	PBAR=(P1+M*P0)/(M+1)
	QBAR=1.-PBAR
	Q1=1.-P1
	Q0=1.-P0
	S1=SQRT(PBAR*QBAR*(1.+1./M)/N)
	S2=SQRT((P0*Q0/M+P1*Q1)/N)
	BETA=REAL(PHI((ZALPHA*S1-(P0-P1))/S2))-
     .		 REAL(PHI((-ZALPHA*S1-(P0-P1))/S2))
C
C	If the test type is Fisher's exact then use the
C	uncorrected beta as a starting point and find
C	the corrected beta.
C
	IF (UORF.EQ.2) THEN
C
C    Set error tolerance values for call to bisection routine.
C
		EPS=0.00002*MIN(P0,1.-P0)
		E1 =0.00002*MIN(P0,1.-P0)
C
C	Set initial end points.  Start at beta.
C
		Y1=BETA
C
C	Power (=1-beta) must not be less than alpha.  Therefore, the
C	largest beta must be 1-alpha.
C
		Y2=1.-ALPHA
C
C	Make the approximation.  Use bisection to find the
C	zero of the function FISHSIZB.
C
		CALL BISEC(Y1,Y2,EPS,E1,FISHSIZB,BETA,ERR)
C
C	Continuity correction is invalid when the power is
C	very low.  Set BETA to 1-ALPHA.
C
		IF (ERR.NE.0.AND.FISHSIZB(Y1).GT.0.AND.FISHSIZB(Y2).GT.0) THEN
C			BETA=Y2
			IPPOWER=-1
			RETURN
	    END IF
	END IF
C
C	Set return value.  Power is 1-beta.
C
	IPPOWER=1.-BETA
C
C	Done.
C
	RETURN
	END

ext: for

max_stars_repo_name: vubiostat/ps

max_issues_repo_name: vubiostat/ps

max_forks_repo_name: vubiostat/ps


text:
      SUBROUTINE KMNS(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
     *    ITRAN, LIVE, ITER, WSS, IFAULT)
C
C     ALGORITHM AS 136  APPL. STATIST. (1979) VOL.28, NO.1
C
C     Divide M points in N-dimensional space into K clusters so that
C     the within cluster sum of squares is minimized.
C
      INTEGER IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K), LIVE(K)
      REAL    A(M,N), D(M), C(K,N), AN1(K), AN2(K), WSS(K), DT(2)
      REAL    ZERO, ONE
C
C     Define BIG to be a very large positive number
C
      DATA BIG /1.E30/, ZERO /0.0/, ONE /1.0/
C
      IFAULT = 3
      IF (K .LE. 1 .OR. K .GE. M) RETURN
C
C     For each point I, find its two closest centres, IC1(I) and
C     IC2(I).     Assign it to IC1(I).
C
      DO 50 I = 1, M
	IC1(I) = 1
	IC2(I) = 2
	DO 10 IL = 1, 2
	  DT(IL) = ZERO
	  DO 10 J = 1, N
	    DA = A(I,J) - C(IL,J)
	    DT(IL) = DT(IL) + DA*DA
   10   CONTINUE
	IF (DT(1) .GT. DT(2)) THEN
	  IC1(I) = 2
	  IC2(I) = 1
	  TEMP = DT(1)
	  DT(1) = DT(2)
	  DT(2) = TEMP
	END IF
	DO 50 L = 3, K
	  DB = ZERO
	  DO 30 J = 1, N
	    DC = A(I,J) - C(L,J)
	    DB = DB + DC*DC
	    IF (DB .GE. DT(2)) GO TO 50
   30     CONTINUE
	  IF (DB .LT. DT(1)) GO TO 40
	  DT(2) = DB
	  IC2(I) = L
	  GO TO 50
   40     DT(2) = DT(1)
	  IC2(I) = IC1(I)
	  DT(1) = DB
	  IC1(I) = L
   50 CONTINUE
C
C     Update cluster centres to be the average of points contained
C     within them.
C
      DO 70 L = 1, K
	NC(L) = 0
	DO 60 J = 1, N
   60   C(L,J) = ZERO
   70 CONTINUE
      DO 90 I = 1, M
	L = IC1(I)
	NC(L) = NC(L) + 1
	DO 80 J = 1, N
   80   C(L,J) = C(L,J) + A(I,J)
   90 CONTINUE
C
C     Check to see if there is any empty cluster at this stage
C
      DO 120 L = 1, K
	IF (NC(L) .EQ. 0) THEN
	  IFAULT = 1
	  RETURN
	END IF
	AA = NC(L)
	DO 110 J = 1, N
  110   C(L,J) = C(L,J) / AA
C
C     Initialize AN1, AN2, ITRAN & NCP
C     AN1(L) = NC(L) / (NC(L) - 1)
C     AN2(L) = NC(L) / (NC(L) + 1)
C     ITRAN(L) = 1 if cluster L is updated in the quick-transfer stage,
C              = 0 otherwise
C     In the optimal-transfer stage, NCP(L) stores the step at which
C     cluster L is last updated.
C     In the quick-transfer stage, NCP(L) stores the step at which
C     cluster L is last updated plus M.
C
	AN2(L) = AA / (AA + ONE)
	AN1(L) = BIG
	IF (AA .GT. ONE) AN1(L) = AA / (AA - ONE)
	ITRAN(L) = 1
	NCP(L) = -1
  120 CONTINUE
      INDX = 0
      DO 140 IJ = 1, ITER
C
C     In this stage, there is only one pass through the data.   Each
C     point is re-allocated, if necessary, to the cluster that will
C     induce the maximum reduction in within-cluster sum of squares.
C
	CALL OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
     *        ITRAN, LIVE, INDX)
C
C     Stop if no transfer took place in the last M optimal transfer
C     steps.
C
	IF (INDX .EQ. M) GO TO 150
C
C     Each point is tested in turn to see if it should be re-allocated
C     to the cluster to which it is most likely to be transferred,
C     IC2(I), from its present cluster, IC1(I).   Loop through the
C     data until no further change is to take place.
C
	CALL QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
     *       ITRAN, INDX)
C
C     If there are only two clusters, there is no need to re-enter the
C     optimal transfer stage.
C
	IF (K .EQ. 2) GO TO 150
C
C     NCP has to be set to 0 before entering OPTRA.
C
	DO 130 L = 1, K
  130   NCP(L) = 0
  140 CONTINUE
C
C     Since the specified number of iterations has been exceeded, set
C     IFAULT = 2.   This may indicate unforeseen looping.
C
      IFAULT = 2
C
C     Compute within-cluster sum of squares for each cluster.
C
  150 DO 160 L = 1, K
	WSS(L) = ZERO
	DO 160 J = 1, N
	  C(L,J) = ZERO
  160 CONTINUE
      DO 170 I = 1, M
	II = IC1(I)
	DO 170 J = 1, N
	  C(II,J) = C(II,J) + A(I,J)
  170 CONTINUE
      DO 190 J = 1, N
	DO 180 L = 1, K
  180   C(L,J) = C(L,J) / FLOAT(NC(L))
	DO 190 I = 1, M
	  II = IC1(I)
	  DA = A(I,J) - C(II,J)
	  WSS(II) = WSS(II) + DA*DA
  190 CONTINUE
C
      RETURN
      END
C
C
      SUBROUTINE OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
     *      ITRAN, LIVE, INDX)
C
C     ALGORITHM AS 136.1  APPL. STATIST. (1979) VOL.28, NO.1
C
C     This is the optimal transfer stage.
C
C     Each point is re-allocated, if necessary, to the cluster that
C     will induce a maximum reduction in the within-cluster sum of
C     squares.
C
      INTEGER IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K), LIVE(K)
      REAL    A(M,N), D(M), C(K,N), AN1(K), AN2(K), ZERO, ONE
C
C     Define BIG to be a very large positive number.
C
      DATA BIG /1.0E30/, ZERO /0.0/, ONE/1.0/
C
C     If cluster L is updated in the last quick-transfer stage, it
C     belongs to the live set throughout this stage.   Otherwise, at
C     each step, it is not in the live set if it has not been updated
C     in the last M optimal transfer steps.
C
      DO 10 L = 1, K
	IF (ITRAN(L) .EQ. 1) LIVE(L) = M + 1
   10 CONTINUE
      DO 100 I = 1, M
	INDX = INDX + 1
	L1 = IC1(I)
	L2 = IC2(I)
	LL = L2
C
C     If point I is the only member of cluster L1, no transfer.
C
	IF (NC(L1) .EQ. 1) GO TO 90
C
C     If L1 has not yet been updated in this stage, no need to
C     re-compute D(I).
C
	IF (NCP(L1) .EQ. 0) GO TO 30
	DE = ZERO
	DO 20 J = 1, N
	  DF = A(I,J) - C(L1,J)
	  DE = DE + DF*DF
   20   CONTINUE
	D(I) = DE * AN1(L1)
C
C     Find the cluster with minimum R2.
C
   30   DA = ZERO
	DO 40 J = 1, N
	  DB = A(I,J) - C(L2,J)
	  DA = DA + DB*DB
   40   CONTINUE
	R2 = DA * AN2(L2)
	DO 60 L = 1, K
C
C     If I >= LIVE(L1), then L1 is not in the live set.   If this is
C     true, we only need to consider clusters that are in the live set
C     for possible transfer of point I.   Otherwise, we need to consider
C     all possible clusters.
C
	  IF (I .GE. LIVE(L1) .AND. I .GE. LIVE(L) .OR. L .EQ. L1 .OR.
     *        L .EQ. LL) GO TO 60
	  RR = R2 / AN2(L)
	  DC = ZERO
	  DO 50 J = 1, N
	    DD = A(I,J) - C(L,J)
	    DC = DC + DD*DD
	    IF (DC .GE. RR) GO TO 60
   50     CONTINUE
	  R2 = DC * AN2(L)
	  L2 = L
   60     CONTINUE
	  IF (R2 .LT. D(I)) GO TO 70
C
C     If no transfer is necessary, L2 is the new IC2(I).
C
	  IC2(I) = L2
	  GO TO 90
C
C     Update cluster centres, LIVE, NCP, AN1 & AN2 for clusters L1 and
C     L2, and update IC1(I) & IC2(I).
C
   70     INDX = 0
	  LIVE(L1) = M + I
	  LIVE(L2) = M + I
	  NCP(L1) = I
	  NCP(L2) = I
	  AL1 = NC(L1)
	  ALW = AL1 - ONE
	  AL2 = NC(L2)
	  ALT = AL2 + ONE
	  DO 80 J = 1, N
	    C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
	    C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
   80     CONTINUE
	  NC(L1) = NC(L1) - 1
	  NC(L2) = NC(L2) + 1
	  AN2(L1) = ALW / AL1
	  AN1(L1) = BIG
	  IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
	  AN1(L2) = ALT / AL2
	  AN2(L2) = ALT / (ALT + ONE)
	  IC1(I) = L2
	  IC2(I) = L1
   90   CONTINUE
	IF (INDX .EQ. M) RETURN
  100 CONTINUE
      DO 110 L = 1, K
C
C     ITRAN(L) = 0 before entering QTRAN.   Also, LIVE(L) has to be
C     decreased by M before re-entering OPTRA.
C
	ITRAN(L) = 0
	LIVE(L) = LIVE(L) - M
  110 CONTINUE
C
      RETURN
      END
C
C
      SUBROUTINE QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
     *    ITRAN, INDX)
C
C     ALGORITHM AS 136.2  APPL. STATIST. (1979) VOL.28, NO.1
C
C     This is the quick transfer stage.
C     IC1(I) is the cluster which point I belongs to.
C     IC2(I) is the cluster which point I is most likely to be
C         transferred to.
C     For each point I, IC1(I) & IC2(I) are switched, if necessary, to
C     reduce within-cluster sum of squares.  The cluster centres are
C     updated after each step.
C
      INTEGER IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K)
      REAL    A(M,N), D(M), C(K,N), AN1(K), AN2(K), ZERO, ONE
C
C     Define BIG to be a very large positive number
C
      DATA BIG /1.0E30/, ZERO /0.0/, ONE /1.0/
C
C     In the optimal transfer stage, NCP(L) indicates the step at which
C     cluster L is last updated.   In the quick transfer stage, NCP(L)
C     is equal to the step at which cluster L is last updated plus M.
C
      ICOUN = 0
      ISTEP = 0
   10 DO 70 I = 1, M
	ICOUN = ICOUN + 1
	ISTEP = ISTEP + 1
	L1 = IC1(I)
	L2 = IC2(I)
C
C     If point I is the only member of cluster L1, no transfer.
C
	IF (NC(L1) .EQ. 1) GO TO 60
C
C     If ISTEP > NCP(L1), no need to re-compute distance from point I to
C     cluster L1.   Note that if cluster L1 is last updated exactly M
C     steps ago, we still need to compute the distance from point I to
C     cluster L1.
C
	IF (ISTEP .GT. NCP(L1)) GO TO 30
	DA = ZERO
	DO 20 J = 1, N
	  DB = A(I,J) - C(L1,J)
	  DA = DA + DB*DB
   20   CONTINUE
	D(I) = DA * AN1(L1)
C
C     If ISTEP >= both NCP(L1) & NCP(L2) there will be no transfer of
C     point I at this step.
C
   30   IF (ISTEP .GE. NCP(L1) .AND. ISTEP .GE. NCP(L2)) GO TO 60
	R2 = D(I) / AN2(L2)
	DD = ZERO
	DO 40 J = 1, N
	  DE = A(I,J) - C(L2,J)
	  DD = DD + DE*DE
	  IF (DD .GE. R2) GO TO 60
   40   CONTINUE
C
C     Update cluster centres, NCP, NC, ITRAN, AN1 & AN2 for clusters
C     L1 & L2.   Also update IC1(I) & IC2(I).   Note that if any
C     updating occurs in this stage, INDX is set back to 0.
C
	ICOUN = 0
	INDX = 0
	ITRAN(L1) = 1
	ITRAN(L2) = 1
	NCP(L1) = ISTEP + M
	NCP(L2) = ISTEP + M
	AL1 = NC(L1)
	ALW = AL1 - ONE
	AL2 = NC(L2)
	ALT = AL2 + ONE
	DO 50 J = 1, N
	  C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
	  C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
   50   CONTINUE
	NC(L1) = NC(L1) - 1
	NC(L2) = NC(L2) + 1
	AN2(L1) = ALW / AL1
	AN1(L1) = BIG
	IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
	AN1(L2) = ALT / AL2
	AN2(L2) = ALT / (ALT + ONE)
	IC1(I) = L2
	IC2(I) = L1
C
C     If no re-allocation took place in the last M steps, return.
C
   60   IF (ICOUN .EQ. M) RETURN
   70 CONTINUE
      GO TO 10
      END

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
      double precision function determ(array,nord,nrows)
c
c     calculate determinate of a square matrix
c        (from bevington "data reduction and error analysis
c         for the physical sciences" pg 294)
c     array: matrix to be analyzed
c     nord: order of matrix
c     nrows:  first dimension of matrix in calling routine
c
      double precision array(nrows,nrows)
      determ = 1.
      do 150 k=1,nord
c
c
        if (array(k,k).ne.0) go to 130
        do 100 j=k,nord
          if (array(k,j).ne.0) go to 110
  100   continue
        determ = 0.
        go to 160
c
  110   do 120 i=k,nord
          saved = array(i,j)
          array(i,j) = array(i,k)
  120   array(i,k) = saved
        determ = -determ
c
  130   determ = determ*array(k,k)
        if (k.ge.nord) go to 150
        k1 = k+1
        do 140 i=k1,nord
          do 140 j=k1,nord
  140   array(i,j) = array(i,j)-array(i,k)*array(k,j)/array(k,k)
  150 continue
  160 return
c end double precision function determ
      end


text:
Im a UC Davis student interested in environmental sustainability, ethical consumption, community development and many other related things. I havent used my real name at this point in time but I plan to create an account with my real name soon.

please consider the importance of using your realname Daubert

20110930 07:31:58 nbsp Welcome to the wiki! Users/CovertProfessor

20111003 08:35:36 nbsp Glad to hear youre enjoying Davis so far!  I was interested to hear about the ways in which West Village was falling short on their promises.  Hopefully, theyll start to improve on those areas soon, especially in the light of feedback from people like you! Users/CovertProfessor

20111005 05:23:15 nbsp To create a new entry, just link to it or search for the title and when you follow the link, or on the search result page, simply click create a new page.  For instance, you could follow the link Occupy Davis and create an entry under that title.  Dont forget to link to it from appropriate entries so people can find it!  Users/JabberWokky Evan JabberWokky Edwards

20111005 08:33:06 nbsp Glad to hear the movement is spreading.  JabberWokky pretty much explained how to create the new page  just click on the link for Occupy Davis that he made, and go from there!  If you need help with editing the page, I or others will be glad to help out. Users/CovertProfessor

20111005 19:03:14 nbsp The FB page for Occupy Davis says Friday, October 7 at 6:00pm, Davis Farmers Market.  Should the wiki go with that??  There could always be a march towards the bank area.
 Users/CovertProfessor

20111005 19:18:23 nbsp I was also thinking that would be a good plan. I think the thing on friday night is just a general assembly to get organised with supplies and discuss agenda. I guess suggestions will come out of that. Feel free to edit it with whatever is most current/relevant. Thanks for your help! Users/ConsciousConsumer

20111005 19:23:34 nbsp Sure, thanks for getting the ball rolling!! Users/CovertProfessor

20111005 19:27:42 nbsp I just tried to create the featured page but I suck at it, still figuring stuff out. Users/ConsciousConsumer

20111005 19:28:42 nbsp Heh... no worries, youll get it. Ill take a look and see if I can help. Users/CovertProfessor

20111006 21:54:55 nbsp I added Occupy Davis to the Activism page! Users/CovertProfessor

20111006 23:02:16 nbsp Nah, just the cleanup crew.  :) Users/CovertProfessor

20111007 13:20:17 nbsp You should use your real name if youre heading up Occupy Davis. It would up your credibility tremendously. Users/ChristyMarsden

20111015 17:15:14 nbsp Thanks for posting  that Rolling Stones link Users/SueHjerpe

20111015 22:27:31 nbsp Those costs will get passed along to the consumer one way or another.  Either through higher food costs, forced gratuity or lower food quality. Thinking anything else is fantasy.   Users/OliviaY
  In most of Western Europe, you dont pay a tip.  The employees are paid a decent wage for working for restaurants.  The quality of the food is about the same or better than here.  In some countries, such as Italy, you pay a different fee instead of a tip.  That fee goes to the restaurant, not to the workers.  (The prices for food in Italy are MUCH lower than they are here, except in the tourist areas, where they are a little bit more expensive at current exchange rates.)  Where you do see some interesting differences are in charges for things like water (free in the US, but you must buy a bottle in most of Europe, which is environmentally horrible), and in the customer service, which ranges from not so great to truly terrible  but this varies significantly by country.   Users/IDoNotExist
    Considering how bad things have been in western europe I wouldnt exactly compare that hellhole to the rapidly developing hell here.  Unless you want faster free fall into total anarchy which hey, lets be bad guys!
  You get some attitude if youre a tourist but I would say that the attitude is hardly any better here by comparison. Water isnt free here, and I would argue that most people buy bottled water here or have it filtered (particularly in Davis with the quality of the water).
A study from a Cornell University professor concluded that consumers assessments of the quality of service correlate weakly to the amount they tip, relying more on superficial factors such as the attractiveness of the server or the amount of the bill.

20111016 19:38:06 nbsp I have been outside the country but its irrelevant.  Businesses will pass the cost on to the consumer so the servers might get a living wage, but youll be footing the bill.  Your comment made it sound like you expect the business to profit less for the benefit of servers and THAT will never happen. Users/OliviaY
  Its not irrelevant if youve seen it work, and that was my point. Why are you so pessimistic? What are your suggestions? Leave things as they are??

20111016 19:44:56 nbsp and who is going to enforce that? No thanks.  More government is never the answer and forcing people to do good by holding a gun to their heads is counterproductive. Users/OliviaY
 It doesnt require more government, just more effective government. Get rid of the oversized military.
  If all we do is raise minimum wage/enforce payment of a living wage that will result only in a pass along of the cost to consumer through one of the three ways I listed.  IF you want to keep food quality up and not see consumers pay for the cost increase, that would require more governing through fixed costs on how each business operates.  THAT is the very definition of more government.  (Military is a total tangent but I agree with you there.)
   Food quality is shit in this country anyway. High fuctose corn syrup in everything is evidence of that. The government needs to stop subsidising big business through paying medicaid to walmart workers (they intentionally pay their workers low enough so that they qualify), and sending U.S. factory jobs overseas under the guise of the USAID program. I want to hear some of your solutions? Because all youve done so far is criticize.

20111016 21:27:49 nbsp Im a libertarian.  My views on how to fix things are shouted down so I dont bother.  Food quality issues start and end with government.  Corn subsidies and the subsequent need for oil because of our government supported monoculture has jacked shit up across the board.   Users/OliviaY

20111016 21:37:38 nbsp Dont be afraid of voicing your opinion. They wont get shouted down if you are respectful. I actually prefer Libertarians to Republicans. Ron Paul makes a whole lot of sense. I just dont see how the poorest get taken care of. I really want to hear some of your ideas for a better United States. I dont align myself with any parties. Have you done the political compass? Users/ConsciousConsumer

20111018 20:12:35 nbsp I think federal minimum wage is a bad idea.  At the very least it should be a state by state issue.  Either way that should be the salary politicians make.  Im a big fan of that idea. Users/OliviaY

20111024 19:52:56 nbsp Ok, I made some small fixes and did some other cleanup on the page.  Glad to be able to help in this small way.  If there are things the protesters need  ways to donate money or materials  that would be good to include on the page, too. Users/CovertProfessor

20111026 15:36:08 nbsp I was wondering if you would be interested in doing an interview for an article I am writing on UC Davis West Village. While I applaud the project I am interested in getting an interview from an especially environmentally conscious student who can talk to some of the short comings of the project. We can conduct the interview online or in person. Let me know if you are interested.  Users/CatherineHawe

20111201 15:37:24 nbsp Thanks for posting on Occupy UC Davis. I like your comment about documentary showings. I dont have a laptop but have access to a projector.  Im an alumni helping out.  I think it might be better if you, as a student, could get access to a campus room with projector (much more likely to attract people to an indoor venue). If you know about copyright issues...we may have to ask permission to show some videos. Inside Job is an important one to share. I already showed Money Masters.  best, Derek....email at hiscomposty@yahoo.com
 Users/DerekDowney

20111208 06:45:21 nbsp This is a quote from the Collected Works of Mahatma Ghandi “Why, of all places in Johannesburg, the Indian Location should be chosen for dumping down all the Kaffirs of the town passes my comprehension…the Town Council must withdraw the Kaffirs from the Location.” (Reference: CWMG, Vol I, pp. 244245) Ghandi found success in South Africa for Indians by stepping on Africans, so I dont think Ghandi is so great as people say. Users/NickSchmalenberger


text:
      subroutine dtimer(ttime)
      double precision ttime
c     *********
c
c     Subroutine dtimer
c
c     This subroutine is used to determine user time. In a typical
c     application, the user time for a code segment requires calls
c     to subroutine dtimer to determine the initial and final time.
c
c     The subroutine statement is
c
c       subroutine dtimer(ttime)
c
c     where
c
c       ttime is a double precision variable.
c         On input ttime need not be specified.
c         On output ttime specifies the user time.
c
c     MINPACK-2 Project. November 1993.
c     Argonne National Laboratory and University of Minnesota.
c     Brett M. Averick and Jorge J. More'.
c
c     **********

c     For machines with etime (sun).

      real tarray(2)
      real etime

      ttime = etime(tarray)
      ttime = dble(tarray(1))

c     For machines with mclock (ibm).

c      integer mclock

c      ttime = dble(mclock())/100.0d0

c     For machines with second (cray).

c     real second

c     ttime = second()

      end

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
module game

use, intrinsic :: iso_fortran_env, only : stdin=>input_unit, stdout=>output_unit
use, intrinsic :: iso_c_binding, only: c_int

implicit none (external)
private
public :: debug, mix, hand

logical :: debug

contains

integer(c_int) function hand(cards) result(win)

integer(c_int), intent(inout) :: cards(52)
integer(c_int) :: P,D,pace,dace,cnt,i,ios
character :: yesno

i=1
P=0
D=0
PACE=0
DACE=0
WIN=0

call hit(P,PACE,i,cards)
call hit(D,DACE,i,cards)
call hit(P,PACE,i,cards)
call hit(D,DACE,i,cards)
cnt = 0

print '(A,I2)', 'Dealer shows: ',cards(i-1)
print '(A,I2,A,I1)', 'Player total: ',P,' incl. # of Aces: ',pace

if (p == 21) then
   print *, 'Player Blackjack!'
   win = 1
else
   cnt = 2
   do
     write(stdout,'(A)', advance='no') 'Hit?  y / n  '
     read(stdin, '(A1)',iostat=ios) yesno
     if(is_iostat_end(ios)) stop 'Goodbye'
     if (yesno /= 'y') exit

     call hit(p,pace,i,cards)
     cnt = cnt+1
     print '(A,I2,A,I1)', 'Player ',P,' # of Aces: ',pace
     if (p > 21) then
       print *, 'player busts -- dealer wins'
       return
     endif
   enddo
endif

if (d == 21) then
  print *, 'Dealer Blackjack'
  if (win == 1) then
    print *,' ---- Push'
    win = 2
    return
  else
    print *, 'Dealer automatically wins'
    return
  endif
else
  if ((p == 21) .or. (cnt >= 5)) then
    print *, 'Player automatically wins'
    win = 1
    return
  endif
endif

print '(A,I2)', 'Dealer has ', D
do while (D <= 16)
  call hit(d,dace,i,cards)
  print '(A,I2)', 'Dealer has ', D
  if (d > 21) then
    print *,'Dealer busts - Player wins'
    win = 1
    return
  endif
enddo

print '(A,I2,A,I2)', 'Player: ',P,' Dealer: ',D
if (p > d) then
  print *,'Player wins'
  win = 1
else
  print *, 'Dealer wins'
endif

end function hand

subroutine hit(total, aces, i, cards) bind(c)

integer(c_int), intent(inout) :: total, aces, i
integer(c_int), intent(inout) :: cards(52)

if (debug) then
  write(stdout,'(A)', advance='no') 'input next card: '
  read(stdin, '(I2)') cards(i)
endif

total = total + cards(i)
if (cards(i) == 11) then
  aces = aces + 1
endif

i = i + 1

if (total > 21 .and. aces >= 1) then
  total = total - 10
  aces = aces - 1
endif

end subroutine hit

subroutine mix(cards) bind(c)
! knuth shuffle
! https://www.rosettacode.org/wiki/Knuth_shuffle#Fortran

integer(c_int), intent(out) :: cards(52)
integer(c_int) :: i, j, temp
real :: r

cards = &
[2,2,2,2, &
 3,3,3,3, &
 4,4,4,4, &
 5,5,5,5, &
 6,6,6,6, &
 7,7,7,7, &
 8,8,8,8, &
 9,9,9,9, &
 10,10,10,10, 10,10,10,10, 10,10,10,10, 10,10,10,10, &
 11,11,11,11]

do i = size(cards), 2, -1
  call random_number(r)
  j = int(r * i) + 1
  temp = cards(j)
  cards(j) = cards(i)
  cards(i) = temp
end do

end subroutine mix

end module game

max_stars_repo_name: fortran-gaming/NIST-Blackjack

max_issues_repo_name: fortran-gaming/NIST-Blackjack

max_forks_repo_name: fortran-gaming/NIST-Blackjack