text:
module mod_init

  use mod_common
  use decomp_2d
  use mod_param
  use mpi
  use mod_common_mpi
  use mod_misc

  implicit none

  private
  public init_program

contains

  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  !
  !

  subroutine init_program(begin, duold,dvold,dwold,pold,curv_restart)

    use mod_loadflds
    use mod_initsolver
    use mod_bound
    use mod_chkdiv
    use mod_param_ls
    use mod_interface
    use mod_ab
    use mod_mom
    use mod_rib
    use mod_contact_line

    integer, intent(in) :: begin

    real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: duold,dvold,dwold
    real, dimension(0:i1,0:j1,0:k1), intent(out) :: pold, curv_restart

    call write_log_header

    !---- Velocity, Pressure, Level Set

    if (begin .eq. 0) then
       time = 0.
       call init_vel_p
       call init_levelset
    else
       call init_vel_p ! still need to re-init the vel in the case of partial inflow
       call loadflds(0,begin)
       if (myid .eq. 0) write(6,*) 'Restarting from step ',begin
    endif

    if (Riblet) then
       call init_static_ls
       call init_rib
    endif

    if (Contact_line) then
       call init_contact_line
    endif

    call init_matrix  ! least squares coef matrix
    call initsolver   ! pressure solver coef (real and Fourier spaces)

    !# before, pp2 => p

    dt = timestep  ! time step

    !---- Boundary Conditions (Physical and Computational)

    call bounduvw(unew,vnew,wnew)
    call bounduvw(dudt,dvdt,dwdt)
    call chkdiv(begin)
    call boundmls(lvset)
    call boundp(pnew)
    if (Riblet) call mask_rib_vel(unew,vnew,wnew)

    call average_vel_at_center  ! cell-center velocity

    !call ls_bound_error  ! check extrap. error

    !---- Phase, Density, Viscosity, and Volume

    call update_materials
    init_volume = volume

    !---- Obtain (n-1) values when (re)starting AB2 (effectively Euler's method)

    call mom_main(duold,dvold,dwold)

    ! only for semi-Lagrangian
    select case (ls_adv_space)
    case('semiLag')
       uccold = uccnew
       vccold = vccnew
       wccold = wccnew
    end select

    select case (surface_tension_method)
    case('CSF')
       call get_curvature(ynew,curv_restart)
       call boundc(curv_restart)
       call surf_ten(curv_restart,surfxold,surfyold,surfzold)

    case('GFM')
       pold = pnew
       p_xold = p_x
       p_yold = p_y
       p_zold = p_z
       dvndn = 0.  ! normal derivative of normal velocity
    end select

    return
  end subroutine init_program

  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! Initialize velocity and pressure
  !
  subroutine init_vel_p
    !
    ! unew,vnew,wnew : Velocity
    ! dudt,dvdt,dwdt : Acceleration
    ! pnew           : Pressure
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    use mod_loadflds

    integer :: i,j,k, iglob,jglob,kglob

    real :: coory,coorz, lowerbound
    real :: v_upper,v_lower,shear_rate
    real :: num, rn1,rn2,rn3
    real :: umean,vmean,wmean
    real :: umean_all,vmean_all,wmean_all
    real :: alpha,beta,pC1,pD1,CsPgrad  ! two-layer flow parameters

    real, dimension(kmax) :: v

    lowerbound =  coords(2)*ly/(dims(2)*1.)

    select case(iniu)

    case('zer')

       forall(i=0:i1,j=0:j1,k=0:k1)
          unew(i,j,k) = 0.
          vnew(i,j,k) = 0.
          wnew(i,j,k) = 0.
          dudt(i,j,k) = 0.
          dvdt(i,j,k) = 0.
          dwdt(i,j,k) = 0.
          pnew(i,j,k) = 0.
       end forall

    case('uni')

       forall(i=0:i1,j=0:j1,k=0:k1)
          unew(i,j,k) = 0.
          vnew(i,j,k) = 1.
          wnew(i,j,k) = 0.
          dudt(i,j,k) = 0.
          dvdt(i,j,k) = 0.
          dwdt(i,j,k) = 0.
          pnew(i,j,k) = 0.
       end forall

    case('poi')

       do k = 0,k1
          coorz = (k-0.5)*dz/lz ! [-dz/2 , 1+dz/2]
          do j = 0,j1
             do i = 0,i1
                unew(i,j,k) = 0.
                vnew(i,j,k) = 6.*v_bulk_desired*coorz*(1.-coorz)  ! max 1.5*bulk
                wnew(i,j,k) = 0.
                dudt(i,j,k) = 0.
                dvdt(i,j,k) = 0.
                dwdt(i,j,k) = 0.
                pnew(i,j,k) = 0.
             enddo
          enddo
       enddo

    case('cou')

       v_upper = 1.5  ! upper plate velocity
       v_lower =-1.5  ! lower plate velocity
       shear_rate = (v_upper - v_lower)/lz

       do k = 0,k1
          coorz = (k-0.5)*dz ! [-dz/2 , lz+dz/2]
          do j = 0,j1
             do i = 0,i1
                unew(i,j,k) = 0.
                vnew(i,j,k) = v_lower + shear_rate*coorz
                wnew(i,j,k) = 0.
                dudt(i,j,k) = 0.
                dvdt(i,j,k) = 0.
                dwdt(i,j,k) = 0.
                pnew(i,j,k) = 0.
             enddo
          enddo
       enddo

    case('sep') !# 2D deformation vortex (Serpentine)

       do k = 0,k1
          coorz = (k-0.5)*dz/lz ! [-dz/2 , 1+dz/2]
          do j = 0,j1
             coory = lowerbound +(j)*dy/lz ! [0 , ly/lz+dy]
             do i = 0,i1
                unew(i,j,k) = 0.
                vnew(i,j,k) = -2. * (sin(pi*coory))**2 * sin(pi*coorz)*cos(pi*coorz)

                dudt(i,j,k) = 0.
                dvdt(i,j,k) = 0.
                dwdt(i,j,k) = 0.
                pnew(i,j,k) = 0.
             enddo
          enddo
       enddo
       do k = 0,k1
          coorz = (k)*dz/lz ! [0 , 1+dz]
          do j = 0,j1
             coory = lowerbound +(j-.5)*dy/lz ! [-dy/2 , ly/lz+dy/2]
             do i = 0,i1

                wnew(i,j,k) = 2. * sin(pi*coory)*cos(pi*coory) * (sin(pi*coorz))**2

             enddo
          enddo
       enddo

    case('rot') ! 2D rigid-body rotation

       do k = 0,k1
          coorz = (k-.5)*dz/lz ! [-dz/2 , 1+dz/2]
          do j = 0,j1
             coory = lowerbound + (j)*dy/lz ! [0 , ly/lz+dy]

             unew(:,j,k) = 0.
             vnew(:,j,k) =  2.*pi*(coorz-.5)
             dudt(:,j,k) = 0.
             dvdt(:,j,k) = 0.
             dwdt(:,j,k) = 0.
             pnew(:,j,k) = 0.
          enddo
       enddo
       do k = 0,k1
          coorz = (k)*dz/lz ! [0 , 1+dz]
          do j = 0,j1
             coory = lowerbound + (j-.5)*dy/lz ! [-dy/2 , ly/lz+dy/2]

             wnew(:,j,k) = -2.*pi*(coory-.5)
          enddo
       enddo

    case('2lc')

       alpha = miu1/miu2
       beta = 0.2  !# interface height or lower layer thickness
       pC1 = - (1.-beta**2*(1.-alpha))/(1.-beta*(1.-alpha))
       pD1 = -1.-pC1
       CsPgrad = -10.*visc  !times visc to ensure velocities are of abuot order 1
       call random_seed ()  !# generating random numbers

       do k=0,k1
          coorz = (k-0.5)/dzi/lz ! normalised with channel height
          do j=0,j
             do i=0,i1
                !# exact solution for two-layer flow
                if (coorz .gt. beta) then
                   vnew(i,j,k) = CsPgrad/miu1/visc*(coorz*coorz+pC1*coorz+pD1)  !# upper layer
                else
                   vnew(i,j,k) = CsPgrad/miu2/visc*(coorz*coorz+pC1*coorz)      !# lower layer
                endif
                call random_number (num)  !#
                unew(i,j,k) = 0.
                wnew(i,j,k) = num*1.e-6  !# small perturbation
                dudt(i,j,k) = 0.
                dvdt(i,j,k) = 0.
                dwdt(i,j,k) = 0.
                pnew(i,j,k) = 0.
             enddo
          enddo
       enddo

       v_bulk = 0.
       do k=1,kmax
          do i=1,itot
             v_bulk = v_bulk+vnew(i,j,k)*dz*dx
          enddo
       enddo
       vnew(:,:,:) = lx*lz/v_bulk*vnew(:,:,:)  !#  make sure v_bulk = 1. initially

    case('log') ! Logarithmic profile

       !## ---- if(isfreeslip) then
       if (.not. ns_np) then
          do k=1,kmax
             v(k) = 2.5*log( (1./visc)*(k-0.5)/(2.*kmax) ) + 5.5
             if ((k-0.5)/kmax/visc .le. 11.6) v(k)=(k-0.5)/kmax/visc
          enddo
       else
          do k=1,kmax/2
             v(k) = 2.5*log( (1./visc)*(k-0.5)/(1.*kmax) ) + 5.5
             if ((k-0.5)/kmax/visc .le. 11.6) v(k)=(k-0.5)/kmax/visc
          enddo
          do k=1,kmax/2
             i = kmax+1-k
             v(kmax+1-k) = 2.5*log( (1./visc)* (k-0.5)/(1.*kmax) ) + 5.5
             if ((k-0.5)/kmax/visc .le. 11.6) v(kmax+1-k)=(k-0.5)/kmax/visc
          enddo
       endif
       !
       ! Below, random numbers are generated such that the initial field
       ! does not depend in the number of mpi tasks.
       !
       umean=0.
       vmean=0.
       wmean=0.
       do kglob=1,ktot
          do jglob=1,jtot
             do iglob=1,itot
                call random_number(rn1)
                call random_number(rn2)
                call random_number(rn3)
                i = iglob-imax*coords(1)
                j = jglob-jmax*coords(2)
                k = kglob
                if( &
                     (1.le.i.and.i.le.imax) .and. &
                     (1.le.j.and.j.le.jmax) &
                     ) then
                   unew(i,j,k)=15.*(rn1-0.5)
                   umean=umean+unew(i,j,k)
                   vnew(i,j,k)=v(k)*(1.+5.*(rn2-0.5))
                   vmean=vmean+vnew(i,j,k)
                   wnew(i,j,k)=15.*(rn3-0.5)
                   wmean=wmean+wnew(i,j,k)
                   dudt(i,j,k)=0.
                   dvdt(i,j,k)=0.
                   dwdt(i,j,k)=0.
                   pnew(i,j,k)=0.
                endif
             enddo
          enddo
       enddo

       call mpi_allreduce(umean,umean_all,1,mpi_real8,mpi_sum,comm_cart,error)
       call mpi_allreduce(vmean,vmean_all,1,mpi_real8,mpi_sum,comm_cart,error)
       call mpi_allreduce(wmean,wmean_all,1,mpi_real8,mpi_sum,comm_cart,error)
       umean_all = umean_all/(1.*itot*jtot*kmax)
       vmean_all = vmean_all/(1.*itot*jtot*kmax)
       wmean_all = wmean_all/(1.*itot*jtot*kmax)
       do k=1,kmax
          do j=1,jmax
             do i=1,imax
                unew(i,j,k) = unew(i,j,k)-umean_all !average = 0
                vnew(i,j,k) = vnew(i,j,k)/vmean_all !average = 1
                wnew(i,j,k) = wnew(i,j,k)-wmean_all !average = 0
             enddo
          enddo
       enddo

    case('fld')  ! load previous fld data

       forall(i=0:i1,j=0:j1,k=0:k1)
          dudt(i,j,k) = 0.
          dvdt(i,j,k) = 0.
          dwdt(i,j,k) = 0.
       end forall

       call init_from_fld  ! initialize unew,vnew,wnew,pnew (only 1:max)

    end select

    !-- Initialize inflow profile (if any) --!

    select case (BC_in_y)
    case('InOut')
       call inflow_profile
    end select

    !-- Initialize pressure jump --!

    p_xold = 0.
    p_yold = 0.
    p_zold = 0.
    p_x = 0.
    p_y = 0.
    p_z = 0.

    return
  end subroutine init_vel_p

  !------------------------------------
  ! Speicify an inflow velocity profile
  !
  subroutine inflow_profile

    integer :: i,j,k, i_glb, num
    integer, parameter :: num_chunk = 1  ! number of inflow chunks
    integer, dimension(num_chunk) :: i_start,i_end, k_start,k_end

    real :: ratio_i,ratio_k
    real :: v_mean, frontbound,coorz

    v_inflow = 0.

    !-- Starting and ending indices for i and k

    i_start(1) = itot/3
    i_end(1)   = itot/3*2

    k_start(1) = 0
    k_end(1)   = 9 !ktot

    !-- Ratios of the non-zero inflow in each Cartesian dir

    ratio_i = 1.0*(i_end(1)-i_start(1))/itot
    ratio_k = 1.0*(k_end(1)-k_start(1))/ktot

    v_mean = v_bulk_desired/(ratio_i*ratio_k)

    frontbound = coords(2)*ly/(1.*dims(2)) ! front boundary of process myid

    if (frontbound .eq. 0.) then
       do k = 1,kmax
          do num = 1,num_chunk
             if (k .gt. k_start(num) .and. k .le. k_end(num)) then ! k-dir is not parallelized
                coorz = (k - 0.5)*dz/(lz*ratio_k)
                do i = 1,imax
                   i_glb = i + coords(1)*imax  ! global i index
                   if (i_glb .gt. i_start(num) .and. i_glb .le. i_end(num)) then
                      v_inflow(i,k) = 6.*v_mean*coorz*(1.-coorz)  ! parabola of max = 1.5 (z from 0 to lz*ratio_k)
                   endif
                enddo ! i loop
             endif
          enddo ! num loop
       enddo! k loop
    endif

    return
  end subroutine inflow_profile

  !~~~~~~~~~~~~~~~~~~~~~
  ! Initialize level set
  !
  subroutine init_levelset
  !
  !     ynew    : Level set function
  !       + for phase 1
  !       - for phase 2
  !
  !     dydt    : No meaning (kept for historical reason)
  !
  !     coor?   : Global coordinates in ? direction (from mpi)
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    integer :: i,j,k,l
    real :: coorx,coory,coorz
    real :: lowerbound,upperbound,leftbound,rightbound

    real :: alpha,beta,gamma
    real :: phi_vc,phi_hc  ! Zalesak's disk
    real :: ypc,zpc,radius, gap,center_dist    ! pancake parameters
    real :: xdp,zdp,ydp, ydp1,ydp2,ydp3,zdp1   ! droplet parameters
    real :: phi_periodic1,phi_periodic2,phi_periodic3
    real :: xs,ys,zs,rs, s_theta,s_phi, sin_the,cos_the,sin_phi,cos_phi    ! Hadamard
    real :: coef_a,coef_b,coef_c,coef_d, vel_r,vel_t

    real, dimension(:), allocatable :: yz_angle, v_angle,h_angle        ! relative angles of droplets/bubbles
    real, dimension(lmax) :: xdp_mls,ydp_mls,zdp_mls,radius_mls  ! multiple level set

    ! sphere packings
    real, dimension(6)  :: x_6pac,y_6pac,z_6pac
    real, dimension(7)  :: x_7pac,y_7pac,z_7pac
    real, dimension(8)  :: x_8pac,y_8pac,z_8pac
    real, dimension(9)  :: x_9pac,y_9pac,z_9pac
    real, dimension(10) :: x_10pac,y_10pac,z_10pac

    real, dimension(-2:i1+2,-2:j1+2,-2:k1+2) :: phi_pancake, phi_slot
    real, dimension(-2:i1+2,-2:j1+2,-2:k1+2) :: phi_drp1, phi_drp2, phi_drp3

    ynew(:,:,:) = 0.
    dydt(:,:,:) = 0.
    volume = 0.

    !---- Bounds of Each Core ----!

    lowerbound =  coords(2)   *ly/(dims(2)*1.)
    upperbound = (coords(2)+1)*ly/(dims(2)*1.)
    leftbound  =  coords(1)   *lx/(dims(1)*1.)
    rightbound = (coords(1)+1)*lx/(dims(1)*1.)

    select case(inil)

    case('pan')

       ! Center (ydp,zdp) and Radius

       ypc = ly*0.5
       zpc = 0.5+0.25*sqrt(3.)
       radius = 0.5

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                ynew(i,j,k) = sqrt((coory-ypc)**2+(coorz-zpc)**2) - radius
                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('2pc')

       xdp = 0.5*lx
       zdp = 0.5*lz
       radius = 0.5
       ydp1= 0.5*ly -(radius + 3.*dy)
       ydp2= 0.5*ly +(radius + 3.*dy)

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                !-- Define 2 pancakes separately

                phi_drp1(i,j,k) = sqrt( (coory-ydp1)**2+(coorz-zdp)**2 ) -radius
                phi_drp2(i,j,k) = sqrt( (coory-ydp2)**2+(coorz-zdp)**2 ) -radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt( (coory-ydp1-ly)**2+(coorz-zdp)**2 )-radius
                phi_periodic2 = sqrt( (coory-ydp2-ly)**2+(coorz-zdp)**2 )-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       do k=-2,k1+2
          do j=-2,j1+2
             do i=-2,i1+2

                !-- Level set as min distance to interface

                if (phi_drp1(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp1(i,j,k)
                if (phi_drp2(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp2(i,j,k)

                if (phi_drp1(i,j,k) .gt. 0. .and. &
                     phi_drp2(i,j,k) .gt. 0.        ) then

                   ynew(i,j,k) = min( phi_drp1(i,j,k),phi_drp2(i,j,k) )

                endif

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

       ! Initialize collision velocity

       where (phi_drp1(0:i1,0:j1,0:k1) .le. 1.5*dy) vnew = .5
       where (phi_drp2(0:i1,0:j1,0:k1) .le. 1.5*dy) vnew =-.5

    case('3pc')

       xdp = .5*lx
       zdp = .5*lz
       radius = 1./6.*lz
       ydp1= .1*ly
       ydp2= ydp1 + 2.*radius+radius/5.
       ydp3= ydp2 + 2.*radius+radius/5.

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                !-- Define 3 pancakes separately

                phi_drp1(i,j,k) = sqrt( (coory-ydp1)**2+(coorz-zdp)**2 ) -radius
                phi_drp2(i,j,k) = sqrt( (coory-ydp2)**2+(coorz-zdp)**2 ) -radius
                phi_drp3(i,j,k) = sqrt( (coory-ydp3)**2+(coorz-zdp)**2 ) -radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt( (coory-ydp1-ly)**2+(coorz-zdp)**2 )-radius
                phi_periodic2 = sqrt( (coory-ydp2-ly)**2+(coorz-zdp)**2 )-radius
                phi_periodic3 = sqrt( (coory-ydp3-ly)**2+(coorz-zdp)**2 )-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)
                if (phi_drp3(i,j,k) .gt. 0.) phi_drp3(i,j,k) = min(phi_drp3(i,j,k), phi_periodic3)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       do k=-2,k1+2
          do j=-2,j1+2
             do i=-2,i1+2

                !-- Level set as min distance to interface

                if (phi_drp1(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp1(i,j,k)
                if (phi_drp2(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp2(i,j,k)
                if (phi_drp3(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp3(i,j,k)

                if (phi_drp1(i,j,k) .gt. 0. .and. &
                     phi_drp2(i,j,k) .gt. 0. .and. &
                     phi_drp3(i,j,k) .gt. 0.        ) then

                   ynew(i,j,k) = min( phi_drp1(i,j,k),phi_drp2(i,j,k),phi_drp3(i,j,k) )

                endif

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('mpc')

       !-- Center of each droplet --!

       radius = 0.5
       xdp_mls = 0.5*lx
       ydp_mls(1) = 0.5*ly
       zdp_mls(1) = 0.25*lz

       gap = 3.*dy
       allocate ( yz_angle(3) )
       yz_angle = (/0., 45., 0./) /180.*pi

       do l = 2,lmax

          ydp_mls(l) = ydp_mls(l-1) + (2.*radius + gap)*cos(yz_angle(l))
          zdp_mls(l) = zdp_mls(l-1) + (2.*radius + gap)*sin(yz_angle(l))

       enddo

       do l = 1,lmax

          !-- Different level set for each pancake

          do k=-2,k1+2
             coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
             do j=-2,j1+2
                coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
                do i=-2,i1+2
                   coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                   !-- Define 3 pancakes separately

                   lvset(i,j,k,l) = sqrt( (coory-ydp_mls(l))**2 + &
                        (coorz-zdp_mls(l))**2   ) -radius

                   !-- Consider periodic BC in y direction

                   phi_periodic1 = sqrt( (coory-ydp_mls(l)-ly)**2+(coorz-zdp_mls(l))**2 )-radius

                   if (lvset(i,j,k,l) .gt. 0.) lvset(i,j,k,l) = min(lvset(i,j,k,l),phi_periodic1)

                enddo
             enddo
          enddo

       enddo

       ynew(:,:,:) =  lvset(:,:,:,1)

    case('drp')

       ! Sphere centered at (xdp,ydp,zdp) with radius

       xdp = 0.5*lx
       ydp = 0.5*ly
       zdp = 0.5*lz
       radius = 0.5

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2
                coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                ynew(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp)**2+(coorz-zdp)**2) - radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt((coorx-xdp)**2+(coory-ydp-ly)**2+(coorz-zdp)**2)-radius

                if (ynew(i,j,k) .gt. 0.) ynew(i,j,k) = min(ynew(i,j,k), phi_periodic1)

                !-- Consider periodic BC in x direction

                phi_periodic1 = sqrt((coorx-xdp-lx)**2+(coory-ydp)**2+(coorz-zdp)**2)-radius

                if (ynew(i,j,k) .gt. 0.) ynew(i,j,k) = min(ynew(i,j,k), phi_periodic1)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('2dp')

       xdp = 0.5*lx
       zdp = 0.5*lz
       radius = 0.5
       ydp1= 0.5*ly -(radius + 1.*dy)
       ydp2= 0.5*ly +(radius + 1.*dy)

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2
                coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                !-- Define 2 spheres separately, differ only in center positions

                phi_drp1(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp1)**2+(coorz-zdp)**2) -radius
                phi_drp2(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp2)**2+(coorz-zdp)**2) -radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt((coorx-xdp)**2+(coory-ydp1-ly)**2+(coorz-zdp)**2)-radius
                phi_periodic2 = sqrt((coorx-xdp)**2+(coory-ydp2-ly)**2+(coorz-zdp)**2)-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)

                !-- Consider periodic BC in x direction

                phi_periodic1 = sqrt((coorx-xdp-lx)**2+(coory-ydp1)**2+(coorz-zdp)**2)-radius
                phi_periodic2 = sqrt((coorx-xdp-lx)**2+(coory-ydp2)**2+(coorz-zdp)**2)-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       do k=-2,k1+2
          do j=-2,j1+2
             do i=-2,i1+2

                !-- Level set as min distance to interface

                if (phi_drp1(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp1(i,j,k)
                if (phi_drp2(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp2(i,j,k)

                if (phi_drp1(i,j,k) .gt. 0. .and. &
                     phi_drp2(i,j,k) .gt. 0.        ) then

                   ynew(i,j,k) = min( phi_drp1(i,j,k),phi_drp2(i,j,k) )

                endif

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

       ! Initialize collision velocity

       !     where (phi_drp1(0:i1,0:j1,0:k1) .le. 1.5*dy) vnew = .5
       !     where (phi_drp2(0:i1,0:j1,0:k1) .le. 1.5*dy) vnew =-.5

    case('3dp')

       xdp = .5*lx
       zdp = .5*lz
       radius = 1./6.*lz
       ydp1= .1*ly
       ydp2= ydp1 + 2.*radius+radius/2.
       ydp3= ydp2 + 2.*radius+radius/2.

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2
                coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                !-- Define 3 spheres separately, differ only in center positions

                phi_drp1(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp1)**2+(coorz-zdp)**2) -radius
                phi_drp2(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp2)**2+(coorz-zdp)**2) -radius
                phi_drp3(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp3)**2+(coorz-zdp)**2) -radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt((coorx-xdp)**2+(coory-ydp1-ly)**2+(coorz-zdp)**2)-radius
                phi_periodic2 = sqrt((coorx-xdp)**2+(coory-ydp2-ly)**2+(coorz-zdp)**2)-radius
                phi_periodic3 = sqrt((coorx-xdp)**2+(coory-ydp3-ly)**2+(coorz-zdp)**2)-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)
                if (phi_drp3(i,j,k) .gt. 0.) phi_drp3(i,j,k) = min(phi_drp3(i,j,k), phi_periodic3)

                !-- Consider periodic BC in x direction

                phi_periodic1 = sqrt((coorx-xdp-lx)**2+(coory-ydp1)**2+(coorz-zdp)**2)-radius
                phi_periodic2 = sqrt((coorx-xdp-lx)**2+(coory-ydp2)**2+(coorz-zdp)**2)-radius
                phi_periodic3 = sqrt((coorx-xdp-lx)**2+(coory-ydp3)**2+(coorz-zdp)**2)-radius

                if (phi_drp1(i,j,k) .gt. 0.) phi_drp1(i,j,k) = min(phi_drp1(i,j,k), phi_periodic1)
                if (phi_drp2(i,j,k) .gt. 0.) phi_drp2(i,j,k) = min(phi_drp2(i,j,k), phi_periodic2)
                if (phi_drp3(i,j,k) .gt. 0.) phi_drp3(i,j,k) = min(phi_drp3(i,j,k), phi_periodic3)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       do k=-2,k1+2
          do j=-2,j1+2
             do i=-2,i1+2

                !-- Level set as min distance to interface

                if (phi_drp1(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp1(i,j,k)
                if (phi_drp2(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp2(i,j,k)
                if (phi_drp3(i,j,k) .le. 0.) ynew(i,j,k) = phi_drp3(i,j,k)

                if (phi_drp1(i,j,k) .gt. 0. .and. &
                    phi_drp2(i,j,k) .gt. 0. .and. &
                    phi_drp3(i,j,k) .gt. 0.        ) then

                   ynew(i,j,k) = min( phi_drp1(i,j,k),phi_drp2(i,j,k),phi_drp3(i,j,k) )

                endif

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('mdp')

       !-- Center of each droplet --!

       radius_mls(1) = 0.5
       gap = 0. !2.0*dx

       xdp_mls(1) = 0.5*lx !-1.0*dx -0.5
       ydp_mls(1) = 1.5 !0.5*ly !-1.0*dy -0.5
       zdp_mls(1) = 0.3*lz !-1.0*dz -0.5

       allocate ( v_angle(3) )
       allocate ( h_angle(3) )
       v_angle = 80./180.*pi !pi/2.
       h_angle = (/0., 60., 75./) /180.*pi

       do l = 2,lmax

          radius_mls(l) = 0.5
          center_dist = radius_mls(l-1) + radius_mls(l) + gap

          xdp_mls(l) = xdp_mls(l-1) + center_dist*sin(v_angle(l))*cos(h_angle(l))
          ydp_mls(l) = ydp_mls(l-1) + center_dist*sin(v_angle(l))*sin(h_angle(l))
          zdp_mls(l) = zdp_mls(l-1) + center_dist*cos(v_angle(l))

       enddo

       !!## -- Up-Down
       !xdp_mls(2) = xdp_mls(1)
       !ydp_mls(2) = ydp_mls(1)
       !zdp_mls(2) = 0.5*lz +1.0*dz +0.5
       !!## --

       !!## -- tetrahedron
       !xdp_mls(4) = (xdp_mls(1)+xdp_mls(2)+xdp_mls(3))/3.
       !ydp_mls(4) = (ydp_mls(1)+ydp_mls(2)+ydp_mls(3))/3.
       !zdp_mls(4) = zdp_mls(3)+sqrt(2./3.)*center_dist
       !!## -- polyhedron 5
       !xdp_mls(5) = (xdp_mls(1)+xdp_mls(2)+xdp_mls(3))/3.
       !ydp_mls(5) = (ydp_mls(1)+ydp_mls(2)+ydp_mls(3))/3.
       !zdp_mls(5) = zdp_mls(3)-sqrt(2./3.)*center_dist
       !!## --

       do l = 1,lmax

          !-- Different level set for each drop

          do k=-2,k1+2
             coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
             do j=-2,j1+2
                coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
                do i=-2,i1+2
                   coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                   !-- Define each drop separately

                   lvset(i,j,k,l) = sqrt( (coorx-xdp_mls(l))**2 + &
                        (coory-ydp_mls(l))**2 + &
                        (coorz-zdp_mls(l))**2   ) -radius_mls(l)

                   !-- Consider periodic BC in y direction

                   phi_periodic1 = sqrt((coorx-xdp_mls(l))**2    + &
                        (coory-ydp_mls(l)-ly)**2 + &
                        (coorz-zdp_mls(l))**2      ) -radius_mls(l)

                   if (lvset(i,j,k,l) .gt. 0.) lvset(i,j,k,l) = min(lvset(i,j,k,l),phi_periodic1)

                   !-- Consider periodic BC in x direction

                   phi_periodic2 = sqrt((coorx-xdp_mls(l)-lx)**2 + &
                        (coory-ydp_mls(l))**2    + &
                        (coorz-zdp_mls(l))**2      ) -radius_mls(l)

                   if (lvset(i,j,k,l) .gt. 0.) lvset(i,j,k,l) = min(lvset(i,j,k,l),phi_periodic2)

                enddo
             enddo
          enddo

       enddo

       ynew(:,:,:) =  lvset(:,:,:,1)

    case('pac') ! Rigid sphere packings of 6-10 drops (no gap)

       ! 6 drops

       xs = 0.6  ! x coordinate shift
       ys = 0.6  ! y coordinate shift
       zs = 0.6  ! z coordinate shift

       x_6pac(1) = xs + 0.0
       y_6pac(1) = ys + 0.0
       z_6pac(1) = zs + 0.0
       x_6pac(2) = xs + 0.5555555555555556
       y_6pac(2) = ys + 1.2830005981991683
       z_6pac(2) = zs + 0.9072184232530289
       x_6pac(3) = xs + 1.0
       y_6pac(3) = ys + 0.0
       z_6pac(3) = zs + 0.0
       x_6pac(4) = xs + 1.3333333333333333
       y_6pac(4) = ys + 0.7698003589195010
       z_6pac(4) = zs + 0.5443310539518174
       x_6pac(5) = xs + 0.5
       y_6pac(5) = ys + 0.8660254037844386
       z_6pac(5) = zs + 0.0
       x_6pac(6) = xs + 0.5
       y_6pac(6) = ys + 0.2886751345948129
       z_6pac(6) = zs + 0.8164965809277260

       ! 7 drops

       xs = 0.6  ! x coordinate shift
       ys = 0.6  ! y coordinate shift
       zs = 0.6  ! z coordinate shift

       x_7pac(1) = xs + 0.0
       y_7pac(1) = ys + 0.0
       z_7pac(1) = zs + 0.0
       x_7pac(2) = xs + 1.0925925925925926
       y_7pac(2) = ys + 0.6949586573578829
       z_7pac(2) = zs + 1.5120307054217148
       x_7pac(3) = xs + 1.0
       y_7pac(3) = ys + 0.0
       z_7pac(3) = zs + 0.0
       x_7pac(4) = xs + 0.5555555555555556
       y_7pac(4) = ys + 1.2830005981991683
       z_7pac(4) = zs + 0.9072184232530289
       x_7pac(5) = xs + 0.5
       y_7pac(5) = ys + 0.8660254037844386
       z_7pac(5) = zs + 0.0
       x_7pac(6) = xs + 1.3333333333333333
       y_7pac(6) = ys + 0.7698003589195010
       z_7pac(6) = zs + 0.5443310539518174
       x_7pac(7) = xs + 0.5
       y_7pac(7) = ys + 0.2886751345948129
       z_7pac(7) = zs + 0.8164965809277260

       ! 8 drops

       xs = 0.6  ! x coordinate shift
       ys = 0.9  ! y coordinate shift
       zs = 1.6  ! z coordinate shift

       x_8pac(1) = xs + 0.0
       y_8pac(1) = ys + 0.0
       z_8pac(1) = zs + 0.0
       x_8pac(2) = xs + 1.3888888888888888
       y_8pac(2) = ys + 0.8018753738744803
       z_8pac(2) = zs + 0.4536092116265145
       x_8pac(3) = xs + 1.3888888888888888
       y_8pac(3) = ys  -0.1603750747748960
       z_8pac(3) = zs  -0.9072184232530289
       x_8pac(4) = xs + 0.5555555555555556
       y_8pac(4) = ys + 1.2830005981991683
       z_8pac(4) = zs  -0.9072184232530289
       x_8pac(5) = xs + 1.0
       y_8pac(5) = ys + 0.0
       z_8pac(5) = zs + 0.0
       x_8pac(6) = xs + 0.5
       y_8pac(6) = ys + 0.8660254037844386
       z_8pac(6) = zs + 0.0
       x_8pac(7) = xs + 0.5
       y_8pac(7) = ys + 0.2886751345948129
       z_8pac(7) = zs  -0.8164965809277260
       x_8pac(8) = xs + 1.3333333333333333
       y_8pac(8) = ys + 0.7698003589195010
       z_8pac(8) = zs  -0.5443310539518174

       ! 9 drops

       xs = 1.1  ! x coordinate shift
       ys = 0.6  ! y coordinate shift
       zs = 0.6  ! z coordinate shift

       x_9pac(1) = xs + 0.0
       y_9pac(1) = ys + 0.0
       z_9pac(1) = zs + 0.0
       x_9pac(2) = xs + 1.0
       y_9pac(2) = ys + 0.0
       z_9pac(2) = zs + 0.0
       x_9pac(3) = xs  -0.5
       y_9pac(3) = ys + 0.8660254037844386
       z_9pac(3) = zs + 0.0
       x_9pac(4) = xs + 1.0
       y_9pac(4) = ys + 1.5396007178390021
       z_9pac(4) = zs + 0.5443310539518174
       x_9pac(5) = xs + 1.5
       y_9pac(5) = ys + 0.8660254037844386
       z_9pac(5) = zs + 0.0
       x_9pac(6) = xs + 0.0
       y_9pac(6) = ys + 1.5396007178390021
       z_9pac(6) = zs + 0.5443310539518174
       x_9pac(7) = xs + 0.0
       y_9pac(7) = ys + 0.5773502691896257
       z_9pac(7) = zs + 0.8164965809277260
       x_9pac(8) = xs + 1.0
       y_9pac(8) = ys + 0.5773502691896257
       z_9pac(8) = zs + 0.8164965809277260
       x_9pac(9) = xs + 0.5
       y_9pac(9) = ys + 0.8660254037844386
       z_9pac(9) = zs + 0.0

       ! 10 drops

       xs = 1.1  ! x coordinate shift
       ys = 0.6  ! y coordinate shift
       zs = 0.6  ! z coordinate shift

       x_10pac(1) = xs + 0.0
       y_10pac(1) = ys + 0.0
       z_10pac(1) = zs + 0.0
       x_10pac(2) = xs + 1.0
       y_10pac(2) = ys + 0.0
       z_10pac(2) = zs + 0.0
       x_10pac(3) = xs  -0.5
       y_10pac(3) = ys + 0.2886751345948129
       z_10pac(3) = zs + 0.8164965809277260
       x_10pac(4) = xs + 0.0
       y_10pac(4) = ys + 0.5773502691896257
       z_10pac(4) = zs + 1.6329931618554521
       x_10pac(5) = xs + 1.5
       y_10pac(5) = ys + 0.2886751345948129
       z_10pac(5) = zs + 0.8164965809277260
       x_10pac(6) = xs + 1.0
       y_10pac(6) = ys + 0.5773502691896257
       z_10pac(6) = zs + 1.6329931618554521
       x_10pac(7) = xs + 0.5
       y_10pac(7) = ys + 0.8660254037844386
       z_10pac(7) = zs + 0.0
       x_10pac(8) = xs + 0.0
       y_10pac(8) = ys + 1.1547005383792515
       z_10pac(8) = zs + 0.8164965809277260
       x_10pac(9) = xs + 1.0
       y_10pac(9) = ys + 1.1547005383792515
       z_10pac(9) = zs + 0.8164965809277260
       x_10pac(10)= xs + 0.5
       y_10pac(10)= ys + 0.2886751345948129
       z_10pac(10)= zs + 0.8164965809277260

       do l = 1,lmax

          if (lmax .eq. 6) then
             xdp_mls(l) = x_6pac(l)
             ydp_mls(l) = y_6pac(l)
             zdp_mls(l) = z_6pac(l)

          elseif (lmax .eq. 7) then
             xdp_mls(l) = x_7pac(l)
             ydp_mls(l) = y_7pac(l)
             zdp_mls(l) = z_7pac(l)

          elseif (lmax .eq. 8) then
             xdp_mls(l) = x_8pac(l)
             ydp_mls(l) = y_8pac(l)
             zdp_mls(l) = z_8pac(l)

          elseif (lmax .eq. 9) then
             xdp_mls(l) = x_9pac(l)
             ydp_mls(l) = y_9pac(l)
             zdp_mls(l) = z_9pac(l)

          elseif (lmax .eq. 10) then
             xdp_mls(l) = x_10pac(l)
             ydp_mls(l) = y_10pac(l)
             zdp_mls(l) = z_10pac(l)
          endif

          !-- Different level set for each drop

          do k=-2,k1+2
             coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
             do j=-2,j1+2
                coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
                do i=-2,i1+2
                   coorx = leftbound + (i-0.5)/dxi ! [-dx*5/2 , lx+dx*5/2]

                   !-- Define each drop separately

                   lvset(i,j,k,l) = sqrt( (coorx-xdp_mls(l))**2 + &
                        (coory-ydp_mls(l))**2 + &
                        (coorz-zdp_mls(l))**2   ) - 0.5

                   !-- Consider periodic BC in y direction

                   phi_periodic1 = sqrt((coorx-xdp_mls(l))**2    + &
                        (coory-ydp_mls(l)-ly)**2 + &
                        (coorz-zdp_mls(l))**2      ) - 0.5

                   if (lvset(i,j,k,l) .gt. 0.) lvset(i,j,k,l) = min(lvset(i,j,k,l),phi_periodic1)

                   !-- Consider periodic BC in x direction

                   phi_periodic2 = sqrt((coorx-xdp_mls(l)-lx)**2 + &
                        (coory-ydp_mls(l))**2    + &
                        (coorz-zdp_mls(l))**2      ) - 0.5

                   if (lvset(i,j,k,l) .gt. 0.) lvset(i,j,k,l) = min(lvset(i,j,k,l),phi_periodic2)

                enddo
             enddo
          enddo

       enddo

       ynew(:,:,:) =  lvset(:,:,:,1)

    case('elp')

       ! Ellipse (distorted)

       coef_a = 4.**2
       coef_b = 2.**2
       ypc = 3.5
       zpc = 2.
       coef_c = 0.1

       do k=-2,k1+2
          coorz = (k-0.5)/dzi - lz/2. ! [lz/2-dz*5/2 , lz/2+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi - ly/2. ! [ly/2-dy*5/2 , ly/2+dy*5/2]
             do i=-2,i1+2

                !-- distortion factor
                coef_d = 1. + ((coory-3.5)**2 + (coorz-2.)**2)*0.05
                !coef_d = 1.
                ynew(i,j,k) = coef_d*(sqrt(coory**2/coef_a + coorz**2/coef_b) - 1.)

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)
       reset = lvset

    case('zal') ! Zalesak's disk

       ypc = ly*0.5
       zpc = lz*0.75
       radius = 0.15*ly

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                phi_pancake(i,j,k) = sqrt((coory-ypc)**2+(coorz-zpc)**2) - radius  ! pancake

                phi_vc = min(abs(coory-0.475),abs(coory-0.525))  ! vertical cut
                phi_hc = min(abs(coorz-0.85),abs(coorz-0.6))   ! horizontal cut
                if (abs(coory-0.5) .le. 0.025 .and. abs(coorz-0.725) .le. .125) then  ! slot
                   phi_slot(i,j,k) = min(phi_vc,phi_hc)
                else
                   if (abs(coory-0.5) .le. 0.025) then
                      phi_slot(i,j,k) = -phi_hc
                   elseif (abs(coorz-0.725) .le. .125) then
                      phi_slot(i,j,k) = -phi_vc
                   else
                      phi_slot(i,j,k) = -sqrt(phi_vc**2+phi_hc**2)
                   endif
                endif

             enddo
          enddo
       enddo

       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                if (phi_pancake(i,j,k) .gt. 0.) ynew(i,j,k) = phi_pancake(i,j,k)
                if (phi_slot(i,j,k) .gt. 0.) ynew(i,j,k) = phi_slot(i,j,k)
                if (phi_pancake(i,j,k) .le. 0. .and. phi_slot(i,j,k) .le. 0.) then
                   ynew(i,j,k) = -min( abs(phi_pancake(i,j,k)),abs(phi_slot(i,j,k)) )
                endif

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('2lc')

       beta = 0.2*lz  !# must be consistent with subroutine init
       do k=-2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2

                !ynew(i,j,k) = coorz-beta ! Interface is a straight line at coorz = beta.
                ynew(i,j,k) = coorz - 0.47
                !ynew(i,j,k) = coory-25.
                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('sin')

       alpha = 0.05  ! amplitude
       beta  = 1.00  ! period
       gamma = 0.50  ! mean value

       do k = -2,k1+2
          coorz = (k-0.5)/dzi
          do j = -2,j1+2
             coory = lowerbound + (j-0.5)/dyi
             do i = -2,i1+2

                ynew(i,j,k) = coorz - alpha*cos(coory*2.*pi/beta) - gamma
                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('any')  ! anything (temporary testing)

       alpha = 0.05  ! amplitude
       beta  = 1.00  ! period
       gamma = 0.50  ! mean value

       do k = -2,k1+2
          coorz = (k-0.5)/dzi
          do j = -2,j1+2
             coory = lowerbound + (j-0.5)/dyi
             do i = -2,i1+2

                !ynew(i,j,k) = coorz - alpha*cos(coory*2.*pi/beta) - gamma  ! wave
                ynew(i,j,k) = 0.5-sqrt((coory-0.5)**2+(coorz-0.5-sqrt(3.)/4.)**2) ! circle
                if (coorz .gt. 1.) ynew(i,j,k) = NarrowBand_2
                phi_pancake(i,j,k) = coorz - gamma  ! straight line
                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       do k = -2,k1+2
          do j = -2,j1+2
             coory = lowerbound + (j-0.5)/dyi
             do i = -2,i1+2

                if (coory .lt. 0.25 .or. coory .gt. 0.75) then
                   ynew(i,j,k) = phi_pancake(i,j,k)
                endif

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

    case('hdm')  ! Hadamard's droplet

       ! Sphere centered at (xdp,ydp,zdp) with radius

       xdp = 0.5*lx
       ydp = 0.5*ly
       zdp = 0.5*lz
       radius = 0.5

       do k=-2,k1+2
          coorz = (k-0.5)*dz ! [-dz*5/2 , lz+dz*5/2]
          do j=-2,j1+2
             coory = lowerbound + (j-0.5)*dy ! [-dy*5/2 , ly+dy*5/2]
             do i=-2,i1+2
                coorx = leftbound + (i-0.5)*dx ! [-dx*5/2 , lx+dx*5/2]

                ynew(i,j,k) = sqrt((coorx-xdp)**2+(coory-ydp)**2+(coorz-zdp)**2) - radius

                !-- Consider periodic BC in y direction

                phi_periodic1 = sqrt((coorx-xdp)**2+(coory-ydp-ly)**2+(coorz-zdp)**2)-radius

                if (ynew(i,j,k) .gt. 0.) ynew(i,j,k) = min(ynew(i,j,k), phi_periodic1)

                !-- Consider periodic BC in x direction

                phi_periodic1 = sqrt((coorx-xdp-lx)**2+(coory-ydp)**2+(coorz-zdp)**2)-radius

                if (ynew(i,j,k) .gt. 0.) ynew(i,j,k) = min(ynew(i,j,k), phi_periodic1)

                dydt(i,j,k) = 0.

             enddo
          enddo
       enddo

       lvset(:,:,:,1) = ynew(:,:,:)

       ! Re-initialize the velocity field based on Hadamard's solution

       coef_d = -(rho2-rho1)*9.81/(12.*miu1+18.*miu2) ! original
       !coef_d = -1./(4.*(1.+3.*miu2/miu1)*radius**2) ! normalized s.t. O(vel)~1
       coef_c = -(3.+2.*miu2/miu1)*coef_d*radius**2
       coef_b = -(2.+3.*miu2/miu1)*coef_d*radius**3
       coef_a = miu2/miu1*coef_d*radius**5

       do k = 0,k1
          coorz = (k-0.5)*dz
          do j = 0,j1
             coory = lowerbound + (j-0.5)*dy
             do i = 0,i1
                coorx = leftbound + (i-0.5)*dx

                xs = coorx-xdp
                ys = coory-ydp
                zs = coorz-zdp
                rs = lvset(i,j,k,1) + radius

                cos_the = ys/rs
                sin_the = sqrt(1.-cos_the**2)  ! always positive
                sin_phi = xs/sqrt(xs**2+zs**2)
                cos_phi = zs/sqrt(xs**2+zs**2)

                if (rs .gt. radius) then  ! outside

                   vel_r = 2.*(coef_a/rs**3 + coef_b/rs)*cos_the
                   vel_t =    (coef_a/rs**3 - coef_b/rs)*sin_the

                else  ! inside

                   vel_r = 2.*(coef_c +    coef_d*rs**2)*cos_the
                   vel_t =-2.*(coef_c + 2.*coef_d*rs**2)*sin_the

                endif

                ! cell-center velocities in the droplet frame (not a bad init for cell-face)

                unew(i,j,k) = (vel_r*sin_the + vel_t*cos_the)*sin_phi
                vnew(i,j,k) =  vel_r*cos_the - vel_t*sin_the - 2.*(coef_c + coef_d*radius**2)
                wnew(i,j,k) = (vel_r*sin_the + vel_t*cos_the)*cos_phi

             enddo
          enddo
       enddo

    end select

    return
  end subroutine init_levelset

  !~~~~~~~~~~~~~~~~~~~~~~~~
  !
  subroutine init_static_ls
  !
  ! Initialize a static level set that is not advected later on.
  ! It might be used as a reference or for visualization.
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    integer :: i,j,k,l
    integer :: l_pil
    parameter (l_pil = 2)

    real :: lowerbound,upperbound,leftbound,rightbound
    real :: ypc,zpc, radius
    real :: coory,coorz
    real, dimension(l_pil) :: ypil,zpil
    real :: phi_vc,phi_hc, w
    real, dimension(-2:i1+2,-2:j1+2,-2:k1+2,l_pil) :: phi_pil
    real :: eps,phi

    character(len=3) :: static_shape
    parameter (static_shape = 'pil')

    !---- Bounds of Each Core ----!

    lowerbound =  coords(2)   *ly/(dims(2)*1.)
    upperbound = (coords(2)+1)*ly/(dims(2)*1.)
    leftbound  =  coords(1)   *lx/(dims(1)*1.)
    rightbound = (coords(1)+1)*lx/(dims(1)*1.)

    !---- Define the static level set ----!

    select case(static_shape)

    case('pan')

       ypc = ly*0.5
       zpc = lz*0.5
       radius = 0.5

       do k = -2,k1+2
          coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
          do j = -2,j1+2
             coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
             do i = -2,i1+2

                fixls(i,j,k) = sqrt((coory-ypc)**2+(coorz-zpc)**2) - radius

             enddo
          enddo
       enddo

    case('sin')

       do k = -2,k1+2
          coorz = (k-0.5)/dzi
          do j = -2,j1+2
             coory = lowerbound + (j-0.5)/dyi
             do i = -2,i1+2

                fixls(i,j,k) = coorz - 0.1*cos(coory*2.*pi)-lz/8.

             enddo
          enddo
       enddo

    case('pil') ! pillars

       do l = 1,l_pil
          ypil(l) = 0. +(l-1)*ly  ! streamwise locations
          zpil(l) = 0.5  ! height of the pillar
       enddo

       w = 0.5 ! width of the pillar

       !---- define each pillar ----!

       do l = 1,l_pil
          do k = -2,k1+2
             coorz = (k-0.5)/dzi ! [-dz*5/2 , lz+dz*5/2]
             do j = -2,j1+2
                coory = lowerbound + (j-0.5)/dyi ! [-dy*5/2 , ly+dy*5/2]
                do i = -2,i1+2

                   phi_vc = min(abs(coory-(ypil(l)-w/2.)),abs(coory-(ypil(l)+w/2.)))  ! vertical cut
                   phi_hc = min(abs(coorz-zpil(l)),abs(coorz))   ! horizontal cut
                   if (abs(coory-ypil(l)) .le. w/2. .and. coorz .le. zpil(l)) then  ! pillar
                      phi_pil(i,j,k,l) = -min(phi_vc,phi_hc)
                   else
                      if (abs(coory-ypil(l)) .le. w/2.) then
                         phi_pil(i,j,k,l) = phi_hc
                      elseif (coorz .le. zpil(l)) then
                         phi_pil(i,j,k,l) = phi_vc
                      else
                         phi_pil(i,j,k,l) = sqrt(phi_vc**2+phi_hc**2)
                      endif
                   endif

                enddo
             enddo
          enddo
       enddo

       fixls = phi_pil(:,:,:,1)
       do l=2,l_pil
          fixls = min(fixls, phi_pil(:,:,:,l))
       enddo

    end select

    !---- Corresponding Heaviside for visualization ----!

    fixhv = Heaviside('4','trig', 1,1.5*dz)  ! the 3rd entry doesn't matter

    return
  end subroutine init_static_ls

  !~~~~~~~~~~~~~~~~~~~~~
  !
  subroutine init_matrix
  !
  ! Pre-compute a least-squares matrix (pseudo_inv_A) for level set/curvature
  ! Ax = b --> x = (pseudo_inv_A)b
  ! Can result in ~90% speedup in get/re_curvature.f90
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    integer :: dist_x, dist_y, dist_z, counter
    real, dimension(7,4) :: Al ! linear coef. matrix
    real, dimension(4,4) :: ATAl, ATAl_inv
    real, dimension(27,10) :: A, Ab ! quadratic coef. matrix
    real, dimension(10,10) :: ATA, ATA_inv, ATAb,ATAb_inv

    !-- Define matrix Al --!

    Al(:,:) = 0.
    Al(:,1) = 1.
    Al(3,2) = -dx
    Al(5,2) =  dx
    Al(2,3) = -dy
    Al(6,3) =  dy
    Al(1,4) = -dz
    Al(7,4) =  dz

    !-- Compute the pseudo inverse of Al --!

    ATAl = matmul(transpose(Al), Al)
    ATAl_inv = inv(ATAl)
    pseudo_inv_Al = matmul(ATAl_inv, transpose(Al))  ! stored in common.f90

    !-- Define matrix A --!

    A(:,1)  = 1.

    counter = 1  ! line counter of matrix A
    do dist_z = -1,1
       do dist_y = -1,1
          do dist_x = -1,1

             A(counter,2) = dx*dist_x
             A(counter,3) = dy*dist_y
             A(counter,4) = dz*dist_z

             counter = counter + 1

          enddo
       enddo
    enddo

    A(:,5) = 0.5 * A(:,2)**2
    A(:,6) = 0.5 * A(:,3)**2
    A(:,7) = 0.5 * A(:,4)**2

    A(:,8)  = A(:,2)*A(:,3)
    A(:,9)  = A(:,2)*A(:,4)
    A(:,10) = A(:,3)*A(:,4)

    !-- Compute the pseudo inverse of A --!

    ATA = matmul(transpose(A), A)
    ATA_inv= inv(ATA)
    pseudo_inv_A = matmul(ATA_inv, transpose(A))  ! stored in common.f90

    !-- Define matrix Ab (sided stencil) --!

    ! lower side

    Ab(:,1)  = 1.

    counter = 1  ! line counter of matrix Ab
    do dist_z = 0,2
       do dist_y = -1,1
          do dist_x = -1,1

             Ab(counter,2) = dx*dist_x
             Ab(counter,3) = dy*dist_y
             Ab(counter,4) = dz*dist_z

             counter = counter + 1

          enddo
       enddo
    enddo

    Ab(:,5) = 0.5 * Ab(:,2)**2
    Ab(:,6) = 0.5 * Ab(:,3)**2
    Ab(:,7) = 0.5 * Ab(:,4)**2

    Ab(:,8)  = Ab(:,2)*Ab(:,3)
    Ab(:,9)  = Ab(:,2)*Ab(:,4)
    Ab(:,10) = Ab(:,3)*Ab(:,4)

    !-- Compute the pseudo inverse of Ab --!

    ATAb = matmul(transpose(Ab), Ab)
    ATAb_inv= inv(ATAb)
    pseudo_inv_Ab_lower = matmul(ATAb_inv, transpose(Ab))  ! stored in common.f90

    ! upper side

    Ab(:,1)  = 1.

    counter = 1  ! line counter of matrix Ab
    do dist_z = -2,0
       do dist_y = -1,1
          do dist_x = -1,1

             Ab(counter,2) = dx*dist_x
             Ab(counter,3) = dy*dist_y
             Ab(counter,4) = dz*dist_z

             counter = counter + 1

          enddo
       enddo
    enddo

    Ab(:,5) = 0.5 * Ab(:,2)**2
    Ab(:,6) = 0.5 * Ab(:,3)**2
    Ab(:,7) = 0.5 * Ab(:,4)**2

    Ab(:,8)  = Ab(:,2)*Ab(:,3)
    Ab(:,9)  = Ab(:,2)*Ab(:,4)
    Ab(:,10) = Ab(:,3)*Ab(:,4)

    !-- Compute the pseudo inverse of Ab --!

    ATAb = matmul(transpose(Ab), Ab)
    ATAb_inv= inv(ATAb)
    pseudo_inv_Ab_upper = matmul(ATAb_inv, transpose(Ab))  ! stored in common.f90

    return
  end subroutine init_matrix

end module mod_init

max_stars_repo_name: GeZhouyang/ICLS-release

max_issues_repo_name: GeZhouyang/ICLS-release

max_forks_repo_name: GeZhouyang/ICLS-release


text:
      SUBROUTINE SB02MT( JOBG, JOBL, FACT, UPLO, N, M, A, LDA, B, LDB,
     $                   Q, LDQ, R, LDR, L, LDL, IPIV, OUFACT, G, LDG,
     $                   IWORK, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To compute the following matrices
C
C                -1
C         G = B*R  *B',
C
C         -          -1
C         A = A - B*R  *L',
C
C         -          -1
C         Q = Q - L*R  *L',
C
C     where A, B, Q, R, L, and G are N-by-N, N-by-M, N-by-N, M-by-M,
C     N-by-M, and N-by-N matrices, respectively, with Q, R and G
C     symmetric matrices.
C
C     When R is well-conditioned with respect to inversion, standard
C     algorithms for solving linear-quadratic optimization problems will
C     then also solve optimization problems with coupling weighting
C     matrix L. Moreover, a gain in efficiency is possible using matrix
C     G in the deflating subspace algorithms (see SLICOT Library routine
C     SB02OD).
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOBG    CHARACTER*1
C             Specifies whether or not the matrix G is to be computed,
C             as follows:
C             = 'G':  Compute G;
C             = 'N':  Do not compute G.
C
C     JOBL    CHARACTER*1
C             Specifies whether or not the matrix L is zero, as follows:
C             = 'Z':  L is zero;
C             = 'N':  L is nonzero.
C
C     FACT    CHARACTER*1
C             Specifies how the matrix R is given (factored or not), as
C             follows:
C             = 'N':  Array R contains the matrix R;
C             = 'C':  Array R contains the Cholesky factor of R;
C             = 'U':  Array R contains the symmetric indefinite UdU' or
C                     LdL' factorization of R.
C
C     UPLO    CHARACTER*1
C             Specifies which triangle of the matrices R and Q (if
C             JOBL = 'N') is stored, as follows:
C             = 'U':  Upper triangle is stored;
C             = 'L':  Lower triangle is stored.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices A, Q, and G, and the number of
C             rows of the matrices B and L.  N >= 0.
C
C     M       (input) INTEGER
C             The order of the matrix R, and the number of columns of
C             the matrices B and L.  M >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, if JOBL = 'N', the leading N-by-N part of this
C             array must contain the matrix A.
C             On exit, if JOBL = 'N', and INFO = 0, the leading N-by-N
C                                                    -          -1
C             part of this array contains the matrix A = A - B*R  L'.
C             If JOBL = 'Z', this array is not referenced.
C
C     LDA     INTEGER
C             The leading dimension of array A.
C             LDA >= MAX(1,N) if JOBL = 'N';
C             LDA >= 1        if JOBL = 'Z'.
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
C             On entry, the leading N-by-M part of this array must
C             contain the matrix B.
C             On exit, if OUFACT = 1, and INFO = 0, the leading N-by-M
C                                                             -1
C             part of this array contains the matrix B*chol(R)  .
C             On exit, B is unchanged if OUFACT = 2 (hence also when
C             FACT = 'U').
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
C             On entry, if JOBL = 'N', the leading N-by-N upper
C             triangular part (if UPLO = 'U') or lower triangular part
C             (if UPLO = 'L') of this array must contain the upper
C             triangular part or lower triangular part, respectively, of
C             the symmetric matrix Q. The stricly lower triangular part
C             (if UPLO = 'U') or stricly upper triangular part (if
C             UPLO = 'L') is not referenced.
C             On exit, if JOBL = 'N' and INFO = 0, the leading N-by-N
C             upper triangular part (if UPLO = 'U') or lower triangular
C             part (if UPLO = 'L') of this array contains the upper
C             triangular part or lower triangular part, respectively, of
C                                  -          -1
C             the symmetric matrix Q = Q - L*R  *L'.
C             If JOBL = 'Z', this array is not referenced.
C
C     LDQ     INTEGER
C             The leading dimension of array Q.
C             LDQ >= MAX(1,N) if JOBL = 'N';
C             LDQ >= 1        if JOBL = 'Z'.
C
C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,M)
C             On entry, if FACT = 'N', the leading M-by-M upper
C             triangular part (if UPLO = 'U') or lower triangular part
C             (if UPLO = 'L') of this array must contain the upper
C             triangular part or lower triangular part, respectively,
C             of the symmetric input weighting matrix R.
C             On entry, if FACT = 'C', the leading M-by-M upper
C             triangular part (if UPLO = 'U') or lower triangular part
C             (if UPLO = 'L') of this array must contain the Cholesky
C             factor of the positive definite input weighting matrix R
C             (as produced by LAPACK routine DPOTRF).
C             On entry, if FACT = 'U', the leading M-by-M upper
C             triangular part (if UPLO = 'U') or lower triangular part
C             (if UPLO = 'L') of this array must contain the factors of
C             the UdU' or LdL' factorization, respectively, of the
C             symmetric indefinite input weighting matrix R (as produced
C             by LAPACK routine DSYTRF).
C             If FACT = 'N', the stricly lower triangular part (if UPLO
C             = 'U') or stricly upper triangular part (if UPLO = 'L') of
C             this array is used as workspace.
C             On exit, if OUFACT = 1, and INFO = 0 (or INFO = M+1),
C             the leading M-by-M upper triangular part (if UPLO = 'U')
C             or lower triangular part (if UPLO = 'L') of this array
C             contains the Cholesky factor of the given input weighting
C             matrix.
C             On exit, if OUFACT = 2, and INFO = 0 (or INFO = M+1),
C             the leading M-by-M upper triangular part (if UPLO = 'U')
C             or lower triangular part (if UPLO = 'L') of this array
C             contains the factors of the UdU' or LdL' factorization,
C             respectively, of the given input weighting matrix.
C             On exit R is unchanged if FACT = 'C' or 'U'.
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= MAX(1,M).
C
C     L       (input/output) DOUBLE PRECISION array, dimension (LDL,M)
C             On entry, if JOBL = 'N', the leading N-by-M part of this
C             array must contain the matrix L.
C             On exit, if JOBL = 'N', OUFACT = 1, and INFO = 0, the
C             leading N-by-M part of this array contains the matrix
C                      -1
C             L*chol(R)  .
C             On exit, L is unchanged if OUFACT = 2 (hence also when
C             FACT = 'U').
C             L is not referenced if JOBL = 'Z'.
C
C     LDL     INTEGER
C             The leading dimension of array L.
C             LDL >= MAX(1,N) if JOBL = 'N';
C             LDL >= 1        if JOBL = 'Z'.
C
C     IPIV    (input/output) INTEGER array, dimension (M)
C             On entry, if FACT = 'U', this array must contain details
C             of the interchanges performed and the block structure of
C             the d factor in the UdU' or LdL' factorization of matrix R
C             (as produced by LAPACK routine DSYTRF).
C             On exit, if OUFACT = 2, this array contains details of
C             the interchanges performed and the block structure of the
C             d factor in the UdU' or LdL' factorization of matrix R,
C             as produced by LAPACK routine DSYTRF.
C             This array is not referenced if FACT = 'C'.
C
C     OUFACT  (output) INTEGER
C             Information about the factorization finally used.
C             OUFACT = 1:  Cholesky factorization of R has been used;
C             OUFACT = 2:  UdU' (if UPLO = 'U') or LdL' (if UPLO = 'L')
C                          factorization of R has been used.
C
C     G       (output) DOUBLE PRECISION array, dimension (LDG,N)
C             If JOBG = 'G', and INFO = 0, the leading N-by-N upper
C             triangular part (if UPLO = 'U') or lower triangular part
C             (if UPLO = 'L') of this array contains the upper
C             triangular part (if UPLO = 'U') or lower triangular part
C                                                                 -1
C             (if UPLO = 'L'), respectively, of the matrix G = B*R  B'.
C             If JOBG = 'N', this array is not referenced.
C
C     LDG     INTEGER
C             The leading dimension of array G.
C             LDG >= MAX(1,N) if JOBG = 'G',
C             LDG >= 1        if JOBG = 'N'.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (M)
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK; if FACT = 'N', DWORK(2) contains the reciprocal
C             condition number of the given matrix R.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= 1              if FACT = 'C';
C             LDWORK >= MAX(2,3*M,N*M) if FACT = 'N';
C             LDWORK >= MAX(1,N*M)     if FACT = 'U'.
C             For optimum performance LDWORK should be larger than 3*M,
C             if FACT = 'N'.
C             The N*M workspace is not needed for FACT = 'N', if matrix
C             R is positive definite.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = i:  if the i-th element (1 <= i <= M) of the d factor is
C                   exactly zero; the UdU' (or LdL') factorization has
C                   been completed, but the block diagonal matrix d is
C                   exactly singular;
C             = M+1:  if the matrix R is numerically singular.
C
C     METHOD
C                            -     -
C     The matrices G, and/or A and Q are evaluated using the given or
C     computed symmetric factorization of R.
C
C     NUMERICAL ASPECTS
C
C     The routine should not be used when R is ill-conditioned.
C
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Algebraic Riccati equation, closed loop system, continuous-time
C     system, discrete-time system, optimal regulator, Schur form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER         FACT, JOBG, JOBL, UPLO
      INTEGER           INFO, LDA, LDB, LDG, LDL, LDQ, LDR, LDWORK, M,
     $                  N, OUFACT
C     .. Array Arguments ..
      INTEGER           IPIV(*), IWORK(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), G(LDG,*),
     $                  L(LDL,*), Q(LDQ,*), R(LDR,*)
C     .. Local Scalars ..
      LOGICAL           LFACTA, LFACTC, LFACTU, LJOBG, LJOBL, LUPLOU
      CHARACTER         TRANS
      INTEGER           I, J, WRKOPT
      DOUBLE PRECISION  EPS, RCOND, RNORM
C     .. External Functions ..
      LOGICAL           LSAME
      DOUBLE PRECISION  DLAMCH, DLANSY
      EXTERNAL          DLAMCH, DLANSY, LSAME
C     .. External Subroutines ..
      EXTERNAL          DCOPY, DGEMM, DGEMV, DPOCON, DPOTRF, DSYCON,
     $                  DSYRK, DSYTRF, DSYTRS, DTRSM, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         INT, MAX
C     .. Executable Statements ..
C
      INFO   = 0
      LJOBG  = LSAME( JOBG, 'G' )
      LJOBL  = LSAME( JOBL, 'N' )
      LFACTC = LSAME( FACT, 'C' )
      LFACTU = LSAME( FACT, 'U' )
      LUPLOU = LSAME( UPLO, 'U' )
      LFACTA = LFACTC.OR.LFACTU
C
C     Test the input scalar arguments.
C
      IF( .NOT.LJOBG .AND. .NOT.LSAME( JOBG, 'N' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LJOBL .AND. .NOT.LSAME( JOBL, 'Z' ) ) THEN
         INFO = -2
      ELSE IF( .NOT.LFACTA .AND. .NOT.LSAME( FACT, 'N' ) ) THEN
         INFO = -3
      ELSE IF( .NOT.LUPLOU .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -4
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( M.LT.0 ) THEN
         INFO = -6
      ELSE IF( ( LDA.LT.1 ) .OR. ( LJOBL .AND. LDA.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -10
      ELSE IF( ( LDQ.LT.1 ) .OR. ( LJOBL .AND. LDQ.LT.N ) ) THEN
         INFO = -12
      ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
         INFO = -14
      ELSE IF( ( LDL.LT.1 ) .OR. ( LJOBL .AND. LDL.LT.N ) ) THEN
         INFO = -16
      ELSE IF( ( LDG.LT.1 ) .OR. ( LJOBG .AND. LDG.LT.N ) ) THEN
         INFO = -20
      ELSE IF( ( LFACTC .AND. LDWORK.LT.1 ) .OR.
     $         ( LFACTU .AND. LDWORK.LT.MAX( 1, N*M ) ) .OR.
     $    ( .NOT.LFACTA .AND. LDWORK.LT.MAX( 2, N*M, 3*M ) ) ) THEN
         INFO = -23
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'SB02MT', -INFO )
         RETURN
      END IF
C
      IF ( LFACTC ) THEN
         OUFACT = 1
      ELSE IF ( LFACTU ) THEN
         OUFACT = 2
      END IF
C
C     Quick return if possible.
C
      IF ( N.EQ.0 .OR. M.EQ.0 .OR. .NOT.( LJOBL.OR.LJOBG ) ) THEN
        DWORK(1) = ONE
        IF ( .NOT.LFACTA ) DWORK(2) = ONE
        RETURN
      END IF
C
C     (Note: Comments in the code beginning "Workspace:" describe the
C     minimal amount of workspace needed at that point in the code,
C     as well as the preferred amount for good performance.
C     NB refers to the optimal block size for the immediately
C     following subroutine, as returned by ILAENV.)
C
      WRKOPT = 1
C
C     Set relative machine precision.
C
      EPS = DLAMCH( 'Epsilon' )
C
      IF ( .NOT.LFACTA ) THEN
C
C        Compute the norm of the matrix R, which is not factored.
C        Then save the given triangle of R in the other strict triangle
C        and the diagonal in the workspace, and try Cholesky
C        factorization.
C        Workspace: need M.
C
         RNORM = DLANSY( '1-norm', UPLO, M, R, LDR, DWORK )
         CALL DCOPY( M, R, LDR+1, DWORK, 1 )
         IF( LUPLOU ) THEN
C
            DO 20 J = 2, M
               CALL DCOPY( J-1, R(1,J), 1, R(J,1), LDR )
   20       CONTINUE
C
         ELSE
C
            DO 40 J = 2, M
               CALL DCOPY( J-1, R(J,1), LDR, R(1,J), 1 )
   40       CONTINUE
C
         END IF
         CALL DPOTRF( UPLO, M, R, LDR, INFO )
         IF( INFO.EQ.0 ) THEN
C
C           Compute the reciprocal of the condition number of R.
C           Workspace: need 3*M.
C
            CALL DPOCON( UPLO, M, R, LDR, RNORM, RCOND, DWORK, IWORK,
     $                   INFO )
C
C           Return if the matrix is singular to working precision.
C
            OUFACT = 1
            DWORK(2) = RCOND
            IF( RCOND.LT.EPS ) THEN
               INFO = M + 1
               RETURN
            END IF
            WRKOPT = MAX( WRKOPT, 3*M )
         ELSE
C
C           Use UdU' or LdL' factorization, first restoring the saved
C           triangle.
C
            CALL DCOPY( M, DWORK, 1, R, LDR+1 )
            IF( LUPLOU ) THEN
C
               DO 60 J = 2, M
                  CALL DCOPY( J-1, R(J,1), LDR, R(1,J), 1 )
   60          CONTINUE
C
            ELSE
C
               DO 80 J = 2, M
                  CALL DCOPY( J-1, R(1,J), 1, R(J,1), LDR )
   80          CONTINUE
C
            END IF
C
C           Compute the UdU' or LdL' factorization.
C           Workspace: need   1,
C                      prefer M*NB.
C
            CALL DSYTRF( UPLO, M, R, LDR, IPIV, DWORK, LDWORK, INFO )
            OUFACT = 2
            IF( INFO.GT.0 ) THEN
               DWORK(2) = ONE
               RETURN
            END IF
            WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
C
C           Compute the reciprocal of the condition number of R.
C           Workspace: need 2*M.
C
            CALL DSYCON( UPLO, M, R, LDR, IPIV, RNORM, RCOND, DWORK,
     $                   IWORK, INFO )
C
C           Return if the matrix is singular to working precision.
C
            DWORK(2) = RCOND
            IF( RCOND.LT.EPS ) THEN
               INFO = M + 1
               RETURN
            END IF
         END IF
      END IF
C
      IF (OUFACT.EQ.1 ) THEN
C
C        Solve positive definite linear system(s).
C
         IF ( LUPLOU ) THEN
            TRANS = 'N'
         ELSE
            TRANS = 'T'
         END IF
C
C        Solve the system X*U = B, overwriting B with X.
C
         CALL DTRSM( 'Right', UPLO, TRANS, 'Non-unit', N, M,
     $               ONE, R, LDR, B, LDB )
C
         IF ( LJOBG ) THEN
C                                      -1
C           Compute the matrix  G = B*R  *B', multiplying X*X' in G.
C
            CALL DSYRK( UPLO, 'No transpose', N, M, ONE, B, LDB, ZERO,
     $                  G, LDG )
         END IF
C
         IF( LJOBL ) THEN
C
C           Update matrices A and Q.
C
C           Solve the system Y*U = L, overwriting L with Y.
C
            CALL DTRSM( 'Right', UPLO, TRANS, 'Non-unit', N, M,
     $                  ONE, R, LDR, L, LDL )
C
C           Compute A <- A - X*Y'.
C
            CALL DGEMM( 'No transpose', 'Transpose', N, N, M, -ONE, B,
     $                  LDB, L, LDL, ONE, A, LDA )
C
C           Compute Q <- Q - Y*Y'.
C
            CALL DSYRK( UPLO, 'No transpose', N, M, -ONE, L, LDL, ONE,
     $                  Q, LDQ )
         END IF
      ELSE
C
C        Solve indefinite linear system(s).
C
C        Solve the system UdU'*X = B' (or LdL'*X = B').
C        Workspace: need N*M.
C
         DO 100 J = 1, M
            CALL DCOPY( N, B(1,J), 1, DWORK(J), M )
  100    CONTINUE
C
         CALL DSYTRS( UPLO, M, N, R, LDR, IPIV, DWORK, M, INFO )
C
         IF ( LJOBG ) THEN
C                                                    -1
C           Compute a triangle of the matrix  G = B*R  *B' = B*X.
C
            IF ( LUPLOU ) THEN
               I = 1
C
               DO 120 J = 1, N
                  CALL DGEMV( 'No transpose', J, M, ONE, B, LDB,
     $                        DWORK(I), 1, ZERO, G(1,J), 1 )
                  I = I + M
  120          CONTINUE
C
            ELSE
C
               DO 140 J = 1, N
                  CALL DGEMV( 'Transpose', M, J, ONE, DWORK, M, B(J,1),
     $                        LDB, ZERO, G(J,1), LDG )
  140          CONTINUE
C
            END IF
         END IF
C
         IF( LJOBL ) THEN
C
C           Update matrices A and Q.
C
C           Solve the system UdU'*Y = L' (or LdL'*Y = L').
C
            DO 160 J = 1, M
               CALL DCOPY( N, L(1,J), 1, DWORK(J), M )
  160       CONTINUE
C
            CALL DSYTRS( UPLO, M, N, R, LDR, IPIV, DWORK, M, INFO )
C
C           A <- A - B*Y.
C
            CALL DGEMM( 'No transpose', 'No transpose', N, N, M, -ONE,
     $                  B, LDB, DWORK, M, ONE, A, LDA )
C                                            -          -1
C           Compute a triangle of the matrix Q = Q - L*R  *L' = Q - L*Y.
C
            IF ( LUPLOU ) THEN
               I = 1
C
               DO 180 J = 1, N
                  CALL DGEMV( 'No transpose', J, M, -ONE, L, LDL,
     $                        DWORK(I), 1, ONE, Q(1,J), 1 )
                  I = I + M
  180          CONTINUE
C
            ELSE
C
               DO 200 J = 1, N
                  CALL DGEMV( 'Transpose', M, J, -ONE, DWORK, M, L(J,1),
     $                        LDL, ONE, Q(J,1), LDQ )
  200          CONTINUE
C
            END IF
         END IF
      END IF
C
      DWORK(1) = WRKOPT
      IF ( .NOT.LFACTA ) DWORK(2) = RCOND
C
C *** Last line of SB02MT ***
      RETURN
      END

ext: f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
      function totdens(r,z)
      common /flags/ idiskflag, ibulgeflag, ihaloflag

      totdens = 0
      psi=pot(r,z)
      psi0=pot(0.,0.)
!      print*,psi,diskdens(r,z,psi),bulgedenspsi(psi),densrpsi(r,psi)
      if( idiskflag .eq. 1 .or. idiskflag.eq.3) then
         totdens = totdens + diskdens(r,z,psi)
      endif
      if( ihaloflag .eq. 1 ) then
         totdens = totdens + densrpsi(r,psi)
      endif
      if( ihaloflag .eq. 3 ) then
         totdens = totdens + halodfdenspsi(psi,psi0)
      endif
      if( ibulgeflag .eq. 1 ) then
         totdens = totdens + bulgedenspsi(psi,psi0)
      endif
      if( ibulgeflag .eq. 3 ) then
         totdens = totdens + bulgedfdenspsi(psi,psi0)
      endif
      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      FUNCTION bnldev(pp,n,idum)
      INTEGER idum,n
      REAL bnldev,pp,PI
CU    USES gammln,ran1
      PARAMETER (PI=3.141592654)
      INTEGER j,nold
      REAL am,em,en,g,oldg,p,pc,pclog,plog,pold,sq,t,y,gammln,ran1
      SAVE nold,pold,pc,plog,pclog,en,oldg
      DATA nold /-1/, pold /-1./
      if(pp.le.0.5)then
        p=pp
      else
        p=1.-pp
      endif
      am=n*p
      if (n.lt.25)then
        bnldev=0.
        do 11 j=1,n
          if(ran1(idum).lt.p)bnldev=bnldev+1.
11      continue
      else if (am.lt.1.) then
        g=exp(-am)
        t=1.
        do 12 j=0,n
          t=t*ran1(idum)
          if (t.lt.g) goto 1
12      continue
        j=n
1       bnldev=j
      else
        if (n.ne.nold) then
          en=n
          oldg=gammln(en+1.)
          nold=n
        endif
        if (p.ne.pold) then
          pc=1.-p
          plog=log(p)
          pclog=log(pc)
          pold=p
        endif
        sq=sqrt(2.*am*pc)
2       y=tan(PI*ran1(idum))
        em=sq*y+am
        if (em.lt.0..or.em.ge.en+1.) goto 2
        em=int(em)
        t=1.2*sq*(1.+y**2)*exp(oldg-gammln(em+1.)-gammln(en-em+1.)+em*
     *plog+(en-em)*pclog)
        if (ran1(idum).gt.t) goto 2
        bnldev=em
      endif
      if (p.ne.pp) bnldev=n-bnldev
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/bnldev.for


text:
C--------------------------------------------------------------------------
        subroutine ftgatp(ibuff,keynam,value,status)

C       Get ASCII Table Parameter
C       test if the keyword is one of the table column definition keywords
C       of an ASCII table. If so, decode it and update the value in the common
C       block

C       ibuff   i sequence number of the data buffer
C       keynam  c name of the keyword
C       value   c value of the keyword
C       status  i returned error status (0=ok)
C
C       written by Wm Pence, HEASARC/GSFC, June 1991

        integer ibuff,status
        character keynam*8,value*70

C-------COMMON BLOCK DEFINITIONS:--------------------------------------------
C       nb = number of file buffers = max. number of FITS file opened at once
C       nf = maximum number of fields allowed in a table
        integer nf,nb
        parameter (nb = 20)
        parameter (nf = 3000)

C       tfield = number of fields in the table
C       tbcol = byte offset in the row of the beginning of the column
C       rowlen = length of one row of the table, in bytes
C       tdtype = integer code representing the datatype of the column
C       trept = the repeat count = number of data values/element in the column
C       tnull = the value used to represent an undefined value in the column
C       tscale = the scale factor for the column
C       tzero = the scaling zero point for the column
C       scount = the total size of the binary table heap (+ gap if any)
C       theap = the starting byte offset for the binary table heap, relative
C               to the start of the binary table data
C       nxheap = the next empty heap location
        integer tfield,tstart,tbcol,rowlen,tdtype,trept,tnull,scount
        integer theap,nxheap
        double precision tscale,tzero
        common/ft0002/tfield(nb),tstart(nb),tbcol(nf),rowlen(nb),
     &  tdtype(nf),trept(nf),tscale(nf),tzero(nf),tnull(nf),scount(nb)
     &  ,theap(nb),nxheap(nb)

C       cnull = character string representing nulls in character columns
C       cform = the Fortran format of the column
        character cnull*16, cform*8
        common/ft0003/cnull(nf),cform(nf)
C-------END OF COMMON BLOCK DEFINITIONS-----------------------------------

        integer nfield,i,c2,bcol,tstat
        character tform*16

        if (status .gt. 0)return
        tstat=status

        if (keynam(1:5) .eq. 'TFORM')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TFORMn keyword
                    status=tstat
                else
C                   get the TFORM character string, without quotes
                    call ftc2s(value,tform,status)
                    if (status .gt. 0)return
                    if  (tform(1:1) .ne. 'A' .and. tform(1:1) .ne. 'I'
     &              .and. tform(1:1) .ne. 'F' .and. tform(1:1) .ne. 'E'
     &              .and. tform(1:1) .ne. 'D')then
                        status=311
                     call ftpmsg('Illegal '//keynam//' format code: '
     &                           //tform)
                        return
                    end if

                    cform(nfield+tstart(ibuff))=tform
C                   set numeric data type code to indicate an ASCII table field
                    tdtype(nfield+tstart(ibuff))=16
C                   set the repeat count to 1
                    trept(nfield+tstart(ibuff))=1
C                   set the TNULL parameter to the width of the field:
                    c2=0
                    do 10 i=2,8
                        if (tform(i:i) .ge. '0' .and. tform(i:i)
     &                     .le. '9')then
                                c2=i
                        else
                                go to 20
                        end if
10                  continue
20                  continue

                    if (status .gt. 0)return
                    if (c2 .eq. 0)then
C                       no explicit field width, so assume width=1 character
                        tnull(nfield+tstart(ibuff))=1
                    else
                        call ftc2ii(tform(2:c2),tnull(nfield+
     &                              tstart(ibuff)),status)
                        if (status .gt. 0)then
C                               error parsing the TFORM value string
                                status=261
           call ftpmsg('Error parsing '//keynam//' field width: '
     &                  //tform)
                        end if
                    end if
                end if
        else if (keynam(1:5) .eq. 'TBCOL')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TBCOLn keyword
                    status=tstat
                else
C                   get the beginning column number
                    call ftc2ii(value,bcol,status)
                     if (status .gt. 0)then
                        call ftpmsg('Error reading value of '//keynam
     &                  //' as an integer: '//value)
                     else
                        tbcol(nfield+tstart(ibuff))=bcol-1
                     end if
                end if
        else if (keynam(1:5) .eq. 'TSCAL')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TSCALn keyword
                    status=tstat
                else
C                   get the scale factor
                    call ftc2dd(value,tscale(nfield+tstart(ibuff)),
     &                          status)
                    if (status .gt. 0)then
                         call ftpmsg('Error reading value of'//keynam
     &                //' as a Double: '//value)
                    end if
                end if
        else if (keynam(1:5) .eq. 'TZERO')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TZEROn keyword
                    status=tstat
                else
C                   get the scaling zero point
                    call ftc2dd(value,tzero(nfield+tstart(ibuff)),
     &                          status)
                    if (status .gt. 0)then
                         call ftpmsg('Error reading value of'//keynam
     &                //' as a Double: '//value)
                    end if
                end if
        else if (keynam(1:5) .eq. 'TNULL')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TNULLn keyword
                    status=tstat
                else
C                   get the Null value flag (character)
                    call ftc2s(value,cnull(nfield+tstart(ibuff)),status)
                    if (status .gt. 0)then
                         call ftpmsg('Error reading value of'//keynam
     &                //' as a character string: '//value)
                    end if
                end if
        end if
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


MODULE HYPERELASTIC

  IMPLICIT NONE
  INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)
  REAL(KIND=DP), PARAMETER :: HALF=.5_DP, MONE=-1._DP
  REAL(KIND=DP), PARAMETER :: ZERO=0._DP, ONE=1._DP, TWO=2._DP, THREE=3._DP
  REAL(KIND=DP), PARAMETER :: FOUR=4._DP
  REAL(KIND=DP), PARAMETER :: THIRD=ONE/THREE, P2THIRD=TWO/THREE, P3HALF=THREE/TWO
  REAL(KIND=DP), PARAMETER :: FOURTH=ONE/FOUR, SIX=6._DP
  REAL(KIND=DP), PARAMETER :: IIMI4(6,6) = RESHAPE(&
                         [ZERO,  ONE,  ONE, ZERO, ZERO, ZERO, &
                           ONE, ZERO,  ONE, ZERO, ZERO, ZERO, &
                           ONE,  ONE, ZERO, ZERO, ZERO, ZERO, &
                          ZERO, ZERO, ZERO, MONE, ZERO, ZERO, &
                          ZERO, ZERO, ZERO, ZERO, MONE, ZERO, &
                          ZERO, ZERO, ZERO, ZERO, ZERO, MONE], SHAPE(IIMI4))

CONTAINS

  SUBROUTINE HYPEREL(NPROPS, PROPS, TEMP, F, NOEL, CMNAME, INCMPFLAG, &
       NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, PK2, CVOL, CISO)
    ! ----------------------------------------------------------------------- !
    ! HYPERELASTIC MATERIAL MODEL
    ! CALLS A UYHPER MODEL AND FORMULATES THE STRESS AND STIFFNESS ARRAYS
    ! ----------------------------------------------------------------------- !
    USE TENSALG, ONLY : DET,INVARS,INV,DYAD,SYMSHUFF,SYMSQ,I6,II1,II5
    ! --- PASSED ARGUMENTS
    CHARACTER*8, INTENT(IN) :: CMNAME
    INTEGER, INTENT(IN) :: NPROPS, NOEL, NSTATV, INCMPFLAG, NFLDV
    REAL(KIND=DP), INTENT(IN) :: PROPS(NPROPS), TEMP, F(3,3)
    REAL(KIND=DP), INTENT(INOUT) :: STATEV(NSTATV), FIELDV(NFLDV), DFIELDV(NFLDV)
    REAL(KIND=DP), INTENT(OUT) :: PK2(6), CVOL(6,6), CISO(6,6)
    ! --- LOCAL VARIABLES
    INTEGER :: I
    REAL(KIND=DP) :: JAC, C(6)
    REAL(KIND=DP) :: U(2), DU(3), D2U(6), D3U(6)
    REAL(KIND=DP) :: I1, I2, I3, I1B, I2B
    REAL(KIND=DP) :: SCALE, CI(6), A(3), DA(3,6), B(3,6), DB(3,6,6)
    REAL(KIND=DP) :: TERM1(6,6), TERM2(6,6)
    REAL(KIND=DP), DIMENSION(6,6) :: CII, ICI, CICI, LCICI, CCI, CIC
    ! ----------------------------------------------------------- HYPEREL --- !

    ! DEFORMATION TENSOR
    ! C = FT.F
    C = SYMSQ(F)

    ! JACOBIAN
    JAC = DET(F)

    ! INVARIANTS OF C
    CALL INVARS(C, I1, I2, I3)

    ! INVARIANTS OF CBAR
    SCALE = JAC ** (-ONE / THREE)
    I1B = I1 * (SCALE ** 2)
    I2B = I2 * (SCALE ** 4)

    CALL UHYPER(I1B, I2B, JAC, U, DU, D2U, D3U, TEMP, NOEL, CMNAME, &
         INCMPFLAG, NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, &
         NPROPS, PROPS)

    ! UPDATE THE STRESS AND MATERIAL STIFFNESS
    ! THE NOTATION BELOWS FOLLOWS THE APPENDIX OF ***TITLE*** WHERE
    ! THE SECOND PIOLA-KIRCHOFF STRESS IS GIVEN BY
    !
    !                           DU   DU  DIB
    !                   PK2 = 2 -- = --- --- = A.B
    !                           DC   DIB DC
    !
    ! WHERE A IS AN ARRAY CONTAINING THE DERIVATIVES OF U WRT ISOCHORIC
    ! INVARIANTS AND B IS AN ARRAY OF TENSORS CONTAINING THE DERIVATIVES OF
    ! THE ISOCHORIC INVARIANTS WRT C

    ! HELPER QUANTITIES
    CI = INV(C)

    A = DU
    B(1,1:6) = SCALE ** 2 * (I6 - THIRD * I1 * CI)
    B(2,1:6) = SCALE ** 4 * (I1 * I6 - C - P2THIRD * I2 * CI)
    B(3,1:6) = HALF * JAC * CI

    ! SECOND PIOLA-KIRCHHOFF STRESS: DERIVATIVE OF ENERGY WRT C
    FORALL(I=1:6) PK2(I) = TWO * SUM(A(1:3) * B(1:3, I))

    ! USING THE NOTATION FROM ABOVE, THE MATERIAL STIFFNESS L IS
    !
    !                      L  = 4(DA.B + A.DB)
    !
    !            DA   DA  DIB   D2U
    ! WHERE DA = -- = --- --- = ----.B
    !            DC   DIB DC    DIB2
    !
    !          DB
    ! AND DB = --
    !          DC
    !
    ! NOTE THAT DA AND B CONTAIN SECOND ORDER TENSORS SO THAT THE ARRAY DA.B
    ! IS AN ARRAY OF DYADIC PRODUCTS OF SECOND ORDER TENSORS (FOURTH ORDER
    ! TENSORS). DB IS AN ARRAY OF FOURTH ORDER TENSORS.

    ! HELPER QUANTITIES
    CICI = DYAD(CI, CI)
    CIC = DYAD(CI, C)
    CCI = DYAD(C, CI)
    LCICI = SYMSHUFF(CI)
    ICI = DYAD(I6, CI)
    CII = DYAD(CI, I6)

    ! DA/DC
    DA=ZERO
    DA(1,1:6) = D2U(1) * B(1,1:6) + D2U(4) * B(2,1:6) + D2U(5) * B(3,1:6)
    DA(2,1:6) = D2U(4) * B(1,1:6) + D2U(2) * B(2,1:6) + D2U(6) * B(3,1:6)
    DA(3,1:6) = D2U(5) * B(1,1:6) + D2U(6) * B(2,1:6) + D2U(3) * B(3,1:6)

    ! DB/DC
    DB = ZERO

    TERM1 = CII + ICI
    TERM2 = I1 * (LCICI + THIRD * CICI)
    DB(1,1:6,1:6) = THIRD * SCALE ** 2 * (-TERM1 + TERM2)

    TERM1 = P3HALF * (II1 - II5) + (CIC + CCI)
    TERM2 = -I1 * (CII + ICI) + I2 * (LCICI + P2THIRD * CICI)
    DB(2,1:6,1:6) = P2THIRD * SCALE ** 4 * (TERM1 + TERM2)

    DB(3,1:6,1:6) = FOURTH * JAC * (CICI - TWO * LCICI)

    ! ISOCHORIC PART
    CISO = ZERO
    DO I=1,2
       TERM1 = FOUR * (DYAD(DA(I,1:6), B(I,1:6)) + DU(I) * DB(I,1:6,1:6))
       CISO = CISO + TERM1
    END DO
    CISO = HALF * (CISO + TRANSPOSE(CISO))

    ! VOLUMETRIC PART
    CVOL = FOUR * (DYAD(DA(3,1:6), B(3,1:6)) + DU(3) * DB(3,1:6,1:6))
    CVOL = HALF * (CVOL + TRANSPOSE(CVOL))

  END SUBROUTINE HYPEREL

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

  SUBROUTINE PUSHALL(F, PK2, CVOL, CISO, CFAC, STRESS, DDSDDE)
    ! ----------------------------------------------------------------------- !
    USE TENSALG, ONLY : PUSH, SYMSHUFF, I6
    ! --- PASSED ARGUMENTS
    REAL(KIND=DP), INTENT(IN) :: F(3,3), PK2(6), CVOL(6,6), CISO(6,6), CFAC(2)
    REAL(KIND=DP), INTENT(OUT) :: STRESS(6), DDSDDE(6,6)
    REAL(KIND=DP) :: TERM1(6,6), TERM2(6,6)
    ! ----------------------------------------------------------- PUSHALL --- !
    STRESS = PUSH(F, PK2)
    TERM1 = (ONE - CFAC(1)) * CISO + (ONE - CFAC(2)) * CVOL
    DDSDDE = PUSH(F, TERM1)

    ! THE CONSTITUTIVE EQUATIONS FOR THE HYPERELASTIC MATERIAL GIVES THE
    ! STRESS IN THE SPATIAL CONFIGURATION (CAUCHY STRESS). HOWEVER, THE
    ! CONSTITITUVE RESPONSE IS IN TERMS OF THE LIE DERIVATIVE, OR TRUESDELL
    ! RATE, OF THE KIRCHOFF STRESS. THE JAUMMANN RATE IS AN ALTERNATIVE LIE
    ! DERIVATIVE, GIVEN BY
    !                              _   .
    !                              S = S + S.W - W.S
    !       _
    ! WHERE S IS THE JAUMMANN RATE OF THE CAUCHY STRESS S. USING THE JAUMMANN
    ! RATE, THE CONSTITUTIVE RESPONSE IN THE COROTATED FRAME IS GIVEN BY
    !                                 _   _
    !                                 S = C:D
    !       _
    ! WHERE C IS THE JAUMMANN TANGENT STIFFNESS TENSOR AND IS DEFINED
    ! IMPLICITY FROM
    !                           _
    !                           S = C:D + D.S + S.D
    ! FROM WHICH
    !           _
    !           C = C + .5 (Iik Sjl + Iil Sjk + Sik Ijl + Sil Ijk)
    TERM1 = SYMSHUFF(I6, STRESS)
    TERM2 = SYMSHUFF(STRESS, I6)
    DDSDDE = DDSDDE + (TERM1 + TERM2)

  END SUBROUTINE PUSHALL

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

  SUBROUTINE NUMERICAL_JACOBIAN(NPROPS, PROPS, TEMP, NOEL, CMNAME, INCMPFLAG, &
       NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, F, DDSDDE)
    ! THIS PROCEDURE NUMERICALLY COMPUTES AND RETURNS THE JACOBIAN MATRIX
    !
    !                      J = (JIJ) = (DSIGI/DEPSJ).
    !
    ! THE COMPONENTS OF JSUB ARE COMPUTED NUMERICALLY USING A CENTERED
    ! DIFFERENCING SCHEME WHICH REQUIRES TWO CALLS TO THE MATERIAL MODEL
    ! SUBROUTINE FOR EACH ELEMENT OF PK2. THE CENTERING IS ABOUT THE POINT
    !                       EPS = EPSOLD + D * DT,
    ! WHERE D IS THE RATE-OF-STRAIN ARRAY.
    USE TENSALG, ONLY: PUSH
    REAL(KIND=DP), PARAMETER :: TOL=1.E-6_DP
    CHARACTER*8, INTENT(IN) :: CMNAME
    INTEGER, INTENT(IN) :: NPROPS, NOEL, NSTATV, INCMPFLAG, NFLDV
    REAL(KIND=DP), INTENT(IN) :: PROPS(NPROPS), TEMP, STATEV(NSTATV)
    REAL(KIND=DP), INTENT(INOUT) :: FIELDV(NFLDV), DFIELDV(NFLDV)
    REAL(KIND=DP), INTENT(IN) :: F(3,3)
    REAL(KIND=DP), INTENT(OUT) :: DDSDDE(6,6)
    INTEGER :: N
    REAL(KIND=DP) :: EPS, DF(3,3), CISO(6,6), CVOL(6,6)
    REAL(KIND=DP) :: SVP(NSTATV), FP(3,3), TP(6), SP(6)
    REAL(KIND=DP) :: SVM(NSTATV), FM(3,3), TM(6), SM(6)
    ! ---------------------------------------------------------- JACOBIAN --- !
    DDSDDE = ZERO
    EPS = SQRT(EPSILON(ONE))
    DO N = 1, 6
       DF = DGEDDG(EPS/TWO, N, F)
       FP = F + DF
       SVP = STATEV
       CALL HYPEREL(NPROPS, PROPS, TEMP, FP, NOEL, CMNAME, INCMPFLAG, NSTATV,&
            SVP, NFLDV, FIELDV, DFIELDV, TP, CVOL, CISO)
       SP = PUSH(FP, TP)
       FM = F - DF
       SVM = STATEV
       CALL HYPEREL(NPROPS, PROPS, TEMP, FM, NOEL, CMNAME, INCMPFLAG, NSTATV,&
            SVM, NFLDV, FIELDV, DFIELDV, TM, CVOL, CISO)
       SM = PUSH(FM, TM)
       DDSDDE(:, N) = (SP - SM) / EPS
    END DO
    DDSDDE = HALF * (DDSDDE + TRANSPOSE(DDSDDE))
  END SUBROUTINE NUMERICAL_JACOBIAN

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

  FUNCTION DGEDDG(EPS, N, F)
    ! -----------------------------------------------------------------------
    !  !
    ! PERTURB THE DEFORMATION GRADIENT
    ! -----------------------------------------------------------------------
    !  !
    USE TENSALG, ONLY: DYAD
    INTEGER, INTENT(IN) :: N
    REAL(KIND=DP), INTENT(IN) :: EPS, F(3,3)
    REAL(KIND=DP) :: DGEDDG(3,3)
    REAL(KIND=DP), PARAMETER :: O=1._DP,Z=0._DP
    REAL(KIND=DP), PARAMETER :: E(3,3)=RESHAPE((/O,Z,Z,Z,O,Z,Z,Z,O/),(/3,3/))
    INTEGER, PARAMETER :: IJ(2,6)=RESHAPE((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/))
    REAL(KIND=DP) :: EI(3), EJ(3), EIJ(3,3), EJI(3,3)
    ! -----------------------------------------------------------------------
    !  !
    EI = E(IJ(1,N),:3)
    EJ = E(IJ(2,N),:3)
    EIJ = DYAD(EI, EJ)
    EJI = DYAD(EJ, EI)
    DGEDDG = EPS / TWO * (MATMUL(EIJ, F) + MATMUL(EJI, F))
  END FUNCTION DGEDDG

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

  SUBROUTINE NEOHOOKE(NPROPS, PROPS, F, SIG, C)
    ! --------------------------------------------------------------------- !
    ! COMPRESSIBLE NEO-HOOKEAN HYPERELASTIC MATERIAL
    !
    ! NOTES
    ! -----
    ! SYMMETRIC TENSOR ORDERING : XX, YY, ZZ, XY, YZ, ZX
    ! --------------------------------------------------------------------- !
    USE TENSALG, ONLY: I6
    INTEGER, INTENT(IN) :: NPROPS
    REAL(DP), INTENT(IN) :: PROPS(NPROPS), F(3,3)
    REAL(DP), INTENT(INOUT) :: SIG(6), C(6,6)
    INTEGER :: I, J
    REAL(DP) :: EE(6), EEP(3), BBP(3), BBN(3,3)
    REAL(DP) :: C10, D1, EG, EK, EG23, PR
    REAL(DP) :: JAC, SCALE, FB(3,3), BB(6), TRBBAR

    ! ELASTIC PROPERTIES
    C10 = PROPS(1)
    D1 = PROPS(2)

    ! JACOBIAN AND DISTORTION TENSOR
    JAC = F(1,1) * F(2,2) * F(3,3) - F(1,2) * F(2,1) * F(3,3) + F(1,2) * F(2&
         &,3) * F(3,1) + F(1,3) * F(3,2) * F(2,1) - F(1,3) * F(3,1) * F(2,2) &
         &- F(2,3) * F(3,2) * F(1,1)
    SCALE = JAC **(-ONE / THREE)
    FB = SCALE * F

    ! DEVIATORIC LEFT CAUCHY-GREEN DEFORMATION TENSOR
    BB(1) = FB(1,1) * FB(1,1) + FB(1,2) * FB(1,2) + FB(1,3) * FB(1,3)
    BB(2) = FB(2,1) * FB(2,1) + FB(2,2) * FB(2,2) + FB(2,3) * FB(2,3)
    BB(3) = FB(3,1) * FB(3,1) + FB(3,2) * FB(3,2) + FB(3,3) * FB(3,3)
    BB(4) = FB(2,1) * FB(1,1) + FB(2,2) * FB(1,2) + FB(2,3) * FB(1,3)
    BB(5) = FB(3,1) * FB(2,1) + FB(3,2) * FB(2,2) + FB(3,3) * FB(2,3)
    BB(6) = FB(3,1) * FB(1,1) + FB(3,2) * FB(1,2) + FB(3,3) * FB(1,3)
    TRBBAR = SUM(BB(1:3)) / THREE
    EG = TWO * C10 / JAC
    EK = TWO / D1 * (TWO * JAC - ONE)
    PR = TWO / D1 * (JAC - ONE)

    ! CAUCHY STRESS
    SIG = EG * (BB - TRBBAR * I6) + PR * I6

    ! SPATIAL STIFFNESS
    EG23 = EG * TWO / THREE
    C(1,1) =  EG23 * (BB(1) + TRBBAR) + EK
    C(1,2) = -EG23 * (BB(1) + BB(2)-TRBBAR) + EK
    C(1,3) = -EG23 * (BB(1) + BB(3)-TRBBAR) + EK
    C(1,4) =  EG23 * BB(4) / TWO
    C(1,5) = -EG23 * BB(5)
    C(1,6) =  EG23 * BB(6) / TWO
    C(2,2) =  EG23 * (BB(2) + TRBBAR) + EK
    C(2,3) = -EG23 * (BB(2) + BB(3)-TRBBAR) + EK
    C(2,4) =  EG23 * BB(4) / TWO
    C(2,5) =  EG23 * BB(5) / TWO
    C(2,6) = -EG23 * BB(6)
    C(3,3) =  EG23 * (BB(3) + TRBBAR) + EK
    C(3,4) = -EG23 * BB(4)
    C(3,5) =  EG23 * BB(5) / TWO
    C(3,6) =  EG23 * BB(6) / TWO
    C(4,4) =  EG * (BB(1) + BB(2)) / TWO
    C(4,5) =  EG * BB(6) / TWO
    C(4,6) =  EG * BB(5) / TWO
    C(5,5) =  EG * (BB(1) + BB(3)) / TWO
    C(5,6) =  EG * BB(4) / TWO
    C(6,6) =  EG * (BB(2) + BB(3)) / TWO
    FORALL(I=1:6,J=1:6,J<I) C(I,J) = C(J,I)
    RETURN
  END SUBROUTINE NEOHOOKE
END MODULE HYPERELASTIC

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

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, RPL,DDSDDT,DRPLDE,DRPLDT,&
     & STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, NDI,NSHR,NTENS&
     &,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, CELENT,F0,F1,NOEL,NPT,LAYER&
     &,KSPT,KSTEP,KINC)
  ! -------------------------------------------------------------------------
  !  !
  !     UMAT INTERFACE FOR FINITE STRAIN PROPELLANT MODEL
  ! -------------------------------------------------------------------------
  !  !
  !     CANNOT BE USED FOR PLANE STRESS
  !     PROPS(1)  - C10
  !     PROPS(2)  - K1
  ! -------------------------------------------------------------------------
  !  !
  USE HYPERELASTIC, ONLY: DP, ZERO, ONE, TWO, THREE, HYPEREL, &
       NUMERICAL_JACOBIAN, NEOHOOKE, PUSHALL
  IMPLICIT NONE

  ! --- PASSED ARGUMENTS
  ! -------------------------------------------------------------------------
  !  !
  CHARACTER*8, INTENT(IN) :: CMNAME
  INTEGER, INTENT(IN) :: NDI, NSHR, NTENS, NSTATV, NPROPS
  INTEGER, INTENT(IN) :: NOEL, NPT, LAYER, KSPT, KSTEP, KINC
  REAL(KIND=DP), INTENT(IN) :: SSE, SPD, SCD, RPL, DRPLDT, TIME(2), DTIME
  REAL(KIND=DP), INTENT(IN) :: TEMP, DTEMP, PNEWDT, CELENT
  REAL(KIND=DP), INTENT(INOUT) :: STRESS(NTENS), STATEV(NSTATV)
  REAL(KIND=DP), INTENT(INOUT) :: DDSDDE(NTENS, NTENS)
  REAL(KIND=DP), INTENT(INOUT) :: DDSDDT(NTENS), DRPLDE(NTENS)
  REAL(KIND=DP), INTENT(IN) :: STRAN(NTENS), DSTRAN(NTENS)
  REAL(KIND=DP), INTENT(IN) :: PREDEF(1), DPRED(1), PROPS(NPROPS), COORDS(3)
  REAL(KIND=DP), INTENT(IN) :: DROT(3, 3), F0(3, 3), F1(3, 3)
  ! --- LOCAL VARIABLES
  INTEGER, PARAMETER :: NFLDV=1, INCMPFLAG=1
  REAL(KIND=DP) :: FIELDV(NFLDV), DFIELDV(NFLDV), ERR
  REAL(KIND=DP) :: PK2(6), CVOL(6,6), CISO(6,6), CFAC(2)
  REAL(KIND=DP), ALLOCATABLE :: SNEO(:), CNEO(:,:), S(:), C(:,:), CNUM(:,:)
  INTEGER :: INTV(1), DBGFLG, I
  CHARACTER*120 :: MSG
  REAL(KIND=DP) :: REALV(1)
  CHARACTER*8 :: CHARV(1)
  ! ---------------------------------------------------------------- UMAT --- !
  DDSDDE = 0._DP
  CALL HYPEREL(NPROPS, PROPS, TEMP, F1, NOEL, CMNAME, INCMPFLAG, &
       NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, PK2, CVOL, CISO)
  CFAC = ZERO
  CALL PUSHALL(F1, PK2, CVOL, CISO, CFAC, STRESS, DDSDDE)

  DBGFLG = 0
  SELECT CASE (DBGFLG)
  CASE (1)
     ! NUMERICAL JACOBIAN
     CALL NUMERICAL_JACOBIAN(NPROPS, PROPS, TEMP, NOEL, CMNAME, INCMPFLAG, &
          NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, F1, DDSDDE)
  CASE (2)
     ALLOCATE(SNEO(6))
     CALL NEOHOOKE(2, (/PROPS(1),PROPS(2)/), F1, SNEO, DDSDDE)
     DEALLOCATE(SNEO)
  CASE (3)

     ALLOCATE(SNEO(6))
     ALLOCATE(CNEO(6,6))
     ALLOCATE(CNUM(6,6))
     ALLOCATE(S(6))
     ALLOCATE(C(6,6))
     C = DDSDDE
     S = STRESS
     WHERE (ABS(C) < EPSILON(ONE))
        C = ONE
     END WHERE
     WHERE (ABS(S) < EPSILON(ONE))
        S = ONE
     END WHERE

     ! NUMERICAL JACOBIAN
     CALL NUMERICAL_JACOBIAN(NPROPS, PROPS, TEMP, NOEL, CMNAME, INCMPFLAG, &
          NSTATV, STATEV, NFLDV, FIELDV, DFIELDV, F1, CNUM)

     ! NEO HOOKE JACOBIAN
     CALL NEOHOOKE(2, (/PROPS(1),PROPS(2)/), F1, SNEO, CNEO)
     ERR = MAXVAL(((SNEO - STRESS) / S))
     IF (ERR > 1.E-8_DP) THEN
        WRITE(MSG, "(A,ES9.3)") "ERROR IN STRESS: ", ERR
        CALL STDB_ABQERR(-1, MSG, INTV, REALV, CHARV)
     END IF

     ERR = MAXVAL(ABS((CNEO - DDSDDE) / C))
     IF (ERR > 1.E-8_DP) THEN
        WRITE(MSG, "(A,ES9.3)") "ERROR IN STIFFNESS: ", ERR
        CALL STDB_ABQERR(-1, MSG, INTV, REALV, CHARV)
     END IF

     DEALLOCATE(SNEO, CNEO, CNUM, S, C)
  END SELECT

END SUBROUTINE UMAT


text:
        subroutine midpoint(a1,a2,a3,b1,b2,b3,m1,m2,m3)
        implicit real*8 (A-Z)

        call normalize(a1,a2,a3,maga)
        call normalize(b1,b2,b3,magb)
        m1 = 0.5d0 * (a1+b1)
        m2 = 0.5d0 * (a2+b2)
        m3 = 0.5d0 * (a3+b3)
        call normalize(m1,m2,m3,magm)

        return
        end

        subroutine anglevectors(a1,a2,a3,b1,b2,b3,angle)
        implicit real*8 (A-Z)

        maga = mag(a1,a2,a3)
        magb = mag(b1,b2,b3)

        call crossproduct(a1,a2,a3,b1,b2,b3,c1,c2,c3)

        magc = mag(c1,c2,c3)

        sinangle = magc/(maga*magb)
        cosangle = dotproduct(a1,a2,a3,b1,b2,b3)/(maga*magb)

        angle = atan2(sinangle,cosangle)

        return
        end

        function dotproduct(a1,a2,a3,b1,b2,b3)
        implicit real*8 (A-Z)

        dotproduct  = a1 * b1 + a2 * b2 + a3 * b3

        return
        end

        function angle_vectors(a1,a2,a3,b1,b2,b3)
        implicit real*8 (A-Z)

        maga = mag(a1,a2,a3)
        magb = mag(b1,b2,b3)

        call crossproduct(a1,a2,a3,b1,b2,b3,c1,c2,c3)

        magc = mag(c1,c2,c3)

        sinangle = magc/(maga*magb)
        cosangle = dotproduct(a1,a2,a3,b1,b2,b3)/(maga*magb)

        angle_vectors = atan2(sinangle,cosangle)

        return
        end

        subroutine normalize(a1,a2,a3,maga)
        implicit real*8 (A-Z)

        maga    = sqrt(a1*a1+a2*a2+a3*a3)

        if(maga .ne. 0.)then

            magainv = 1.d0/maga

            a1 = a1 * magainv
            a2 = a2 * magainv
            a3 = a3 * magainv

        endif

        return
        end

        function mag(a1,a2,a3)
        implicit real*8 (A-Z)

        mag    = sqrt(a1*a1+a2*a2+a3*a3)

        return
        end

        subroutine xyztolatlon(x,y,z,lat,lon)
        implicit real*8 (A-Z)

        lat = atan(z/sqrt(x**2 + y**2))
        lon = atan3(y,x)

        return
        end

        subroutine latlontoxyz(lat,lon,x,y,z)
        implicit real*8 (A-Z)

        x = cos(lon)*cos(lat)
        y = sin(lon)*cos(lat)
        z = sin(lat)

        return
        end

        subroutine crossproduct(a1,a2,a3,b1,b2,b3,c1,c2,c3)
        implicit real*8 (A-Z)

        c1   = a2 * b3 - a3 * b2
        c2   = a3 * b1 - a1 * b3
        c3   = a1 * b2 - a2 * b1

        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense