text:
! Copyright 2019 Pascal Audet
!
! This file is part of PlateFlex.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
!===========================================================================
!
! MODULE conf_cpwt
!
! Configuration module defining global variables used in modules to
! interface with the Python codes.
!
!===========================================================================

      MODULE conf_cpwt

      IMPLICIT NONE

      REAL, PARAMETER :: pi = 3.141592653589793
      INTEGER, PARAMETER :: na = 11
!
! Wavelet parameter
!
      REAL :: k0

      END MODULE conf_cpwt

!===========================================================================
!
! MODULE cpwt
!
! Contains subroutines to calculate the wavelet transform, scalogram and
! cross-spectral quantities with errors
!
!===========================================================================

      MODULE cpwt

      CONTAINS

!---------------------------------------------------------------------------
! Subroutine wlet_transform
!
! Subroutine to compute the complex-valued wavelet transform coefficients
! obtained from directional Morlet wavelets
!---------------------------------------------------------------------------

      SUBROUTINE wlet_transform(grid, nx, ny, dx, dy, kf, ns, wt_grid)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, nnx, nny, ns, nn(2)
      REAL    :: dx, dy
      REAL    :: grid(nx,ny), grid_mirror(2*nx, 2*ny), kf(ns)
      COMPLEX :: wt_grid(nx,ny,na,ns)

      INTEGER :: is, ia
      REAL    :: da, scales, angle, lam, kk
      CHARACTER(LEN=1) :: trail
      CHARACTER(LEN=30) :: progress

! Allocatable work arrays
      REAL, ALLOCATABLE :: grid_pad(:,:), kx(:), ky(:)
      COMPLEX, ALLOCATABLE :: ft_grid(:,:), ft2_grid(:,:), daughter(:,:)
!
! Python bindings
!
!f2py REAL, intent(in) :: grid
!f2py INTEGER, intent(hide),depend(grid) :: nx=shape(grid,0), ny=shape(grid,1)
!f2py REAL, intent(in) :: dx, dy
!f2py REAL, intent(in) :: kf
!f2py INTEGER, intent(hide),depend(kf) :: ns=shape(kf,0)
!f2py COMPLEX, intent(out) :: wt_grid

        CALL npow2(2*nx, nnx)
        CALL npow2(2*ny, nny)
        ALLOCATE(grid_pad(nnx,nny))
!
! Mirror data
!
        grid_mirror = 0.
        grid_mirror(1:nx,1:ny) = grid
        CALL mirror_data(grid_mirror,nx,ny,2*nx,2*ny)
!
! Taper data
!
        ! CALL taper_data(grid_mirror,2*nx,2*ny,6)
!
! Pad with zeros and remove mean value
!
        grid_pad = 0.
        grid_pad(1:2*nx,1:2*ny) = grid_mirror - sum(grid_mirror)/2./2./nx/ny
!
! Complex form
!
        ALLOCATE(ft_grid(nnx,nny), ft2_grid(nnx,nny), daughter(nnx,nny))
        ft_grid = CMPLX(grid_pad)
!
! Fourier transform
!
        nn(1) = nnx
        nn(2) = nny
        CALL fourn(ft_grid,nn,2,1)
!
! Define wavenumbers
!
        ALLOCATE(kx(nnx), ky(nny))
        CALL defk(nnx,nny,dx,dy,kx,ky)
!
! Define angle parameters as per Kirby (2005)
!
        da = 2.*SQRT(-2.*LOG(0.75))/k0
!_______________________________________________
!
!     MAIN WAVELET TANSFORM LOOP
!_______________________________________________
!
!
! Loop through scales
!
        trail = '-'
        progress = ''
        WRITE(*, fmt="(1x,a,i2,a)", advance="no")'#loops = ',ns,':'
        DO is = 1,ns
!
! Define wavenumbers and scales
!
          trail = '_'
          kk = kf(is)*1.e3
          lam = 2.*pi/kk
          scales = k0/kk
          ! WRITE(*, fmt="(1x,a)", advance="no")TRIM(progress)//trail
          WRITE(*, fmt="(1x,i2)", advance="no")is
          progress = TRIM(progress)//trail
!
! Loop through angles
!
          DO ia = 1,na
!
! Define angle increments
!
            angle = REAL(ia-1)*da-pi/2.e0
!
! Calculate daughter wavelet
!
            CALL wave_function(nnx,nny,kx,ky,k0,scales,angle,daughter)
!
! Compute wavelet transform in Fourier space
!
            daughter = daughter*CONJG(ft_grid)
!
! Back to physical space
!
            CALL fourn(daughter,nn,2,-1)
!
! Normalization
!
            ft2_grid = (daughter)/FLOAT(nnx*nny)
!
! Quadrant swap
!
            CALL cfftswap(ft2_grid,nnx,nny)
!
! Store array
!
            wt_grid(:,:,ia,is) = ft2_grid(1:nx,1:ny)

          END DO
        END DO
        WRITE(*,*)

      END SUBROUTINE wlet_transform

!---------------------------------------------------------------------------
! Subroutine wlet_scalogram
!
! Subroutine to compute the real-valued wavelet scalogram obtained from
! azimuthal-averaging of wavelet transform coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_scalogram(wl_trans, nx, ny, ns, wl_sg, ewl_sg)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans(nx,ny,na,ns)
      REAL    :: wl_sg(nx,ny,ns), ewl_sg(nx,ny,ns)
      REAL    :: sg(nx,ny,na,ns)

      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans
!f2py INTEGER, intent(hide),depend(wl_trans) :: nx=shape(wl_trans,0)
!f2py INTEGER, intent(hide),depend(wl_trans) :: ny=shape(wl_trans,1)
!f2py INTEGER, intent(hide),depend(wl_trans) :: ns=shape(wl_trans,3)
!f2py REAL, intent(out) :: wl_sg, ewl_sg

        wl_sg = 0.e0
        DO is = 1,ns
          DO ia = 1,na
            sg(:,:,ia,is) = CABS(wl_trans(:,:,ia,is)*CONJG(wl_trans(:,:,ia,is)))
!
! Azimuthal average
!
            wl_sg(:,:,is) = wl_sg(:,:,is) + sg(:,:,ia,is)
          END DO
        END DO
!
! Normalization
!
        wl_sg = wl_sg/FLOAT(na)
!
! Jackknife estimation
!
        CALL jackknife_1d_sgram(sg,nx,ny,na,ns,ewl_sg)

        RETURN
      END SUBROUTINE wlet_scalogram

!---------------------------------------------------------------------------
! Subroutine wlet_xscalogram
!
! Subroutine to compute the complex-valued wavelet cross-scalogram obtained
! from azimuthal-averaging of wavelet transform coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_xscalogram(wl_trans1, wl_trans2, nx, ny, ns, wl_xsg, ewl_xsg)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans1(nx,ny,na,ns), wl_trans2(nx,ny,na,ns)
      COMPLEX :: wl_xsg(nx,ny,ns), xsg(nx,ny,na,ns)
      REAL :: ewl_xsg(nx,ny,ns), ewl_xsg_real(nx,ny,ns), ewl_xsg_imag(nx,ny,ns)

      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans1, wl_trans2
!f2py INTEGER, intent(hide),depend(wl_trans1) :: nx=shape(wl_trans1,0)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ny=shape(wl_trans1,1)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ns=shape(wl_trans1,3)
!f2py REAL, intent(out) :: wl_xsg, ewl_xsg

        wl_xsg = CMPLX(0.e0,0.e0)
        DO is = 1,ns
          DO ia = 1, na
            xsg(:,:,ia,is) = wl_trans1(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is))
!
! Azimuthal average
!
            wl_xsg(:,:,is) = wl_xsg(:,:,is) + xsg(:,:,ia,is)
          END DO
        END DO
!
! Normalization
!
        wl_xsg = wl_xsg/FLOAT(na)
!
! Jackknife estimation
!
        CALL jackknife_1d_sgram(REAL(xsg),nx,ny,na,ns,ewl_xsg_real)
        CALL jackknife_1d_sgram(IMAG(xsg),nx,ny,na,ns,ewl_xsg_imag)
        ewl_xsg = SQRT(ewl_xsg_real**2 + ewl_xsg_imag**2)

        RETURN
      END SUBROUTINE wlet_xscalogram

!---------------------------------------------------------------------------
! Subroutine wlet_admit_coh
!
! Subroutine to compute the complex-valued wavelet admittance and real-valued
! wavelet coherence obtained from azimuthal-averaging of wavelet transform
! coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_admit_coh(wl_trans1, wl_trans2, nx, ny, ns, &
        wl_admit, ewl_admit, wl_coh, ewl_coh)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans1(nx,ny,na,ns), wl_trans2(nx,ny,na,ns)
      REAL    :: wl_admit(nx,ny,ns), wl_coh(nx,ny,ns)
      REAL    :: ewl_admit(nx,ny,ns), ewl_coh(nx,ny,ns)

      REAL    :: sg1(nx,ny,na,ns), sg2(nx,ny,na,ns)
      COMPLEX :: xsg(nx,ny,na,ns)

      REAL    :: wl_sg1(nx,ny,ns), wl_sg2(nx,ny,ns)
      COMPLEX :: wl_xsg(nx,ny,ns)
      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans1, wl_trans2
!f2py INTEGER, intent(hide),depend(wl_trans1) :: nx=shape(wl_trans1,0)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ny=shape(wl_trans1,1)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ns=shape(wl_trans1,3)
!f2py REAL, intent(out) :: wl_admit, ewl_admit, wl_coh, ewl_coh

        DO is = 1,ns
          DO ia = 1, na

            sg1(:,:,ia,is) = CABS(wl_trans1(:,:,ia,is)*CONJG(wl_trans1(:,:,ia,is)))
            sg2(:,:,ia,is) = CABS(wl_trans2(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is)))
            xsg(:,:,ia,is) = wl_trans1(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is))
!
! Azimuthal average
!
            wl_sg1(:,:,is) = wl_sg1(:,:,is) + sg1(:,:,ia,is)
            wl_sg2(:,:,is) = wl_sg2(:,:,is) + sg2(:,:,ia,is)
            wl_xsg(:,:,is) = wl_xsg(:,:,is) + xsg(:,:,ia,is)

          END DO
        END DO
!
! Normalization
!
        wl_sg1 = wl_sg1/FLOAT(na)
        wl_sg2 = wl_sg2/FLOAT(na)
        wl_xsg = wl_xsg/FLOAT(na)
!
! Admittance and coherence (squared-real coherency)
!
        wl_admit = REAL(wl_xsg)/wl_sg1
        wl_coh = REAL(wl_xsg)**2/(wl_sg1*wl_sg2)
        ! wl_coh = wl_xsg/SQRT(wl_sg1*wl_sg2)
!
! Jackknife estimation
!
        PRINT*,'Calculation jackknife error on admittance and coherence'
        CALL jackknife_admit_coh(xsg,sg1,sg2,nx,ny,na,ns,ewl_admit,ewl_coh)

        RETURN
      END SUBROUTINE wlet_admit_coh

      END MODULE cpwt

max_stars_repo_name: CHEN-Zhaohui/geoist

max_issues_repo_name: CHEN-Zhaohui/geoist

max_forks_repo_name: CHEN-Zhaohui/geoist


text:
      double precision function dgamma (x)
c jan 1984 edition.  w. fullerton, c3, los alamos scientific lab.
C     double precision x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax,
C    1  xmin, y, d9lgmc, dcsevl, d1mach, dexp, dint, dlog,
C    2  dsin, dsqrt
      double precision x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax,
     1  xmin, y, d9lgmc, dcsevl
      double precision temp
C     external d1mach, d9lgmc, dcsevl, dexp, dint, dlog, dsin, dsqrt,
C    1  initds
      external d9lgmc, dcsevl, initds

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
c series for gam        on the interval  0.          to  1.00000e+00
c                                        with weighted error   5.79e-32
c                                         log weighted error  31.24
c                               significant figures required  30.00
c                                    decimal places required  32.05
c
      data gam cs(  1) / +.8571195590 9893314219 2006239994 2 d-2      /
      data gam cs(  2) / +.4415381324 8410067571 9131577165 2 d-2      /
      data gam cs(  3) / +.5685043681 5993633786 3266458878 9 d-1      /
      data gam cs(  4) / -.4219835396 4185605010 1250018662 4 d-2      /
      data gam cs(  5) / +.1326808181 2124602205 8400679635 2 d-2      /
      data gam cs(  6) / -.1893024529 7988804325 2394702388 6 d-3      /
      data gam cs(  7) / +.3606925327 4412452565 7808221722 5 d-4      /
      data gam cs(  8) / -.6056761904 4608642184 8554829036 5 d-5      /
      data gam cs(  9) / +.1055829546 3022833447 3182350909 3 d-5      /
      data gam cs( 10) / -.1811967365 5423840482 9185589116 6 d-6      /
      data gam cs( 11) / +.3117724964 7153222777 9025459316 9 d-7      /
      data gam cs( 12) / -.5354219639 0196871408 7408102434 7 d-8      /
      data gam cs( 13) / +.9193275519 8595889468 8778682594 0 d-9      /
      data gam cs( 14) / -.1577941280 2883397617 6742327395 3 d-9      /
      data gam cs( 15) / +.2707980622 9349545432 6654043308 9 d-10     /
      data gam cs( 16) / -.4646818653 8257301440 8166105893 3 d-11     /
      data gam cs( 17) / +.7973350192 0074196564 6076717535 9 d-12     /
      data gam cs( 18) / -.1368078209 8309160257 9949917230 9 d-12     /
      data gam cs( 19) / +.2347319486 5638006572 3347177168 8 d-13     /
      data gam cs( 20) / -.4027432614 9490669327 6657053469 9 d-14     /
      data gam cs( 21) / +.6910051747 3721009121 3833697525 7 d-15     /
      data gam cs( 22) / -.1185584500 2219929070 5238712619 2 d-15     /
      data gam cs( 23) / +.2034148542 4963739552 0102605193 2 d-16     /
      data gam cs( 24) / -.3490054341 7174058492 7401294910 8 d-17     /
      data gam cs( 25) / +.5987993856 4853055671 3505106602 6 d-18     /
      data gam cs( 26) / -.1027378057 8722280744 9006977843 1 d-18     /
      data gam cs( 27) / +.1762702816 0605298249 4275966074 8 d-19     /
      data gam cs( 28) / -.3024320653 7353062609 5877211204 2 d-20     /
      data gam cs( 29) / +.5188914660 2183978397 1783355050 6 d-21     /
      data gam cs( 30) / -.8902770842 4565766924 4925160106 6 d-22     /
      data gam cs( 31) / +.1527474068 4933426022 7459689130 6 d-22     /
      data gam cs( 32) / -.2620731256 1873629002 5732833279 9 d-23     /
      data gam cs( 33) / +.4496464047 8305386703 3104657066 6 d-24     /
      data gam cs( 34) / -.7714712731 3368779117 0390152533 3 d-25     /
      data gam cs( 35) / +.1323635453 1260440364 8657271466 6 d-25     /
      data gam cs( 36) / -.2270999412 9429288167 0231381333 3 d-26     /
      data gam cs( 37) / +.3896418998 0039914493 2081663999 9 d-27     /
      data gam cs( 38) / -.6685198115 1259533277 9212799999 9 d-28     /
      data gam cs( 39) / +.1146998663 1400243843 4761386666 6 d-28     /
      data gam cs( 40) / -.1967938586 3451346772 9510399999 9 d-29     /
      data gam cs( 41) / +.3376448816 5853380903 3489066666 6 d-30     /
      data gam cs( 42) / -.5793070335 7821357846 2549333333 3 d-31     /
c
      data pi / 3.1415926535 8979323846 2643383279 50 d0 /
c sq2pil is 0.5*alog(2*pi) = alog(sqrt(2*pi))
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
      data ngam, xmin, xmax, xsml, dxrel / 0, 4*0.d0 /
c
      if (ngam.ne.0) go to 10
C     ngam = initds (gamcs, 42, 0.1*sngl(  d1mach) )
      ngam = initds (gamcs, 42, 0.1*sngl(  EPSMIN ) )
c
      call d9gaml (xmin, xmax)
      if (IGAMMA .ne. 0) return
C     xsml = dexp (dmax1 (dlog(d1mach(1)), -dlog(d1mach(2)))+0.01d0)
      xsml =  exp ( max  ( log( FLTMIN  ), - log( FLTMAX  ))+0.01d0)
C     dxrel = dsqrt (d1mach(4))
      dxrel =  sqrt (  EPSMAX )
c
C10   y = dabs(x)
 10   y =  abs(x)
      if (y.gt.10.d0) go to 50
c
c compute gamma(x) for -xbnd .le. x .le. xbnd.  reduce interval and find
c gamma(1+y) for 0.0 .le. y .lt. 1.0 first of all.
c
      n = int(x)
      if (x.lt.0.d0) n = n - 1
      y = x - dble(float(n))
      n = n - 1
C     dgamma = 0.9375d0 + dcsevl (2.d0*y-1.d0, gamcs, ngam)
      temp = dcsevl (2.d0*y-1.d0, gamcs, ngam)
      if (IGAMMA .ne. 0) return
      dgamma = 0.9375d0 + temp
      if (n.eq.0) return
c
      if (n.gt.0) go to 30
c
c compute gamma(x) for x .lt. 1.0
c
      n = -n

C     if (x.eq.0.d0) call seteru (14hdgamma  x is 0, 14, 4, 2)
C     if (x.lt.0.0d0 .and. x+dble(float(n-2)).eq.0.d0) call seteru (
C    1  31hdgamma  x is a negative integer, 31, 4, 2)
C     if (x.lt.(-0.5d0) .and. dabs((x-dint(x-0.5d0))/x).lt.dxrel) call
C    1  seteru (68hdgamma  answer lt half precision because x too near n
C    2egative integer, 68, 1, 1)
C     if (y.lt.xsml) call seteru (
C    1  54hdgamma  x is so close to 0.0 that the result overflows,
C    2  54, 5, 2)

      if (x.eq.0.d0) then
C       write(6,*) 'dgamma : x is 0'
        IGAMMA = 11
        return
      end if

      if (x.lt.0.0d0 .and. x+dble(float(n-2)).eq.0.d0) then
C       write( 6, *) 'dgamma : x is a negative integer'
        IGAMMA = 12
        return
      end if

      if (x.lt.(-0.5d0) .and. abs((x-dble(int(x-0.5d0)))/x).lt.dxrel)
C    1  write(6,*) 'dgamma : answer lt half precision because
C    2                       x too near a negative integer'
     *  JGAMMA = 11

      if (y.lt.xsml) then
c       write(6,*)  'dgamma :,
c    1               x is so close to 0.0 that the result overflows'
        IGAMMA = 13
        return
      end if
c
      do 20 i=1,n
        dgamma = dgamma/(x+dble(float(i-1)) )
 20   continue
      return
c
c gamma(x) for x .ge. 2.0 and x .le. 10.0
c
 30   do 40 i=1,n
        dgamma = (y+dble(float(i))) * dgamma
 40   continue
      return
c
c gamma(x) for dabs(x) .gt. 10.0.  recall y = dabs(x).
c
C50   if (x.gt.xmax) call seteru (32hdgamma  x so big gamma overflows,
C    1  32, 3, 2)

 50   if (x.gt.xmax) then
c       write(6,*) 'dgamma : x so big gamma overflows'
        IGAMMA = 14
        return
      end if
c
      dgamma = 0.d0
C     if (x.lt.xmin) call seteru (35hdgamma  x so small gamma underflows
C    1  , 35, 2, 0)
C     if (x.lt.xmin) return

      if (x.lt.xmin) then
c       write(6,*) 'dgamma : x so small gamma underflows'
        JGAMMA = 12
        return
      end if
c
C     dgamma = dexp ((y-0.5d0)*dlog(y) - y + sq2pil + d9lgmc(y) )
      temp = d9lgmc(y)
      if (IGAMMA .ne. 0) return
      dgamma =  exp ((y-0.5d0)* log(y) - y + sq2pil + temp)
      if (x.gt.0.d0) return
c
C     if (dabs((x-dint(x-0.5d0))/x).lt.dxrel) call seteru (
C    1  61hdgamma  answer lt half precision, x too near negative integer
C    2  , 61, 1, 1)

      if (abs((x-dble(int(x-0.5d0)))/x).lt.dxrel) JGAMMA = 11
c
C     sinpiy = dsin (pi*y)
      sinpiy =  sin (pi*y)
C     if (sinpiy.eq.0.d0) call seteru (
C    1  31hdgamma  x is a negative integer, 31, 4, 2)

      if (sinpiy.eq.0.d0) then
C       write(6,*) 'dgamma : x is a negative integer'
        IGAMMA = 12
        return
      end if
c
      dgamma = -pi/(y*sinpiy*dgamma)
c
      return
      end
      double precision function dgamr (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c this routine, not dgamma(x), should be the fundamental one.
c
C     double precision x, alngx, sgngx, dgamma, dint, dexp, d1mach
      double precision x, alngx, sgngx, dgamma

C     external dexp, dgamma, dint, d1mach
      external dgamma

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
      dgamr = 0.0d0
C     if (x.le.0.0d0 .and. dint(x).eq.x) return
      if (x.le.0.0d0 .and. dble(int(x)).eq.x) return
c
C     call entsrc (irold, 1)
      if (dabs(x).gt.10.0d0) go to 10
C     dgamr = 1.0d0/dgamma(x)
C     call erroff
C     call entsrc (ir, irold)
      temp = dgamma(x)
      if (IGAMMA .ne. 0) then
C       dgamr = d1mach(2)
        dgamr = FLTMAX
        return
      end if
      dgamr = 1.0d0/temp
      return
c
 10   call dlgams (x, alngx, sgngx)
      if (IGAMMA .ne. 0) return
C     call erroff
C     call entsrc (ir, irold)
C     dgamr = sgngx * dexp(-alngx)
      dgamr = sgngx *  exp(-alngx)
      return
c
      end
      subroutine dlgams (x, dlgam, sgngam)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate log abs (gamma(x)) and return the sign of gamma(x) in sgngam.
c sgngam is either +1.0 or -1.0.
c
C     double precision x, dlgam, sgngam, dint, dlngam
      double precision x, dlgam, sgngam, dlngam
C     external dint, dlngam
      external dlngam

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
      dlgam = dlngam(x)
      if (IGAMMA .ne. 0) return
      sgngam = 1.0d0
      if (x.gt.0.d0) return
c
C     int = dmod (-dint(x), 2.0d0) + 0.1d0
C     if (int.eq.0) sgngam = -1.0d0
      intx =  mod (-dble(int(x)), 2.0d0) + 0.1d0
      if (intx.eq.0) sgngam = -1.0d0
c
      return
      end
      function initds (dos, nos, eta)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
c
c initialize the double precision orthogonal series dos so that initds
c is the number of terms needed to insure the error is no larger than
c eta.  ordinarily eta will be chosen to be one-tenth machine precision.
c
c             input arguments --
c dos    dble prec array of nos coefficients in an orthogonal series.
c nos    number of coefficients in dos.
c eta    requested accuracy of series.
c
      double precision dos(nos)

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
C     if (nos.lt.1) call seteru (
C    1  35hinitds  number of coefficients lt 1, 35, 2, 2)
      if (nos.lt.1) JGAMMA = 31
c
      err = 0.
      do 10 ii=1,nos
        i = nos + 1 - ii
        err = err + abs(sngl(dos(i)))
        if (err.gt.eta) go to 20
 10   continue
c
C20   if (i.eq.nos) call seteru (28hinitds  eta may be too small, 28,
C    1  1, 2)
 20   continue
C     if (i.eq.nos) write(6,*) 'initds : eta may be too small'
      if (i.eq.nos) JGAMMA = 32
      initds = i
c
      return
      end
      subroutine d9gaml (xmin, xmax)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
c
c calculate the minimum and maximum legal bounds for x in gamma(x).
c xmin and xmax are not the only bounds, but they are the only non-
c trivial ones to calculate.
c
c             output arguments --
c xmin   dble prec minimum legal value of x in gamma(x).  any smaller
c        value of x might result in underflow.
c xmax   dble prec maximum legal value of x in gamma(x).  any larger
c        value of x might cause overflow.
c
C     double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach,
C    1  dlog
      double precision xmin, xmax, alnbig, alnsml, xln, xold
C     external d1mach, dlog

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
C     alnsml = dlog(d1mach(1))
      alnsml =  log( FLTMIN  )
      xmin = -alnsml
      do 10 i=1,10
        xold = xmin
C       xln = dlog(xmin)
        xln =  log(xmin)
        xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml)
     1    / (xmin*xln+0.5d0)
C       if (dabs(xmin-xold).lt.0.005d0) go to 20
        if ( abs(xmin-xold).lt.0.005d0) go to 20
 10   continue
C     call seteru (27hd9gaml  unable to find xmin, 27, 1, 2)
C     write(6,*) 'd9gaml : unable to find xmin'
      IGAMMA = 21
      return

c
 20   xmin = -xmin + 0.01d0
c
C     alnbig = dlog (d1mach(2))
      alnbig =  log ( FLTMAX  )
      xmax = alnbig
      do 30 i=1,10
        xold = xmax
C       xln = dlog(xmax)
        xln =  log(xmax)
        xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig)
     1    / (xmax*xln-0.5d0)
C       if (dabs(xmax-xold).lt.0.005d0) go to 40
        if ( abs(xmax-xold).lt.0.005d0) go to 40
 30   continue
C     call seteru (27hd9gaml  unable to find xmax, 27, 2, 2)
C     write(6,*) 'd9gaml : unable to find xmax'
      IGAMMA = 22
      return
c
 40   xmax = xmax - 0.01d0
      xmin = dmax1 (xmin, -xmax+1.d0)
c
      return
      end
      double precision function d9lgmc (x)
c august 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c compute the log gamma correction factor for x .ge. 10. so that
c dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x)
c
C     double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach,
C    1  dexp, dlog, dsqrt
      double precision x, algmcs(15), xbig, xmax, dcsevl
      double precision temp
C     external d1mach, dcsevl, dexp, dlog, dsqrt, initds
      external dcsevl, initds

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
c series for algm       on the interval  0.          to  1.00000e-02
c                                        with weighted error   1.28e-31
c                                         log weighted error  30.89
c                               significant figures required  29.81
c                                    decimal places required  31.48
c
      data algmcs(  1) / +.1666389480 4518632472 0572965082 2 d+0      /
      data algmcs(  2) / -.1384948176 0675638407 3298605913 5 d-4      /
      data algmcs(  3) / +.9810825646 9247294261 5717154748 7 d-8      /
      data algmcs(  4) / -.1809129475 5724941942 6330626671 9 d-10     /
      data algmcs(  5) / +.6221098041 8926052271 2601554341 6 d-13     /
      data algmcs(  6) / -.3399615005 4177219443 0333059966 6 d-15     /
      data algmcs(  7) / +.2683181998 4826987489 5753884666 6 d-17     /
      data algmcs(  8) / -.2868042435 3346432841 4462239999 9 d-19     /
      data algmcs(  9) / +.3962837061 0464348036 7930666666 6 d-21     /
      data algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23     /
      data algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24     /
      data algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26     /
      data algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27     /
      data algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29     /
      data algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30     /
c
      data nalgm, xbig, xmax / 0, 2*0.d0 /
c
      if (nalgm.ne.0) go to 10
C     nalgm = initds (algmcs, 15, sngl(d1mach(3)) )
      nalgm = initds (algmcs, 15, sngl(  EPSMIN ) )
C     xbig = 1.0d0/dsqrt(d1mach(3))
      xbig = 1.0d0/ sqrt(  EPSMIN )
C     xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1))))
      xmax =  exp ( min ( log(FLTMAX   /12.d0), - log(12.d0*FLTMIN   )))
c
C10   if (x.lt.10.d0) call seteru (23hd9lgmc  x must be ge 10, 23, 1, 2)
c
 10   if (x.lt.10.d0) then
c       write(6,*) 'd9lgmc : x must be ge 10'
        IGAMMA = 51
C       d9lgmc = d1mach(2)
        d9lgmc = FLTMAX
        return
      end if

      if (x.ge.xmax) go to 20
c
      d9lgmc = 1.d0/(12.d0*x)
C     if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs,
C    1  nalgm) / x

      if (x.lt.xbig) then
        temp   = dcsevl(2.0d0*(10.d0/x)**2-1.d0, algmcs, nalgm)
        if (IGAMMA .ne. 0) then
C         d9gmlc = d1mach(2)
          d9gmlc = FLTMAX
        else
          d9lgmc = temp / x
        end if
      end if
      return
c
 20   d9lgmc = 0.d0
C     call seteru (34hd9lgmc  x so big d9lgmc underflows, 34, 2, 0)
c     write(6,*) 'd9lgmc : x so big d9lgmc underflows'
      JGAMMA = 51
      return
c
      end
      double precision function dcsevl (x, a, n)
c
c evaluate the n-term chebyshev series a at x.  adapted from
c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).
c
c             input arguments --
c x      dble prec value at which the series is to be evaluated.
c a      dble prec array of n terms of a chebyshev series.  in eval-
c        uating a, only half the first coef is summed.
c n      number of terms in array a.
c
      double precision a(n), x, twox, b0, b1, b2
C     double precision d1mach
C     external         d1mach

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/

c
C     if (n.lt.1) call seteru (28hdcsevl  number of terms le 0, 28, 2,2)
C     if (n.gt.1000) call seteru (31hdcsevl  number of terms gt 1000,
C    1  31, 3, 2)
C     if (x.lt.(-1.1d0) .or. x.gt.1.1d0) call seteru (
C    1  25hdcsevl  x outside (-1,+1), 25, 1, 1)
c
      if (n.lt.1) then
C       write(6,*) 'dcsevl : number of terms le 0'
        IGAMMA = 41
C       dcsevl = d1mach(2)
        dcsevl = FLTMAX
        return
      end if

      if (n.gt.1000) then
C       write(6,*) 'dcsevl : number of terms gt 1000'
        IGAMMA = 42
C       dcsevl = d1mach(2)
        dcsevl = FLTMAX
        return
      end if

      if (x.lt.(-1.1d0) .or. x.gt.1.1d0) then
C       write(6,*) 'dcsevl : x outside (-1,+1)'
        IGAMMA = 43
C       dcsevl = d1mach(2)
        dcsevl = FLTMAX
        return
      end if
c
      twox = 2.0d0*x
      b1 = 0.d0
      b0 = 0.d0
      do 10 i=1,n
        b2 = b1
        b1 = b0
        ni = n - i + 1
        b0 = twox*b1 - b2 + a(ni)
 10   continue
c
      dcsevl = 0.5d0 * (b0-b2)
c
      return
      end
      double precision function dlngam (x)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
C     double precision x, dxrel, pi, sinpiy, sqpi2l, sq2pil,
C    1  y, xmax, dint, dgamma, d9lgmc, d1mach, dlog, dsin, dsqrt
      double precision x, dxrel, pi, sinpiy, sqpi2l, sq2pil,
     1  y, xmax, dgamma, d9lgmc
      double precision   temp
C     external d1mach, d9lgmc, dgamma, dint, dlog, dsin, dsqrt
      external d9lgmc, dgamma

      double precision   FLTMIN, FLTMAX, EPSMIN, EPSMAX
      common /MACHFD/    FLTMIN, FLTMAX, EPSMIN, EPSMAX
      save   /MACHFD/

      integer            IGAMMA, JGAMMA
      common /GAMMFD/    IGAMMA, JGAMMA
      save   /GAMMFD/
c
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
c sq2pil = alog (sqrt(2*pi)),  sqpi2l = alog(sqrt(pi/2))
      data sqpi2l / +.2257913526 4472743236 3097614947 441 d+0    /
      data pi / 3.1415926535 8979323846 2643383279 50 d0 /
c
      data xmax, dxrel / 2*0.d0 /
c
      if (xmax.ne.0.d0) go to 10
C     xmax = d1mach(2)/dlog(d1mach(2))
      xmax =  FLTMAX  / log( FLTMAX  )
C     dxrel = dsqrt (d1mach(4))
      dxrel =  sqrt ( FLTMAX  )
c
C10   y = dabs (x)
 10   y =  abs (x)
      if (y.gt.10.d0) go to 20
c
c dlog (dabs (dgamma(x)) ) for dabs(x) .le. 10.0
c
C     dlngam = dlog (dabs (dgamma(x)) )
      temp   = dgamma(x)
      if (IGAMMA .ne. 0) then
C       dlngam = d1mach(2)
        dlngam = FLTMAX
        return
      end if
      dlngam = log (abs (temp) )
      return
c
c dlog ( dabs (dgamma(x)) ) for dabs(x) .gt. 10.0
c
C20   if (y.gt.xmax) call seteru (
C    1  39hdlngam  dabs(x) so big dlngam overflows, 39, 2, 2)

 20   if (y.gt.xmax) then
c       write(6,*) 'dlngam : abs(x) so big dlngam overflows'
        IGAMMA = 61
C       dlngam = d1mach(2)
        dlngam = FLTMAX
        return
      end if
c
C     if (x.gt.0.d0) dlngam = sq2pil + (x-0.5d0)*dlog(x) - x + d9lgmc(y)

      temp = d9lgmc(y)
      if (IGAMMA .ne. 0) then
C       dlngam = d1mach(2)
        dlngam = FLTMAX
        return
      end if

      if (x.gt.0.d0) dlngam = sq2pil + (x-0.5d0)*log(x) - x + temp
      if (x.gt.0.d0) return
c
C     sinpiy = dabs (dsin(pi*y))
      sinpiy =  abs ( sin(pi*y))
C     if (sinpiy.eq.0.d0) call seteru (
C    1  31hdlngam  x is a negative integer, 31, 3, 2)

      if (sinpiy.eq.0.d0) then
c       write(6,*) 'dlngam : x is a negative integer'
        IGAMMA = 62
C       dlngam = d1mach(2)
        dlngam = FLTMAX
        return
      end if
c
C     dlngam = sqpi2l + (x-0.5d0)*dlog(y) - x - dlog(sinpiy) - d9lgmc(y)

      temp = d9lgmc(y)
      if (IGAMMA .ne. 0) then
C       dlngam = d1mach(2)
        dlngam = FLTMAX
        return
      end if

      dlngam = sqpi2l + (x-0.5d0)*log(y) - x - log(sinpiy) - temp
c
C     if (dabs((x-dint(x-0.5d0))*dlngam/x).lt.dxrel) call seteru (
C    1  68hdlngam  answer lt half precision because x too near negative
C    2integer, 68, 1, 1)
      if ( abs((x-dble(int(x-0.5d0)))*dlngam/x).lt.dxrel) JGAMMA = 61

      return
c
      end

ext: f

max_stars_repo_name: PrescienceLab/RPS

max_issues_repo_name: PrescienceLab/RPS

max_forks_repo_name: PrescienceLab/RPS


text:
program Unidades
implicit none

integer:: x

!STATUS='new' crea un nuevo archivo dentro de la carpeta donde estamos trabajando
! si dentro de la carpeta ya exite manda error, porque no lo puede guardar
!STATUS='old' revisa y trabaja con el archivo, si el archivo no existe manda error
!STATUS='UNKNOWN' revisa si existe y crea uno nuevo si es que no y si existe lo reescribe
!STATUS='SCRATCH' CREA un archivo temporal que se borra al momento que se cierra tu programa

open(unit=10,file='Unidades.txt', Action='write',status='replace')
open(unit=20,file='funcion.txt', Action='read',status='old')

read(20,*) x

write(10,*)x

pause
end program Unidades

max_stars_repo_name: Angelpacman/fundamentos-de-fortran

max_issues_repo_name: Angelpacman/fundamentos-de-fortran

max_forks_repo_name: Angelpacman/fundamentos-de-fortran


text:
module mod_gfm

  use mod_param

  implicit none

  private
  public gfm4p, gfm4diff

contains

  !------------------------------------------------------
  !-- The ghost fluid method (GFM) for the computation of
  !-- surface tenseion across a sharp interface.
  !
  subroutine gfm4p(p_jump)

    use mod_common
    use mod_interface
    use mod_bound
    use mod_common_mpi

    integer :: i,j,k,l, nx,ny,nz,flag, cnt

    real :: nondi, factor,idvndn, t1,tmax
    real :: icurv, jump_x, jump_y, jump_z
    real :: cmax,cmin, cmax_all,cmin_all !#

    real, dimension(-1:1) :: ph,cv,ps, jp_x, jp_y, jp_z

    real, dimension(0:i1,0:j1,0:k1) :: curvature, p_exjump, cux,cuy,cuz, re_curv !#
    real, dimension(0:i1,0:j1,0:k1,1:lmax) :: px_sep,py_sep,pz_sep, pj_sep

    real, dimension(0:i1,0:j1,0:k1), intent(out) :: p_jump

    nondi = 1./We  ! non-dimensional surface tension coef

    select case (NS_diff_discretization)
    case('del')
       factor = 0.

    case('GFM')
       factor = 2.*visc*(miu2-miu1)
       if (myid .eq. 0) then
          write(*,*) 'This method seems unrobust. Use at your own risk.'
          write(*,*) 'Program aborted. (Comment these if you insist.)'
       endif
       call mpi_finalize(error)
       stop
    end select

    select case (surface_tension_method)
    case('CSF') ! add nothing
       p_jump = 0.

    case('GFM')

       call ipot  ! get minimal distances between interfaces

       !-- Update jump at time level (n-1) --!

       p_xold = p_x
       p_yold = p_y
       p_zold = p_z

       !-- Compute jump at time level (n) --!

       px_sep = 0.
       py_sep = 0.
       pz_sep = 0.
       pj_sep = 0.
       p_jump = 0.
       p_exjump = 0.
       p_x = 0.
       p_y = 0.
       p_z = 0.

       do l=1,lmax

          !-- Obtain cell-center curvature of each drop --!

          call get_curvature(lvset(:,:,:,l),curv_cmn)
          call boundc(curv_cmn)
          !#
          if (.not. TwoD) then
             call re_curvature(lvset(:,:,:,l),curv_cmn,re_curv)
             curv_cmn = re_curv
             call boundc(curv_cmn)
          endif
          !#

          do k=1,kmax
             do j=1,jmax
                do i=1,imax

                   if ( abs(lvset(i,j,k,l)) .le. 2.*dz) then ! jump only occurs across cut cells

                      jump_x = 0.
                      jump_y = 0.
                      jump_z = 0.
                      jp_x = 0.
                      jp_y = 0.
                      jp_z = 0.

                      !-- x component --!

                      do nx = -1,1,2

                         ph(-1:1) =    lvset(i-1:i+1,j,k,l)
                         cv(-1:1) = curv_cmn(i-1:i+1,j,k)

                         if ( ph(0)*ph(nx) .lt. 0.) then

                            call icurvature(nx, ph,cv, icurv)

                            !-- jump due to surface tension

                            jp_x(nx) = nondi*(-icurv)

                            !-- jump due to discontinuous viscosity

                            !#idvndn = ( dvndn(i,j,k)*abs(ph(nx)) + dvndn(i+nx,j,k)*abs(ph(0)) ) &
                            !#     /( abs(ph(0)) + abs(ph(nx)) )
                            !#
                            !#jump_x = jump_x + factor*idvndn

                            if (nx .eq. 1) px_sep(i,j,k,l) = jp_x(nx)

                         endif
                      enddo

                      jump_x = sum(jp_x)

                      !-- y component --!

                      do ny = -1,1,2

                         ph(-1:1) =    lvset(i,j-1:j+1,k,l)
                         cv(-1:1) = curv_cmn(i,j-1:j+1,k)

                         if ( ph(0)*ph(ny) .lt. 0.) then

                            call icurvature(ny, ph,cv, icurv)

                            !-- jump due to surface tension

                            jp_y(ny) = nondi*(-icurv)

                            !-- jump due to discontinuous viscosity

                            !#idvndn = ( dvndn(i,j,k)*abs(ph(ny)) + dvndn(i,j+ny,k)*abs(ph(0)) ) &
                            !#     /( abs(ph(0)) + abs(ph(ny)) )
                            !#
                            !#jump_y = jump_y + factor*idvndn

                            if (ny .eq. 1) py_sep(i,j,k,l) = jp_y(ny)

                         endif
                      enddo

                      jump_y = sum(jp_y)

                      !-- z component --!

                      do nz = -1,1,2

                         ph(-1:1) =    lvset(i,j,k-1:k+1,l)
                         cv(-1:1) = curv_cmn(i,j,k-1:k+1)

                         if ( ph(0)*ph(nz) .lt. 0.) then

                            call icurvature(nz, ph,cv, icurv)

                            !-- jump due to surface tension

                            jp_z(nz) = nondi*(-icurv)

                            !-- jump due to discontinuous viscosity

                            !#idvndn = ( dvndn(i,j,k)*abs(ph(nz)) + dvndn(i,j,k+nz)*abs(ph(0)) ) &
                            !#     /( abs(ph(0)) + abs(ph(nz)) )
                            !#
                            !#jump_z = jump_z + factor*idvndn

                            if (nz .eq. 1) pz_sep(i,j,k,l) = jp_z(nz)

                         endif
                      enddo

                      jump_z = sum(jp_z)

                      !-- "Laplacian" of pressure jump (separately) --!

                      if (lvset(i,j,k,l) .gt. 0.) then  ! jump higher

                         pj_sep(i,j,k,l) =   (jump_x/dx**2 + jump_y/dy**2 + jump_z/dz**2)

                      elseif (lvset(i,j,k,l) .lt. 0.) then  ! drop lower

                         pj_sep(i,j,k,l) = - (jump_x/dx**2 + jump_y/dy**2 + jump_z/dz**2)
                         px_sep(i,j,k,l) = - px_sep(i,j,k,l)
                         py_sep(i,j,k,l) = - py_sep(i,j,k,l)
                         pz_sep(i,j,k,l) = - pz_sep(i,j,k,l)

                      else
                         pj_sep(i,j,k,l) = 0.
                         px_sep(i,j,k,l) = 0.
                         py_sep(i,j,k,l) = 0.
                         pz_sep(i,j,k,l) = 0.
                         write(*,*) 'Happens to be at the interface (highly unlikely)!'
                      endif

                      !-- Account for close neighbors --!

                      if (l .eq. 1) then  ! initialize

                         p_jump(i,j,k) = pj_sep(i,j,k,l)
                         p_x(i,j,k) = px_sep(i,j,k,l)
                         p_y(i,j,k) = py_sep(i,j,k,l)
                         p_z(i,j,k) = pz_sep(i,j,k,l)

                      else ! superimpose each separate jump

                         p_jump(i,j,k) = p_jump(i,j,k) + pj_sep(i,j,k,l)
                         p_x(i,j,k) = p_x(i,j,k) + px_sep(i,j,k,l)
                         p_y(i,j,k) = p_y(i,j,k) + py_sep(i,j,k,l)
                         p_z(i,j,k) = p_z(i,j,k) + pz_sep(i,j,k,l)

                      endif

                   endif ! jump shell
                enddo ! i
             enddo ! j
          enddo ! k
       enddo ! l

       !!###########
       !
       !      cmax = maxval(cuz(:,:,:))
       !      call mpi_allreduce(cmax,cmax_all,1,mpi_real8,mpi_max,comm_cart,error)
       !      cmin = minval(cuz(:,:,:))
       !      call mpi_allreduce(cmin,cmin_all,1,mpi_real8,mpi_min,comm_cart,error)
       !
       !
       !      if (myid .eq. 2) then
       !        open(60,file=datadir//'curv.txt',position='append')
       !        write(60,'(2ES12.4)') cmax_all,cmin_all
       !        close(60)
       !      endif
       !
       !!###########

       !-- Extral jump due to depletion --!

       if (attractive) then
          call exjump(p_exjump)
          p_jump = p_jump + p_exjump
       endif
    end select

    return
  end subroutine gfm4p

  !-----------------------------------------------
  !-- A hydrodynamic model for the depletion force
  !-- using MLS and GFM
  !
  subroutine exjump(p_jump)

    use mod_common
    use mod_common_mpi
    use mod_bound

    integer :: i,j,k,l, m,n, nx,ny,nz
    real :: phi, phx,phy,phz
    real :: jump_x, jump_y, jump_z

    real :: del_r, dvol,vol_ex, pot_tot, cp, c0,c1, osp

    real, dimension(-1:1) :: fo, pf, jp_x, jp_y, jp_z
    real, dimension(0:i1,0:j1,0:k1) :: phim,phin, flago, pfactor, px_osp,py_osp,pz_osp

    !-- Model Parameters --!

    real, parameter :: rsm = 2.0*dz    ! radius of surfactant micelles
    real, parameter :: rs1 = 0.1*rsm   ! reduce osp when less than this distance
    real, parameter :: rs2 = 0.0*rsm   ! negative osp when less than this distance
    real, parameter :: osp0 = -40.     ! osmotic pressure due to surfactant micelles

    real, dimension(0:i1,0:j1,0:k1), intent(out) :: p_jump

    px_osp = 0.
    py_osp = 0.
    pz_osp = 0.
    p_jump = 0.

    dvol = dx*dy*dz

    do m = 1,lmax-1
       phim = lvset(0:i1,0:j1,0:k1,m)  ! droplet m

       do n = m+1,lmax
          phin = lvset(0:i1,0:j1,0:k1,n)  ! droplet n

          !-- Reduce magnitude if too close --!

          if (dmin(m,n) .gt. rs1) then
             osp = osp0
          else
             del_r = (rs1-dmin(m,n))/(rs1-rs2)
             osp = osp0*(1.-del_r)
          endif

          !-- Identify overlap of surfactant shell --!

          where (phim .le. rsm .and. phin .le. rsm)
             flago = 1.  ! inside overlap
          elsewhere
             flago = 0.  ! outside overlap
          end where

          call boundc(flago)

          !-- Compute overlap volume --!

          vol_ex = sum( flago(1:imax,1:jmax,1:kmax) )*dvol  ! a rough count
          call mpi_allreduce(mpi_in_place,vol_ex,1,mpi_real8,mpi_sum,comm_cart,error)

          pot_tot = osp*vol_ex  ! total potential energy

          !-- Compute the spatially varying osmotic pressure --!

          pfactor = 0.

          do k=1,kmax
             do j=1,jmax
                do i=1,imax
                   if (flago(i,j,k) .eq. 1.) then

                      pfactor(i,j,k) = (phim(i,j,k) + phin(i,j,k))/(2.*rsm) -1.

                   endif
                enddo
             enddo
          enddo

          call boundc(pfactor)

          c0 = sum( pfactor(1:imax,1:jmax,1:kmax) )*dvol
          call mpi_allreduce(mpi_in_place,c0,1,mpi_real8,mpi_sum,comm_cart,error)

          cp = pot_tot/c0  ! a constant prefactor

          !-- Impose the osmotic pressure jump --!

          px_osp = 0.
          py_osp = 0.
          pz_osp = 0.

          do k=1,kmax
             do j=1,jmax
                do i=1,imax

                   if (abs(slset(i,j,k)) .le. 2.*rsm) then ! jump within this shell

                      jump_x = 0.
                      jump_y = 0.
                      jump_z = 0.
                      jp_x = 0.
                      jp_y = 0.
                      jp_z = 0.

                      !-- x component

                      do nx = -1,1,2

                         fo(-1:1) =   flago(i-1:i+1,j,k)
                         pf(-1:1) = pfactor(i-1:i+1,j,k)

                         if ( fo(0) .ne. fo(nx) ) then

                            if ( fo(0) .eq. 1.) then  ! (i,j,k) is inside overlap
                               c1 = pf(0)
                            else
                               c1 = pf(nx)
                            endif

                            jp_x(nx) = cp*c1
                            if (nx .eq. 1) px_osp(i,j,k) = jp_x(nx)

                         endif
                      enddo

                      jump_x = sum(jp_x)

                      !-- y component

                      do ny = -1,1,2

                         fo(-1:1) =   flago(i,j-1:j+1,k)
                         pf(-1:1) = pfactor(i,j-1:j+1,k)

                         if ( fo(0) .ne. fo(ny) ) then

                            if ( fo(0) .eq. 1.) then  ! (i,j,k) is inside overlap
                               c1 = pf(0)
                            else
                               c1 = pf(ny)
                            endif

                            jp_y(ny) = cp*c1
                            if (ny .eq. 1) py_osp(i,j,k) = jp_y(ny)

                         endif
                      enddo

                      jump_y = sum(jp_y)

                      !-- z component

                      do nz = -1,1,2

                         fo(-1:1) =   flago(i,j,k-1:k+1)
                         pf(-1:1) = pfactor(i,j,k-1:k+1)

                         if ( fo(0) .ne. fo(nz) ) then

                            if ( fo(0) .eq. 1.) then  ! (i,j,k) is inside overlap
                               c1 = pf(0)
                            else
                               c1 = pf(nz)
                            endif

                            jp_z(nz) = cp*c1
                            if (nz .eq. 1) pz_osp(i,j,k) = jp_z(nz)

                         endif
                      enddo

                      jump_z = sum(jp_z)

                      !-- "Laplacian" of osmotic pressure jump (superimposed) --!

                      if (flago(i,j,k) .gt. 0.) then  ! jump from overlap to outside

                         p_jump(i,j,k) = p_jump(i,j,k) + (jump_x/dx**2 + jump_y/dy**2 + jump_z/dz**2)

                      else  ! jump from outside to overlap

                         p_jump(i,j,k) = p_jump(i,j,k) - (jump_x/dx**2 + jump_y/dy**2 + jump_z/dz**2)
                         px_osp(i,j,k) = - px_osp(i,j,k)
                         py_osp(i,j,k) = - py_osp(i,j,k)
                         pz_osp(i,j,k) = - pz_osp(i,j,k)

                      endif

                      !-- Correction to the pressure gradient --!

                      p_x(i,j,k) = p_x(i,j,k) + px_osp(i,j,k)
                      p_y(i,j,k) = p_y(i,j,k) + py_osp(i,j,k)
                      p_z(i,j,k) = p_z(i,j,k) + pz_osp(i,j,k)

                   endif ! jump shell
                enddo ! i
             enddo ! j
          enddo ! k
       enddo  ! n
    enddo ! m

    return
  end subroutine exjump

  !---------------------------------------------------
  !-- The lengthy code below is for the computation of
  !-- a discontinuous viscosity using GFM
  !
  subroutine ls_c2f(dir,phi)

    ! Average level set values from cell centers to faces

    use mod_common
    use mod_bound

    character(len=6), intent(in) :: dir
    real, dimension(0:i1,0:j1,0:k1), intent(out) :: phi

    integer :: i,j,k, cnt

    phi = lvset(0:i1,0:j1,0:k1, 1)  !## now only consider the first level set ##temp

    select case (dir)

    case('u-grid')

       do cnt = 1,nb_cnt
          i = nb_i(cnt)
          j = nb_j(cnt)
          k = nb_k(cnt)
          phi(i,j,k) = ( lvset(i,j,k,1) +lvset(i+1,j,k,1) )/2.  !##temp
       enddo

    case('v-grid')

       do cnt = 1,nb_cnt
          i = nb_i(cnt)
          j = nb_j(cnt)
          k = nb_k(cnt)
          phi(i,j,k) = ( lvset(i,j,k,1) +lvset(i,j+1,k,1) )/2.  !##temp
       enddo

    case('w-grid')

       do cnt = 1,nb_cnt
          i = nb_i(cnt)
          j = nb_j(cnt)
          k = nb_k(cnt)
          phi(i,j,k) = ( lvset(i,j,k,1) +lvset(i,j,k+1,1) )/2.  !##temp
       enddo

    end select

    call boundc(phi) !# wall halos untouched

    return
  end subroutine ls_c2f

  !------------------------
  !
  !
  subroutine gradjump(jump)

    ! Compute the jumpy velocity-gradient matrix

    use mod_common
    use mod_bound
    use mod_interface
    use mod_newtypes

    type (real_3by3_matrix), dimension(-2:i1+2,-2:j1+2,-2:k1+2), intent(out) :: jump !temp bigger than needed

    integer :: i,j,k, cnt, m

    real :: mag

    real, dimension(3) :: dummy, tan_vec1,tan_vec2
    real, dimension(1,3) :: nor_vec, temp
    real, dimension(3,3) :: grad_vel, nor_sq,tan_sq, subt, term1,term2,term3
    real, dimension(0:i1,0:j1,0:k1) :: u,v,w

    u = 0.
    v = 0.
    w = 0.

    subt = 0.

    dvndn = 0.

    ! initialize gradient jump matrix

    forall(i=-2:i1+2,j=-2:j1+2,k=-2:k1+2)
       jump(i,j,k)%grad = subt
    end forall

    if (miu1 .eq. miu2) return

    ! obtain cell-center velocities

    do k = 1,kmax
       do j = 1,jmax
          do i = 1,imax
             u(i,j,k) = (unew(i,j,k) + unew(i-1,j,k))/2.
             v(i,j,k) = (vnew(i,j,k) + vnew(i,j-1,k))/2.
             w(i,j,k) = (wnew(i,j,k) + wnew(i,j,k-1))/2.
          enddo
       enddo
    enddo

    call bounduvw(u,v,w)

    ! obtain normals (also curvature but not used now, can improve efficiency later)

    call get_curvature(lvset(:,:,:,1),curv_cmn)  !## now only the first level set ##temp
    call boundc(curv_cmn)

    ! construct matrices within a narrow band

    do cnt = 1,nb_cnt

       i = nb_i(cnt)
       j = nb_j(cnt)
       k = nb_k(cnt)

       ! velocity gradient matrix

       grad_vel(1,1) = (u(i+1,j,k) - u(i-1,j,k))/2./dx
       grad_vel(1,2) = (u(i,j+1,k) - u(i,j-1,k))/2./dy
       grad_vel(1,3) = (u(i,j,k+1) - u(i,j,k-1))/2./dz
       grad_vel(2,1) = (v(i+1,j,k) - v(i-1,j,k))/2./dx
       grad_vel(2,2) = (v(i,j+1,k) - v(i,j-1,k))/2./dy
       grad_vel(2,3) = (v(i,j,k+1) - v(i,j,k-1))/2./dz
       grad_vel(3,1) = (w(i+1,j,k) - w(i-1,j,k))/2./dx
       grad_vel(3,2) = (w(i,j+1,k) - w(i,j-1,k))/2./dy
       grad_vel(3,3) = (w(i,j,k+1) - w(i,j,k-1))/2./dz

       ! normal and tangential vectors

       nor_vec(1,1) = normal(i,j,k)%x
       nor_vec(1,2) = normal(i,j,k)%y
       nor_vec(1,3) = normal(i,j,k)%z

       dummy(:) = abs(nor_vec(1,:))
       m = minloc(dummy, dim=1)
       dummy(:) = 0.
       dummy(m) = 1.

       tan_vec1(1) =   nor_vec(1,2)*dummy(3) - nor_vec(1,3)*dummy(2)
       tan_vec1(2) = - nor_vec(1,1)*dummy(3) + nor_vec(1,3)*dummy(1)
       tan_vec1(3) =   nor_vec(1,1)*dummy(2) - nor_vec(1,2)*dummy(1)
       mag = dot_product(tan_vec1, tan_vec1)
       tan_vec1 = tan_vec1/sqrt(mag)

       tan_vec2(1) =   nor_vec(1,2)*tan_vec1(3) - nor_vec(1,3)*tan_vec1(2)
       tan_vec2(2) = - nor_vec(1,1)*tan_vec1(3) + nor_vec(1,3)*tan_vec1(1)
       tan_vec2(3) =   nor_vec(1,1)*tan_vec1(2) - nor_vec(1,2)*tan_vec1(1)

       ! corresponding matrices

       nor_sq = matmul(transpose(nor_vec),nor_vec)

       tan_sq(1,:) = 0.
       tan_sq(2,:) = tan_vec1(:)
       tan_sq(3,:) = tan_vec2(:)

       ! first-derivative jump matrix (w/o viscosity term)

       subt = matmul(grad_vel, transpose(tan_sq))
       term1 = matmul(subt, tan_sq)

       subt = matmul(nor_sq, grad_vel)
       term2 = matmul(subt, nor_sq)

       term3 = matmul(transpose(tan_sq), tan_sq)
       subt = matmul(term3, transpose(grad_vel))
       term3 = matmul(subt, nor_sq)

       jump(i,j,k)%grad = term1 + term2 - term3

       ! normal derivative of normal velocity (part of pressure jump)

       temp = matmul(nor_vec,grad_vel)
       dvndn(i,j,k) = dot_product( temp(1,:) ,nor_vec(1,:) )

    enddo

    call boundj(jump)
    call boundc(dvndn)

    return
  end subroutine gradjump

  !---------------------------------------------------------
  !
  !
  subroutine gfm4diff(RDu,RDv,RDw, u,v,w, rho_u,rho_v,rho_w)

    use mod_newtypes

    real, dimension(0:,0:,0:), intent(in) :: u,v,w
    real, dimension(1:imax,1:jmax,1:kmax), intent(in) :: rho_u, rho_v, rho_w
    real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: RDu,RDv,RDw

    integer :: i,j,k, im,ip,jm,jp,km,kp, nb

    real :: s1, theta, miu_pls,miu_mns,miuh
    real :: dudxp,dudxm,dudyp,dudym,dudzp,dudzm
    real :: dvdxp,dvdxm,dvdyp,dvdym,dvdzp,dvdzm
    real :: dwdxp,dwdxm,dwdyp,dwdym,dwdzp,dwdzm
    real :: RD1,RD2,RD3,RD4,RD5
    real :: J_l,J_m,J_r, Ji

    real, dimension(-1:1) :: ph, jmean
    real, dimension(-1:1) :: dudx,dudy,dudz, dvdx,dvdy,dvdz, dwdx,dwdy,dwdz
    real, dimension(0:i1,0:j1,0:k1) :: phi

    type (real_3by3_matrix), dimension(-2:i1+2,-2:j1+2,-2:k1+2) :: jump

    !== Obtain velocity-gradient-jump matrix at cell centers ==!

    call gradjump(jump)  ! visc jump is NOT multiplied

    !-- v-component --! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    call ls_c2f('v-grid',phi)  ! get level set on the v grid
    ph = 0.                    ! level set buffer
    s1 = 3.*dz                 ! tube half width

    dvdx = 0.
    dudy = 0.
    dvdy = 0.
    dvdz = 0.
    dwdy = 0.

    do k = 1,kmax
       do j = 1,jmax
          do i = 1,imax

             ip = i + 1
             jp = j + 1
             kp = k + 1
             im = i - 1
             jm = j - 1
             km = k - 1

             ! continuous first derivatives

             dvdxp = (v(ip,j,k)-v(i ,j,k))/dx
             dvdxm = (v(i ,j,k)-v(im,j,k))/dx

             dudyp = (u(i ,jp,k)-u(i ,j,k))/dy
             dudym = (u(im,jp,k)-u(im,j,k))/dy

             dvdyp = (v(i,jp,k)-v(i,j ,k))/dy
             dvdym = (v(i,j ,k)-v(i,jm,k))/dy

             dvdzp = (v(i,j,kp)-v(i,j,k ))/dz
             dvdzm = (v(i,j,k )-v(i,j,km))/dz

             dwdyp = (w(i,jp,k )-w(i,j,k ))/dy
             dwdym = (w(i,jp,km)-w(i,j,km))/dy

             RD1 = (dvdxp - dvdxm)/dx     ! d/dx(dv/dx)
             RD3 = (dvdyp - dvdym)/dy     ! d/dy(dv/dy)
             RD4 = (dvdzp - dvdzm)/dz     ! d/dz(dv/dz)

             ph(0) = phi(i,j,k)

             if (ph(0) .gt. s1) then       ! safely in fluid 1

                RDv(i,j,k) = miu1*(RD1+RD3+RD4)/rho1

             elseif (ph(0) .lt. -s1) then  ! safely in fluid 2

                RDv(i,j,k) = miu2*(RD1+RD3+RD4)/rho2

             else

                !- 2nd x-derivative -!

                if (TwoD) then

                   RD1 = 0.
                   RD2 = 0.

                else

                   ph(-1) = phi(i-1,j,k)
                   ph( 1) = phi(i+1,j,k)

                   if (minval(ph) .gt. 0.) then      ! in fluid 1

                      RD1 = miu1*RD1
                      RD2 = 0.

                   elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                      RD1 = miu2*RD1
                      RD2 = 0.

                   else

                      if (ph(0) .gt. 0.) then
                         miu_pls = miu1
                         miu_mns = miu2
                      else
                         miu_pls = miu2
                         miu_mns = miu1
                      endif

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i-1,j,k)%grad(2,1) + jump(i-1,j+1,k)%grad(2,1))/2.!#
                      else
                         nb = 1
                         jmean(1) = (jump(i+1,j,k)%grad(2,1) + jump(i+1,j+1,k)%grad(2,1))/2.!#
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(2,1) + jump(i,j+1,k)%grad(2,1))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dvdx(-1) = dvdxm
                      dvdx( 1) = dvdxp

                      dvdx( nb) = miuh*dvdx(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dvdx(-nb) = miu_pls*dvdx(-nb)

                      RD1 = (dvdx(1)-dvdx(-1))/dx  ! d/dx(miu*dv/dx)

                      ! Below L and R locates on the vertices

                      ph(-1) = (phi(i,j,k)+phi(i-1,j,k))/2.
                      ph( 1) = (phi(i,j,k)+phi(i+1,j,k))/2.

                      if (minval(ph) .gt. 0.) then  ! in fluid 1

                         RD2 = miu1*(dudyp-dudym)/dx

                      elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                         RD2 = miu2*(dudyp-dudym)/dx

                      else  ! between L and R

                         if ( ph(-1)*ph(0) .le. 0.) then
                            nb = -1
                            jmean(-1) = (jump(i-1,j,k)%grad(1,2) + jump(i-1,j+1,k)%grad(1,2) &
                                 +jump(i,j  ,k)%grad(1,2) + jump(i,j+1  ,k)%grad(1,2))/4.
                         elseif ( ph(1)*ph(0) .lt. 0.) then
                            nb = 1
                            jmean(1) = (jump(i+1,j,k)%grad(1,2) + jump(i+1,j+1,k)%grad(1,2) &
                                 +jump(i,j  ,k)%grad(1,2) + jump(i,j+1  ,k)%grad(1,2))/4.
                         endif

                         theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                         jmean(0) = (jump(i,j,k)%grad(1,2) + jump(i,j+1,k)%grad(1,2))/2.
                         Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                         Ji = Ji*(miu_pls-miu_mns)

                         miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                         dudy(-1) = dudym
                         dudy( 1) = dudyp

                         dudy( nb) = miuh*dudy(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                         dudy(-nb) = miu_pls*dudy(-nb)

                         RD2 = (dudy(1)-dudy(-1))/dx  ! d/dx(miu*du/dy)

                      endif

                   endif
                endif

                !- 2nd y-derivative -!

                ph(-1) = phi(i,j-1,k)
                ph( 1) = phi(i,j+1,k)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD3 = miu1*RD3

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD3 = miu2*RD3

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j,k)%grad(2,2) + jump(i,j-1,k)%grad(2,2))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j+1,k)%grad(2,2) + jump(i,j+2,k)%grad(2,2))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(2,2) + jump(i,j+1,k)%grad(2,2))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dvdy(-1) = dvdym
                   dvdy( 1) = dvdyp

                   dvdy( nb) = miuh*dvdy(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dvdy(-nb) = miu_pls*dvdy(-nb)

                   RD3 = 2.*(dvdy(1)-dvdy(-1))/dy  ! d/dy(2.*miu*dv/dy)

                endif

                !- 2nd z-derivative -!

                ph(-1) = phi(i,j,k-1)
                ph( 1) = phi(i,j,k+1)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD4 = miu1*RD4
                   RD5 = 0.

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD4 = miu2*RD4
                   RD5 = 0.

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j,k-1)%grad(2,3) + jump(i,j+1,k-1)%grad(2,3))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j,k+1)%grad(2,3) + jump(i,j+1,k+1)%grad(2,3))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(2,3) + jump(i,j+1,k)%grad(2,3))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dvdz(-1) = dvdzm
                   dvdz( 1) = dvdzp

                   dvdz( nb) = miuh*dvdz(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dvdz(-nb) = miu_pls*dvdz(-nb)

                   RD4 = (dvdz(1)-dvdz(-1))/dz  ! d/dz(miu*dv/dz)

                   ! Below L and R locates on the vertices, separated by dz (not 2dz)

                   ph(-1) = (phi(i,j,k)+phi(i,j,k-1))/2.
                   ph( 1) = (phi(i,j,k)+phi(i,j,k+1))/2.

                   if (minval(ph) .gt. 0.) then  ! in fluid 1

                      RD5 = miu1*(dwdyp-dwdym)/dz

                   elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                      RD5 = miu2*(dwdyp-dwdym)/dz

                   else  ! between L and R

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i,j,k-1)%grad(2,3) + jump(i,j+1,k-1)%grad(2,3) &
                              +jump(i,j  ,k)%grad(2,3) + jump(i,j+1  ,k)%grad(2,3))/4.
                      elseif ( ph(1)*ph(0) .lt. 0.) then
                         nb = 1
                         jmean(1) = (jump(i,j,k+1)%grad(2,3) + jump(i,j+1,k+1)%grad(2,3) &
                              +jump(i,j  ,k)%grad(2,3) + jump(i,j+1  ,k)%grad(2,3))/4.
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(2,3) + jump(i,j+1,k)%grad(2,3))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dwdy(-1) = dwdym
                      dwdy( 1) = dwdyp

                      dwdy( nb) = miuh*dwdy(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dwdy(-nb) = miu_pls*dwdy(-nb)

                      RD5 = (dwdy(1)-dwdy(-1))/dz  ! d/dz(miu*dw/dy)

                   endif
                endif

                RDv(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_v(i,j,k)

             endif

          enddo
       enddo
    enddo

    !-- w-component --! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    call ls_c2f('w-grid',phi)  ! get level set on the w grid
    ph = 0.                    ! level set buffer
    s1 = 3.*dz                 ! tube half width

    dwdx = 0.
    dudz = 0.
    dwdy = 0.
    dvdz = 0.
    dwdz = 0.

    do k = 1,kmax
       do j = 1,jmax
          do i = 1,imax

             ip = i + 1
             jp = j + 1
             kp = k + 1
             im = i - 1
             jm = j - 1
             km = k - 1

             ! continuous first derivatives

             dwdxp = (w(ip,j,k)-w(i ,j,k))/dx
             dwdxm = (w(i ,j,k)-w(im,j,k))/dx

             dudzp = (u(i ,j,kp)-u(i ,j,k))/dz
             dudzm = (u(im,j,kp)-u(im,j,k))/dz

             dwdyp = (w(i,jp,k)-w(i,j ,k))/dy
             dwdym = (w(i,j ,k)-w(i,jm,k))/dy

             dvdzp = (v(i,j ,kp)-v(i,j ,k))/dz
             dvdzm = (v(i,jm,kp)-v(i,jm,k))/dz

             dwdzp = (w(i,j,kp)-w(i,j,k ))/dz
             dwdzm = (w(i,j,k )-w(i,j,km))/dz

             RD1 = (dwdxp - dwdxm)/dx  ! d/dx(dw/dx)
             RD3 = (dwdyp - dwdym)/dy  ! d/dy(dw/dy)
             RD5 = (dwdzp - dwdzm)/dz  ! d/dz(dw/dz)

             ph(0) = phi(i,j,k)

             if (ph(0) .gt. s1) then       ! (most likely) safely in fluid 1

                RDw(i,j,k) = miu1*(RD1+RD3+RD5)/rho1

             elseif (ph(0) .lt. -s1) then  ! (most likely) safely in fluid 2

                RDw(i,j,k) = miu2*(RD1+RD3+RD5)/rho2

             else

                !- 2nd x-derivative -!

                if (TwoD) then

                   RD1 = 0.
                   RD2 = 0.

                else

                   ph(-1) = phi(i-1,j,k)
                   ph( 1) = phi(i+1,j,k)

                   if (minval(ph) .gt. 0.) then      ! in fluid 1

                      RD1 = miu1*RD1
                      RD2 = 0.

                   elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                      RD1 = miu2*RD1
                      RD2 = 0.

                   else

                      if (ph(0) .gt. 0.) then
                         miu_pls = miu1
                         miu_mns = miu2
                      else
                         miu_pls = miu2
                         miu_mns = miu1
                      endif

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i-1,j,k)%grad(3,1) + jump(i-1,j,k+1)%grad(3,1))/2.!#
                      else
                         nb = 1
                         jmean(1) = (jump(i+1,j,k)%grad(3,1) + jump(i+1,j,k+1)%grad(3,1))/2.!#
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(3,1) + jump(i,j,k+1)%grad(3,1))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dwdx(-1) = dwdxm
                      dwdx( 1) = dwdxp

                      dwdx( nb) = miuh*dwdx(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dwdx(-nb) = miu_pls*dwdx(-nb)

                      RD1 = (dwdx(1)-dwdx(-1))/dx  ! d/dx(miu*dw/dx)

                      ! Below L and R locates on the vertices

                      ph(-1) = (phi(i,j,k)+phi(i-1,j,k))/2.
                      ph( 1) = (phi(i,j,k)+phi(i+1,j,k))/2.

                      if (minval(ph) .gt. 0.) then  ! in fluid 1

                         RD2 = miu1*(dudzp-dudzm)/dx

                      elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                         RD2 = miu2*(dudzp-dudzm)/dx

                      else  ! between L and R

                         if ( ph(-1)*ph(0) .le. 0.) then
                            nb = -1
                            jmean(-1) = (jump(i-1,j,k)%grad(1,3) + jump(i-1,j,k+1)%grad(1,3) &
                                 +jump(i,j  ,k)%grad(1,3) + jump(i,j  ,k+1)%grad(1,3))/4.
                         elseif ( ph(1)*ph(0) .lt. 0.) then
                            nb = 1
                            jmean(1) = (jump(i+1,j,k)%grad(1,3) + jump(i+1,j,k+1)%grad(1,3) &
                                 +jump(i,j  ,k)%grad(1,3) + jump(i,j  ,k+1)%grad(1,3))/4.
                         endif

                         theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                         jmean(0) = (jump(i,j,k)%grad(1,3) + jump(i,j,k+1)%grad(1,3))/2.
                         Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                         Ji = Ji*(miu_pls-miu_mns)

                         miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                         dudz(-1) = dudzm
                         dudz( 1) = dudzp

                         dudz( nb) = miuh*dudz(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                         dudz(-nb) = miu_pls*dudz(-nb)

                         RD2 = (dudz(1)-dudz(-1))/dx  ! d/dx(miu*du/dz)

                      endif

                   endif
                endif

                !- 2nd y-derivative -!

                ph(-1) = phi(i,j-1,k)
                ph( 1) = phi(i,j+1,k)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD3 = miu1*RD3
                   RD4 = 0.

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD3 = miu2*RD3
                   RD4 = 0.

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j-1,k)%grad(3,2) + jump(i,j-1,k+1)%grad(3,2))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j+1,k)%grad(3,2) + jump(i,j+1,k+1)%grad(3,2))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(3,2) + jump(i,j,k+1)%grad(3,2))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dwdy(-1) = dwdym
                   dwdy( 1) = dwdyp

                   dwdy( nb) = miuh*dwdy(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dwdy(-nb) = miu_pls*dwdy(-nb)

                   RD3 = (dwdy(1)-dwdy(-1))/dy  ! d/dy(miu*dw/dy)

                   ! Below L and R locates on the vertices

                   ph(-1) = (phi(i,j,k)+phi(i,j-1,k))/2.
                   ph( 1) = (phi(i,j,k)+phi(i,j+1,k))/2.

                   if (minval(ph) .gt. 0.) then  ! in fluid 1

                      RD4 = miu1*(dvdzp-dvdzm)/dy

                   elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                      RD4 = miu2*(dvdzp-dvdzm)/dy

                   else  ! between L and R

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i,j-1,k)%grad(2,3) + jump(i,j-1,k+1)%grad(2,3) &
                              +jump(i,j  ,k)%grad(2,3) + jump(i,j  ,k+1)%grad(2,3))/4.
                      elseif ( ph(1)*ph(0) .lt. 0.) then
                         nb = 1
                         jmean(1) = (jump(i,j+1,k)%grad(2,3) + jump(i,j+1,k+1)%grad(2,3) &
                              +jump(i,j  ,k)%grad(2,3) + jump(i,j  ,k+1)%grad(2,3))/4.
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(2,3) + jump(i,j,k+1)%grad(2,3))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dvdz(-1) = dvdzm
                      dvdz( 1) = dvdzp

                      dvdz( nb) = miuh*dvdz(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dvdz(-nb) = miu_pls*dvdz(-nb)

                      RD4 = (dvdz(1)-dvdz(-1))/dy  ! d/dy(miu*dv/dz)

                   endif

                endif

                !- 2nd z-derivative -!

                ph(-1) = phi(i,j,k-1)
                ph( 1) = phi(i,j,k+1)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD5 = miu1*RD5

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD5 = miu2*RD5

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j,k)%grad(3,3) + jump(i,j,k-1)%grad(3,3))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j,k+1)%grad(3,3) + jump(i,j,k+2)%grad(3,3))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(3,3) + jump(i,j,k+1)%grad(3,3))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dwdz(-1) = dwdzm
                   dwdz( 1) = dwdzp

                   dwdz( nb) = miuh*dwdz(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dwdz(-nb) = miu_pls*dwdz(-nb)

                   RD5 = 2.*(dwdz(1)-dwdz(-1))/dz  ! d/dz(2.*miu*dw/dz)

                endif

                RDw(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_w(i,j,k)

             endif

          enddo
       enddo
    enddo

    !-- u-component --! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    if (TwoD) then
       RDu = 0.
       return
    endif

    call ls_c2f('u-grid',phi)  ! get level set on the u grid
    ph = 0.                    ! level set buffer
    s1 = 3.*dz                 ! tube half width

    dudx = 0.
    dudy = 0.
    dvdx = 0.
    dudz = 0.
    dwdx = 0.

    do k = 1,kmax
       do j = 1,jmax
          do i = 1,imax

             ip = i + 1
             jp = j + 1
             kp = k + 1
             im = i - 1
             jm = j - 1
             km = k - 1

             ! continuous first derivatives

             dudxp = (u(ip,j,k)-u(i ,j,k))/dx
             dudxm = (u(i ,j,k)-u(im,j,k))/dx

             dudyp = (u(i,jp,k)-u(i,j ,k))/dy
             dudym = (u(i,j ,k)-u(i,jm,k))/dy

             dvdxp = (v(ip,j ,k)-v(i,j ,k))/dx
             dvdxm = (v(ip,jm,k)-v(i,jm,k))/dx

             dudzp = (u(i,j,kp)-u(i,j,k ))/dz
             dudzm = (u(i,j,k )-u(i,j,km))/dz

             dwdxp = (w(ip,j,k )-w(i,j,k ))/dx
             dwdxm = (w(ip,j,km)-w(i,j,km))/dx

             RD1 = (dudxp -dudxm)/dx  ! d/dx(du/dx)
             RD2 = (dudyp -dudym)/dy  ! d/dy(du/dy)
             RD4 = (dudzp -dudzm)/dz  ! d/dz(du/dz)

             ph(0) = phi(i,j,k)

             if (ph(0) .gt. s1) then       ! (most likely) safely in fluid 1

                RDu(i,j,k) = miu1*(RD1+RD2+RD4)/rho1

             elseif (ph(0) .lt. -s1) then  ! (most likely) safely in fluid 2

                RDu(i,j,k) = miu2*(RD1+RD2+RD4)/rho2

             else

                !- 2nd x-derivative -!

                ph(-1) = phi(i-1,j,k)
                ph( 1) = phi(i+1,j,k)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD1 = miu1*RD1

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD1 = miu2*RD1

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j,k)%grad(1,1) + jump(i-1,j,k)%grad(1,1))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i+1,j,k)%grad(1,1) + jump(i+2,j,k)%grad(1,1))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(1,1) + jump(i+1,j,k)%grad(1,1))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dudx(-1) = dudxm
                   dudx( 1) = dudxp

                   dudx( nb) = miuh*dudx(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dudx(-nb) = miu_pls*dudx(-nb)

                   RD1 = 2.*(dudx(1)-dudx(-1))/dx  ! d/dx(2.*miu*du/dx)

                endif

                !- 2nd y-derivative -!

                ph(-1) = phi(i,j-1,k)
                ph( 1) = phi(i,j+1,k)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD2 = miu1*RD2
                   RD3 = 0.

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD2 = miu2*RD2
                   RD3 = 0.

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j-1,k)%grad(1,2) + jump(i+1,j-1,k)%grad(1,2))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j+1,k)%grad(1,2) + jump(i+1,j+1,k)%grad(1,2))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(1,2) + jump(i+1,j,k)%grad(1,2))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dudy(-1) = dudym
                   dudy( 1) = dudyp

                   dudy( nb) = miuh*dudy(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dudy(-nb) = miu_pls*dudy(-nb)

                   RD2 = (dudy(1)-dudy(-1))/dy  ! d/dy(miu*du/dy)

                   ! Below L and R locates on the vertices

                   ph(-1) = (phi(i,j,k)+phi(i,j-1,k))/2.
                   ph( 1) = (phi(i,j,k)+phi(i,j+1,k))/2.

                   if (minval(ph) .gt. 0.) then  ! in fluid 1

                      RD3 = miu1*(dvdxp-dvdxm)/dy

                   elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                      RD3 = miu2*(dvdxp-dvdxm)/dy

                   else  ! between L and R

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i,j-1,k)%grad(2,1) + jump(i+1,j-1,k)%grad(2,1) &
                              +jump(i,j  ,k)%grad(2,1) + jump(i+1,j  ,k)%grad(2,1))/4.
                      elseif ( ph(1)*ph(0) .lt. 0.) then
                         nb = 1
                         jmean(1) = (jump(i,j+1,k)%grad(2,1) + jump(i+1,j+1,k)%grad(2,1) &
                              +jump(i,j  ,k)%grad(2,1) + jump(i+1,j  ,k)%grad(2,1))/4.
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(2,1) + jump(i+1,j,k)%grad(2,1))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dvdx(-1) = dvdxm
                      dvdx( 1) = dvdxp

                      dvdx( nb) = miuh*dvdx(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dvdx(-nb) = miu_pls*dvdx(-nb)

                      RD3 = (dvdx(1)-dvdx(-1))/dy  ! d/dy(miu*dv/dx)

                   endif
                endif

                !- 2nd z-derivative -!

                ph(-1) = phi(i,j,k-1)
                ph( 1) = phi(i,j,k+1)

                if (minval(ph) .gt. 0.) then      ! in fluid 1

                   RD4 = miu1*RD4
                   RD5 = 0.

                elseif (maxval(ph) .lt. 0.) then  ! in fluid 2

                   RD4 = miu2*RD4
                   RD5 = 0.

                else

                   if (ph(0) .gt. 0.) then
                      miu_pls = miu1
                      miu_mns = miu2
                   else
                      miu_pls = miu2
                      miu_mns = miu1
                   endif

                   if ( ph(-1)*ph(0) .le. 0.) then
                      nb = -1
                      jmean(-1) = (jump(i,j,k-1)%grad(1,3) + jump(i+1,j,k-1)%grad(1,3))/2.
                   else
                      nb = 1
                      jmean(1) = (jump(i,j,k+1)%grad(1,3) + jump(i+1,j,k+1)%grad(1,3))/2.
                   endif

                   theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                   jmean(0) = (jump(i,j,k)%grad(1,3) + jump(i+1,j,k)%grad(1,3))/2.
                   Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                   Ji = Ji*(miu_pls-miu_mns)

                   miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                   dudz(-1) = dudzm
                   dudz( 1) = dudzp

                   dudz( nb) = miuh*dudz(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                   dudz(-nb) = miu_pls*dudz(-nb)

                   RD4 = (dudz(1)-dudz(-1))/dz  ! d/dz(miu*du/dz)

                   ! Below L and R locates on the vertices

                   ph(-1) = (phi(i,j,k)+phi(i,j,k-1))/2.
                   ph( 1) = (phi(i,j,k)+phi(i,j,k+1))/2.

                   if (minval(ph) .gt. 0.) then  ! in fluid 1

                      RD5 = miu1*(dwdxp-dwdxm)/dz

                   elseif ( maxval(ph) .lt. 0. ) then  ! in fluid 2

                      RD5 = miu2*(dwdxp-dwdxm)/dz

                   else  ! between L and R

                      if ( ph(-1)*ph(0) .le. 0.) then
                         nb = -1
                         jmean(-1) = (jump(i,j,k-1)%grad(3,1) + jump(i+1,j,k-1)%grad(3,1) &
                              +jump(i,j  ,k)%grad(3,1) + jump(i+1,j  ,k)%grad(3,1))/4.
                      elseif ( ph(1)*ph(0) .lt. 0.) then
                         nb = 1
                         jmean(1) = (jump(i,j,k+1)%grad(3,1) + jump(i+1,j,k+1)%grad(3,1) &
                              +jump(i,j  ,k)%grad(3,1) + jump(i+1,j  ,k)%grad(3,1))/4.
                      endif

                      theta = abs(ph(nb))/( abs(ph(nb)) +abs(ph(0)) )  ! linear interpln.

                      jmean(0) = (jump(i,j,k)%grad(3,1) + jump(i+1,j,k)%grad(3,1))/2.
                      Ji = theta*jmean(0) +(1.-theta)*jmean(nb)  ! interpolated jump at the interface
                      Ji = Ji*(miu_pls-miu_mns)

                      miuh = miu_pls*miu_mns/(miu_pls*theta+miu_mns*(1.-theta))  ! effective viscosity

                      dwdx(-1) = dwdxm
                      dwdx( 1) = dwdxp

                      dwdx( nb) = miuh*dwdx(nb) + miuh/miu_mns*theta*Ji  ! interface btw 0 and nb
                      dwdx(-nb) = miu_pls*dwdx(-nb)

                      RD5 = (dwdx(1)-dwdx(-1))/dz  ! d/dz(miu*dw/dx)

                   endif
                endif

                RDu(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_u(i,j,k)

             endif

          enddo
       enddo
    enddo

    return
  end subroutine gfm4diff

end module mod_gfm

max_stars_repo_name: GeZhouyang/ICLS-release

max_issues_repo_name: GeZhouyang/ICLS-release

max_forks_repo_name: GeZhouyang/ICLS-release


text:
program test
  character(4)             :: char_length
  integer                  :: line_length

  line_length = 100

  write(char_length, '(I4)') line_length

  print *, char_length

end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
program testcap
    use defcapsol
    implicit none
    integer :: n, m, l_js, verbose, d_steps, test
    real(dp) :: h0, rho_max, Z_max, rtip, theta_deg, hcone, rcant, hsam, dcant, eps_r
    real(dp) :: C_unitless   ! finest grid spacing and simulation box size (truncation lengths)

    character*10 method_in

    real(dp), allocatable :: hess(:, :) ! grid-spacing and coordinates of each grid point

    ! Nuni = 5
    n=500
    m=500
    l_js=20
    h0=0.5_dp
    rho_max=1d6
    Z_max=1d6
    d_steps=20
    rtip=20.0
    theta_deg=15.0
    hcone=15000.0
    rcant=15000.0
    dcant=500.0
    eps_r=5.0
    hsam=10.0
    method_in='LAPACK'
    test=0
    verbose=2
    C_unitless = 0.0_dp

    allocate (Hess(0:min(m+ l_js + d_steps + 1, n+1), (n+1)*( l_js + d_steps  + m + 1)))

    CALL CapSolCyl(n,m,l_js, h0,rho_max,z_max, d_steps, rtip,theta_deg, hcone, rcant, &
    dcant, eps_r,hsam, method_in, test, verbose, C_unitless, Hess)

    deallocate (Hess)
end program testcap

max_stars_repo_name: maroniea/xsede-spm

max_issues_repo_name: maroniea/xsede-spm

max_forks_repo_name: maroniea/xsede-spm