text:
**==cthei.spg  processed by SPAG 4.50J  at 14:49 on 30 Jun 1995
      FUNCTION CTHEI(N,J,T,Dene)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL CTHEI , d1 , Dene , p3 , qdw , qup , r21 , T , t4
      INTEGER J , N
C*** End of declarations inserted by SPAG
C  SCOTT'S RATES FOR HE+  + J  TO HE0  + J+1
C  THROUGH ARGON  0.1 < T4 < 10.
      CTHEI = 0.
      t4 = 0.0001*T
      IF ( J.EQ.1 .OR. J.GE.4 ) RETURN
      IF ( J.NE.3 ) THEN
         IF ( N.EQ.6 ) CTHEI = 5.3E-10*EXP(-13.2/t4)
C  ALBERT'S RATES FOR N+ FROM 1D LEVEL
         IF ( N.EQ.7 ) THEN
            qdw = 5.14E-06*Dene/SQRT(T)
            qup = 2.86E-6*EXP(-21826./T)
            r21 = qup/(qdw+.0041)
            d1 = r21/(1.+r21)
            p3 = 1. - d1
            CTHEI = p3*4.1E-10*EXP(-10.4/t4) + d1*1.1E-10*EXP(-3.59/t4)
         ENDIF
         IF ( N.EQ.14 ) CTHEI = 2.9E-10*EXP(-8.7/t4)
         IF ( N.EQ.18 ) CTHEI = 1.1E-10*EXP(-3.57/t4)
         RETURN
      ENDIF
      IF ( N.EQ.14 ) CTHEI = 3.4E-9*EXP(-10.5/t4)
      IF ( N.EQ.16 ) CTHEI = 4.8E-10*EXP(-15.7/t4)
      IF ( N.EQ.26 ) CTHEI = 1.2E-9*SQRT(t4)*EXP(-12./t4)
      RETURN
      END

ext: f

max_stars_repo_name: DougBurke/xspeclmodels

max_issues_repo_name: DougBurke/xspeclmodels

max_forks_repo_name: DougBurke/xspeclmodels


text:
      subroutine vadd(n, x, w, y, z)
      integer n, incx, incy
      double precision one
      double precision yy(n)
      double precision w(*), x(*), y(*), z(*)
c
c   vadd
c   in
c   w : vector, y : vector, z : vector
c   out
c   x : vector
c   {
c    x = w + y + z
c   }
c
      incx = 1
      incy = 1
      one = 1.0d0
c
c  Copy y into yy so that input is not overwritten.
c
      call dcopy(n, y, incx, yy, incy)
c
c  y = w + y
c
      call daxpy(n, one, w, incx, yy, incy)
c
c  y = z + y
c
      call daxpy(n, one, z, incx, yy, incy)
c
c  Copy result into x.
c
      call dcopy(n, yy, incx, x, incy)
c
      return
      end

max_stars_repo_name: parsabee/Orio

max_issues_repo_name: parsabee/Orio

max_forks_repo_name: parsabee/Orio


MODULE phycst
  USE par_oce
  USE in_out_manager
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: phy_cst
  REAL(KIND = wp), PUBLIC :: rpi = 3.141592653589793_wp
  REAL(KIND = wp), PUBLIC :: rad = 3.141592653589793_wp / 180._wp
  REAL(KIND = wp), PUBLIC :: rsmall = 0.5 * EPSILON(1.E0)
  REAL(KIND = wp), PUBLIC :: rday = 24. * 60. * 60.
  REAL(KIND = wp), PUBLIC :: rsiyea
  REAL(KIND = wp), PUBLIC :: rsiday
  REAL(KIND = wp), PUBLIC :: raamo = 12._wp
  REAL(KIND = wp), PUBLIC :: rjjhh = 24._wp
  REAL(KIND = wp), PUBLIC :: rhhmm = 60._wp
  REAL(KIND = wp), PUBLIC :: rmmss = 60._wp
  REAL(KIND = wp), PUBLIC :: omega
  REAL(KIND = wp), PUBLIC :: ra = 6371229._wp
  REAL(KIND = wp), PUBLIC :: grav = 9.80665_wp
  REAL(KIND = wp), PUBLIC :: rt0 = 273.15_wp
  REAL(KIND = wp), PUBLIC :: rau0
  REAL(KIND = wp), PUBLIC :: r1_rau0
  REAL(KIND = wp), PUBLIC :: rcp
  REAL(KIND = wp), PUBLIC :: r1_rcp
  REAL(KIND = wp), PUBLIC :: rau0_rcp
  REAL(KIND = wp), PUBLIC :: r1_rau0_rcp
  REAL(KIND = wp), PUBLIC :: emic = 0.97_wp
  REAL(KIND = wp), PUBLIC :: sice = 6.0_wp
  REAL(KIND = wp), PUBLIC :: soce = 34.7_wp
  REAL(KIND = wp), PUBLIC :: rLevap = 2.5E+6_wp
  REAL(KIND = wp), PUBLIC :: vkarmn = 0.4_wp
  REAL(KIND = wp), PUBLIC :: stefan = 5.67E-8_wp
  REAL(KIND = wp), PUBLIC :: rhos = 330._wp
  REAL(KIND = wp), PUBLIC :: rhoi = 917._wp
  REAL(KIND = wp), PUBLIC :: rhow = 1000._wp
  REAL(KIND = wp), PUBLIC :: rcnd_i = 2.034396_wp
  REAL(KIND = wp), PUBLIC :: rcpi = 2067.0_wp
  REAL(KIND = wp), PUBLIC :: rLsub = 2.834E+6_wp
  REAL(KIND = wp), PUBLIC :: rLfus = 0.334E+6_wp
  REAL(KIND = wp), PUBLIC :: rTmlt = 0.054_wp
  REAL(KIND = wp), PUBLIC :: r1_rhoi
  REAL(KIND = wp), PUBLIC :: r1_rhos
  REAL(KIND = wp), PUBLIC :: r1_rcpi
  CONTAINS
  SUBROUTINE phy_cst
    rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp
    rsiday = rday / (1._wp + rday / rsiyea)
    omega = 2._wp * rpi / rsiday
    r1_rhoi = 1._wp / rhoi
    r1_rhos = 1._wp / rhos
    r1_rcpi = 1._wp / rcpi
    IF (lwp) THEN
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) 'phy_cst : initialization of ocean parameters and constants'
      WRITE(numout, FMT = *) '~~~~~~~'
      WRITE(numout, FMT = *) '      mathematical constant                 rpi = ', rpi
      WRITE(numout, FMT = *) '      day                                rday   = ', rday, ' s'
      WRITE(numout, FMT = *) '      sideral year                       rsiyea = ', rsiyea, ' s'
      WRITE(numout, FMT = *) '      sideral day                        rsiday = ', rsiday, ' s'
      WRITE(numout, FMT = *) '      omega                              omega  = ', omega, ' s^-1'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      nb of months per year               raamo = ', raamo, ' months'
      WRITE(numout, FMT = *) '      nb of hours per day                 rjjhh = ', rjjhh, ' hours'
      WRITE(numout, FMT = *) '      nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
      WRITE(numout, FMT = *) '      nb of seconds per minute            rmmss = ', rmmss, ' s'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      earth radius                         ra   = ', ra, ' m'
      WRITE(numout, FMT = *) '      gravity                              grav = ', grav, ' m/s^2'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      freezing point of water              rt0  = ', rt0, ' K'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '   reference density and heat capacity now defined in eosbn2.f90'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      thermal conductivity of pure ice          = ', rcnd_i, ' J/s/m/K'
      WRITE(numout, FMT = *) '      thermal conductivity of snow is defined in a namelist '
      WRITE(numout, FMT = *) '      fresh ice specific heat                   = ', rcpi, ' J/kg/K'
      WRITE(numout, FMT = *) '      latent heat of fusion of fresh ice / snow = ', rLfus, ' J/kg'
      WRITE(numout, FMT = *) '      latent heat of subl.  of fresh ice / snow = ', rLsub, ' J/kg'
      WRITE(numout, FMT = *) '      density of sea ice                        = ', rhoi, ' kg/m^3'
      WRITE(numout, FMT = *) '      density of snow                           = ', rhos, ' kg/m^3'
      WRITE(numout, FMT = *) '      density of freshwater (in melt ponds)     = ', rhow, ' kg/m^3'
      WRITE(numout, FMT = *) '      salinity of ice (for pisces)              = ', sice, ' psu'
      WRITE(numout, FMT = *) '      salinity of sea (for pisces and isf)      = ', soce, ' psu'
      WRITE(numout, FMT = *) '      latent heat of evaporation (water)        = ', rLevap, ' J/m^3'
      WRITE(numout, FMT = *) '      von Karman constant                       = ', vkarmn
      WRITE(numout, FMT = *) '      Stefan-Boltzmann constant                 = ', stefan, ' J/s/m^2/K^4'
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      conversion: degre ==> radian          rad = ', rad
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) '      smallest real computer value       rsmall = ', rsmall
    END IF
  END SUBROUTINE phy_cst
END MODULE phycst


text:
program prepro
   use TermClass
   use DictionaryClass
   use PreprocessingClass
   use PostProcessingClass
   implicit none
   type(Dictionary_)       :: InfileList
   type(MPI_)              :: MPIData
   type(PreProcessing_)    :: Root,Soil
   type(Term_)             :: term
   integer :: ans
   character * 200 :: name,ElemType,SolverName
   call MPIData%Start()
   call term%Init()
   ElemType = 'LinearRectangularGp4'

   call InfileList%Init(1)
   call InfileList%Input(1,'Tutorial/LeafImageData/grass_traced.png')
   name = InfileList%Get(MPIData%MyRank+1)
   print *, "My_rank : ",MPIData%MyRank,"InfileName : ",trim(name)
   call Root%Init(Default=.true.)
   call Root%ImportPictureName(name)
   call Root%GetPixcelSize(MPIData)
   call Root%SetColor(0,128,0)
   call Root%ConvertMesh2Scf(MPIData,ElementType=ElemType)
   call Root%FEMDomain%checkconnectivity(fix=.true.)
   call Root%Convert3Dto2D()
   call Root%SetScale(scalex=7.10d0,scaley=113.350d0)
   SolverName="DiffusionEq_"
   call Root%SetSolver(InSolverType=SolverName)
   call Root%SetUp(NoFacetMode=.true.)
   call Root%Reverse()
   call Root%Convert2Dto3D(Thickness=0.250d0,division=4)
   call Root%SetMatID(MaterialID=1,Xmin=-1000000000000.00d0,&
Xmax=1000000000000.00d0,&
Ymin=-1000000000000.00d0,&
Ymax=1000000000000.00d0,&
Zmin=-1000000000000.00d0,&
Zmax=1000000000000.00d0,&
Tmin=-1000000000000.00d0,&
Tmax=1000000000000.00d0)
   call Root%SetSizeOfBC(Dirichlet=.true.,NumOfValue=1)
   call Root%SetBC(Dirichlet=.true., val=10.00d0,val_id=1,Xmin=-1000000000000.00d0,&
Xmax=1000000000000.00d0,&
Ymin=-25.00d0,&
Ymax=-0.00d0,&
Zmin=-1000000000000.00d0,&
Zmax=1000000000000.00d0,&
Tmin=-1000000000000.00d0,&
Tmax=1000000000000.00d0)
   call Root%SetBC(Dirichlet=.true., val=40.00d0,val_id=1,Xmin=-1000000000000.00d0,&
Xmax=1000000000000.00d0,&
Ymin=-1000000000000.00d0,&
Ymax=-100.00d0,&
Zmin=-1000000000000.00d0,&
Zmax=1000000000000.00d0,&
Tmin=-1000000000000.00d0,&
Tmax=1000000000000.00d0)
   call Root%SetControlPara(OptionalItrTol=100,OptionalTimestep=100,OptionalSimMode=1)
   call Root%Export(Name="Tutorial/InputData/test")

   call MPIData%End()
end program prepro

max_stars_repo_name: kazulagi/plantFEM

max_issues_repo_name: kazulagi/plantFEM

max_forks_repo_name: kazulagi/plantFEM


text:
module mod_monolis_precond_ROM
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  use mod_monolis_util
  use mod_monolis_precond_rom_gs

  implicit none

  private

  integer(kint), save :: n_get_eigen
  real(kdouble), save, allocatable :: val(:), vec(:,:)

  public :: monolis_precond_ROM_setup
  public :: monolis_precond_ROM_apply
  public :: monolis_precond_ROM_clear

contains

  subroutine monolis_precond_ROM_setup(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    integer(kint) :: maxiter, NPNDOF, n_bc
    real(kdouble) :: ths
    logical, allocatable :: is_bc(:)

    call monolis_precond_rom_gs_setup(monoPRM, monoCOM, monoMAT)

    allocate(is_bc(monoMAT%NP*monoMAT%NDOF), source = .false.)
    call get_is_bc(monoPRM, monoCOM, monoMAT, is_bc, n_bc)
    !write(*,*)"n_bc", n_bc

    !maxiter = monoMAT%NP*monoMAT%NDOF
    maxiter = min(50, monoMAT%NP*monoMAT%NDOF)
    n_get_eigen = 10 !monoMAT%NP*monoMAT%NDOF - n_bc !min(10, monoMAT%NP*monoMAT%NDOF)
    ths = 1.0d-6
    allocate(val(n_get_eigen), source = 0.0d0)
    allocate(vec(monoMAT%NP*monoMAT%NDOF, n_get_eigen), source = 0.0d0)

    monoPRM%precond = monolis_prec_MUMPS
    monoPRM%show_summary = .false.
    monoPRM%show_time = .false.
    monoPRM%show_iterlog = .false.

    !write(*,*)"n_get_eigen", n_get_eigen
    !write(*,*)"maxiter", maxiter
    !write(*,*)"ths", ths

    NPNDOF = monoMAT%NP*monoMAT%NDOF
!    call interface_monolis_eigen_inverted_standard_lanczos_ &
!    & (monoPRM, monoCOM, monoMAT, NPNDOF, n_get_eigen, ths, maxiter, val, vec, is_bc)

    monoPRM%show_summary = .true.
    monoPRM%show_time = .true.
    monoPRM%show_iterlog = .true.

    monoPRM%precond = monolis_prec_ROM
  end subroutine monolis_precond_ROM_setup

  subroutine get_is_bc(monoPRM, monoCOM, monoMAT, is_bc, n_bc)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    integer(kint) :: i, j, jS, jE, k, l, in, jn, NDOF, NDOF2, n_bc
    integer(kint), allocatable :: flag(:)
    logical :: is_bc(:)

    NDOF  = monoMAT%NDOF
    NDOF2 = NDOF*NDOF
    allocate(flag(NDOF))

    n_bc = 0
    do i = 1, monoMAT%N
      jS = monoMAT%index(i-1) + 1
      jE = monoMAT%index(i  )
      flag = 0
      do j = jS, jE
        in = monoMAT%item(j)
        do k = 1, NDOF
          do l = 1, NDOF
            jn = NDOF2*(j-1)+NDOF*(k-1)+l
            if(monoMAT%A(jn) /= 0.0d0) flag(k) = flag(k) + 1
          enddo
        enddo
      enddo
      do j = 1, NDOF
        if(flag(j) == 1)then
          is_bc(NDOF*(i-1)+j) = .true.
          n_bc = n_bc + 1
        endif
      enddo
    enddo
  end subroutine get_is_bc

  subroutine monolis_precond_ROM_apply(monoPRM, monoCOM, monoMAT, X, Y)
    use mod_monolis_matvec
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    integer(kint), save :: rom_gs_iter = 1
    integer(kint) :: i, j
    real(kdouble) :: X(:), Y(:), phi, t1, t2
    real(kdouble), allocatable :: coef(:), dU(:), U(:), G(:), T(:)

    allocate(coef(n_get_eigen), source = 0.0d0)
    allocate(dU(monoMAT%NP*monoMAT%NDOF), source = 0.0d0)
    allocate(U(monoMAT%NP*monoMAT%NDOF), source = 0.0d0)
    allocate(G(monoMAT%NP*monoMAT%NDOF), source = 0.0d0)
    allocate(T(monoMAT%NP*monoMAT%NDOF), source = 0.0d0)

    Y = X
!    U = X
!    !U = 0.0d0
!    call monolis_precond_rom_gs_apply(monoPRM, monoCOM, monoMAT, X, U, rom_gs_iter)
!    call monolis_residual(monoCOM, monoMAT, U, X, G, t1, t2)
!
!    do i = 1, n_get_eigen
!      phi = 0.0d0
!      do j = 1, monoMAT%NP*monoMAT%NDOF
!        phi = phi + vec(j,i)*G(j)
!      enddo
!      coef(i) = phi/val(i)
!    enddo
!
!    dU = 0.0d0
!    do i = 1, n_get_eigen
!      do j = 1, monoMAT%NP*monoMAT%NDOF
!        dU(j) = dU(j) + vec(j,i)*coef(i)
!      enddo
!    enddo
!
!    U = U + dU
!    call monolis_precond_rom_gs_apply(monoPRM, monoCOM, monoMAT, X, U, rom_gs_iter)
!    Y = U
  end subroutine monolis_precond_ROM_apply

  subroutine monolis_precond_ROM_clear(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
  end subroutine monolis_precond_ROM_clear
end module mod_monolis_precond_ROM

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
      SUBROUTINE SW_R DERIV

      implicit real*8 (a-h,o-z)

c***********************************************************************

c  Change History:
c     26 Apr 96     Changed to generic functions.

c*******************!***************************************************

      complex  * 16 c,s,csq,ssq,f,dfdtht,
     &              hg,norm11,norm22,norm12,rbar11,rbar22

      common/sw_mode/
     &              omega,wn,thetar,thetai,c,s,csq,ssq,f,dfdtht,
     &              hg,norm11,norm22,norm12,rbar11,rbar22,
     &              nriter,adjflg,isotrp

      complex  * 16 m11,m12,m13,m21,m22,m23,m31,m32,m33

      common/sw_mmtx/
     &              m11,m12,m13,m21,m22,m23,m31,m32,m33,ht

      complex  * 16 r11,r22,r12,r21,
     &              logr11,logr22,logr12,logr21,
     &              dl11dh,dl22dh,dl12dh,dl21dh

      common/sw_rmtx/
     &              r11,r22,r12,r21,
     &              logr11,logr22,logr12,logr21,
     &              dl11dh,dl22dh,dl12dh,dl21dh,delh

      complex  * 16 k2i,capd,
     &              t11,t31,t42,t44,t12vrc,t14vrc,t32vrc,t34vrc,ct41,
     &              s11a,d11a,s11b,d11b,c12,c21,
     &              s12,d12,s21,d21,s22,d22,b11,b22,b12,b21

      k2i=DCMPLX(0.d0,-0.5d0*wn)

      if (DREAL(logr11) .gt.  15.d0) then
         logr11= 15.d0
      else
     &if (DREAL(logr11) .lt. -15.d0) then
         logr11=-15.d0
      end if
      r11=EXP(logr11)
      if (DREAL(logr22) .gt.  15.d0) then
         logr22= 15.d0
      else
     &if (DREAL(logr22) .lt. -15.d0) then
         logr22=-15.d0
      end if
      r22=EXP(logr22)
      if (DREAL(logr12) .gt.  15.d0) then
         logr12= 15.d0
      else
     &if (DREAL(logr12) .lt. -15.d0) then
         logr12=-15.d0
      end if
      r12=EXP(logr12)
      if (DREAL(logr21) .gt.  15.d0) then
         logr21= 15.d0
      else
     &if (DREAL(logr21) .lt. -15.d0) then
         logr21=-15.d0
      end if
      r21=EXP(logr21)

      capd=1.d0/(1.d0+m33)

      t11=-s*m31*capd
      t12vrc=s*m32*capd/c
      t14vrc=(csq+m33)*capd/c
      t31=m23*m31*capd-m21
      t32vrc=c+(m22-m23*m32*capd)/c
      t34vrc=s*m23*capd/c
      ct41=(1.d0+m11-m13*m31*capd)*c
      t42=m32*m13*capd-m12
      t44=-s*m13*capd

      s11a=t11+t44
      d11a=t11-t44
      s11b=t14vrc+ct41
      d11b=t14vrc-ct41
      s12=t12vrc+t42
      d12=t12vrc-t42
      s21=t34vrc+t31
      d21=t34vrc-t31
      s22=c+t32vrc
      d22=c-t32vrc

      b11=r11*(d11a-d11b)
      b22=r22*d22
      b12=r12*d21
      b21=r21*s12
      c12=r12*s21
      c21=r21*d12

      dl11dh=k2i*
     &      (b11+b12+b21-s11b-s11b+(r12*r21*d22+c12+c21-d11a-d11b)/r11)
      dl22dh=k2i*
     &      (b12+b21+b22-s22-s22+(r12*r21*(d11a-d11b)+b12+b21+d22)/r22)
      dl12dh=k2i*
     &      (b11+b12+b22+s11a-s11b-s22+(r11*s12+d12)*(1.d0+r22)/r12)
      dl21dh=k2i*
     &      (b11+b21+b22-s11a-s11b-s22+(r11*d21+s21)*(1.d0+r22)/r21)
      RETURN
      END      ! SW_DRV EQU

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC