! -*- f90 -*-
! file leastsq.f90
! Main Geodesic-Bold-BroydenUpdate-Levenberg-Marquardt routine
! version 1.1

SUBROUTINE geodesiclm(func, jacobian, Avv, &
     & x, fvec, fjac, n, m, &
     & callback, info, &
     & analytic_jac, analytic_Avv, &
     & center_diff, h1, h2,&
     & dtd, damp_mode, &
     & niters, nfev, njev, naev, &
     & maxiter, maxfev, maxjev, maxaev, maxlam, minlam, &
     & artol, Cgoal, gtol, xtol, xrtol, ftol, frtol, &
     & converged, &
     & print_level, print_unit, &
     & imethod, iaccel, ibold, ibroyden, &
     & initialfactor, factoraccept, factorreject, avmax)

!*****************************************************************
!
!    subroutine geodesicLM
!
!    The purpose of geolevmar is to minimize the sum of the squares
!    of m nonlinear functions of n variables by a modification of
!    the Levenberg-Marquardt algorithm that utilizes the geodesic
!    acceleration step correction, bold acceptance criterion, and
!    a Broyden update of the jacobian matrix.  The method employs one
!    of several possible schemes for updating the Levenberg-Marquardt
!    parameter.  The user must provide a subroutine which calcualtes
!    the functions, and optionally the jacobian and a directional
!    derivative of the functions.  The latter two will be estimated
!    by finite differences if not supplied.
!
!    If you use this code, please acknowledge such by referencing one
!    one of the following papers in any published work:
!
!    Transtrum M.K., Machta B.B., and Sethna J.P, Why are nonlinear
!    fits to data so challenging?  Phys. Rev. Lett. 104, 060201 (2010)
!
!    Transtrum M.K., Machta B.B., and Sethna J.P., The geometry of
!    nonlinear least squares with applications to sloppy model and
!    optimization.  Phys. Rev. E. 80, 036701 (2011)
!
!
!    The subroutine statement is:
!
!    geodesicLM(func, jacobian, Avv, x, fvec, fjac, n, m, callback, info,
!              analytic_jac, analytic_Avv, center_diff, h1, h2,
!              dtd, damp_mode, niteres, nfev, njev, naev,
!              maxiters, maxfev, maxjev, maxaev, maxlam, minlam,
!              artol, Cgoal, gtol, xtol, xrtol, ftol, frtol,
!              converged, print_level, print_unit,
!              imethod, iaccel, ibold, ibroyden,
!              initialfactor, factoraccept, factorreject, avmax)
!
!    where
!
!    func is a user supplied subroutine which calculates the functions and
!    should be written as follows:
!
!      subroutine func(m, n, x, fvec)
!      integer m, n
!      double precision x(n), fvec(m)
!      --------------------------------------------------------------------
!      calculates the function at x and returns their values in fvec
!      x, m, and n should be left unchanged
!      --------------------------------------------------------------------
!      end subroutine func
!
!    jacobian is a user supplied subroutine which calculates the jacobian of
!    of the functions if analytic_jac is .TRUE.
!    jacobian should be writen as follows
!
!      subroutine jacobian(m, n, x, fjac)
!      integer m, n
!      double precision x(n), fjac(m,n)
!      --------------------------------------------------------------------
!      calculates the jacobian at x and returns their values in fjac
!      x, m, and n should be left unchanged
!      --------------------------------------------------------------------
!      end subroutine jacobian
!
!    Avv is a user supplied subroutine which calculates the directional
!    second derivative of the functions if analytic_Avv is .TRUE.
!    Avv should be writen as follows
!
!      subroutine Avv(m, n, x, v, acc)
!      integer m, n
!      double precision x(n), v(n), acc(m)
!      --------------------------------------------------------------------
!      calculates the directional second derivative at x in the direction
!      of v and returns the values in acc
!      x, v, m, and n should be left unchanged
!      --------------------------------------------------------------------
!      end subroutine Avv
!
!    x is an array of length n.  On input it contains an initial estimate of
!    the solution.  On exit, it contains the final estimate of the solution.
!
!    fvec is an output array of length m containing the funtion evaluation at
!    the final solution
!
!    fjac is an output array of dimension(m,n) containing the jacobian evaluation
!    the final solution.  The array MATMUL( TRANSPOSE(fjac), fjac) is an estimate
!    of the covariance matrix of parameters at the final solution
!
!    n an input integer set to the number of parameters
!
!    m an input integer set to the number of functions
!
!    callback a user supplied subroutine to be called after each iteration of the
!    algorithm.
!    callback should be written as follows
!
!      subroutine callback(m,n,x,v,a,fvec,fjac,acc,lam,dtd,fvec_new,accepted,info)
!      integer m, n, accepted, info
!      double precision x(n), v(n), a(n), fvec(m), fjac(m,n), acc(m), lam, dtd(n,n), fvec_new(m)
!      --------------------------------------------------------------------
!      m, n, x, v, a, fvec, fjac, acc, lam, dtd, fvec_new, accepted, should be left unchanged
!      On input, info = 0 and should be changed to a nonzero value if the user
!      wishes to terminate calculation
!      --------------------------------------------------------------------
!      end subroutine callback
!
!    info an output integer set to a nonzero value if the user terminated the routine
!    (see callback).
!
!    analytic_jac an input boolean set to .TRUE. if the subroutine jacobian calculates
!    the jacobian.  If .FALSE. then a finite difference estimate will be used.
!
!    analytic_Avv an input boolean set to .TRUE. if the subroutine Avv calculates
!    the directional second derivative.  If .FALSE. then a finite difference estimate
!    will be used.
!
!    center_diff an input boolean.  If finite differences are used to estimate the jacobian
!    then center differences will used if center_diff is .TRUE., otherwise, forward
!    differences will be used.  Note that center differences are more accurate by require
!    more function evaluations.
!
!    h1 an input double precision specifying the step size for the finite difference estimates
!    of the jacobian.
!
!    h2 an input double precision specifying the steps ize for the finite difference estiamtes
!    of the directional second derivative.
!
!    dtd a double precision array of dimension(n,n).  dtd is used as the damping matrix in the
!    Levenberg-Marquardt routine.  It's exact treatment is specified by the damp_mode input.
!
!    damp_mode an input integer specifying the details of the LM damping as follows:
!      damp_mode = 0: dtd is set to the identity.
!      damp_mode = 1: dtd should be a positive definite, diagonal matrix whose entries are dynamically
!                updated based on the elements of the jacobian.
!
!    niters an output integer specifying the number of iterations of the algorithm.
!
!    nfev an output integer specifying the number of calls to func.
!
!    njev an output integer specifying the number of calls to jacobian.
!
!    naev an output integer specifying the number of calls to Avv.
!
!    maxiter an input integer specifying the maximum number of routine iterations.
!
!    maxfev an input integer specifying the maximum number of function calls
!    if maxfev = 0, then there is no limit to the number of function calls.
!
!    maxjev an input integer specifying the maximum number of jacobian calls
!    if maxjev = 0, then there is no limit to the number of jacobian calls.
!
!    maxaev an input integer specifying the maximum number of Avv calls
!    if maxaev = 0, then there is no limit to the number of Avv calls.
!
!    maxlam an input double precision specifying the maximum allowed value of
!    the damping term lambda. If this is negative, then there is no limit.
!
!    minlam an input double precision specifying the minimum allowed value of
!    the damping term lambda. If lambda is smaller than this value for three consecutive steps
!    the routine terminates.  If this is negative, then there is no limit.
!
!    artol an input double precision.  The method will terminate when the cosine of the
!    angle between the residual vector and the range of the jacobian is less than artol.
!
!    Cgoal an input double precision.  The method will terminate when the cost (one half
!    the sum of squares of the function) falls below Cgoal.
!
!    gtol an input double precision.  The method will terminate when norm of Cost gradient
!    falls below gtol.
!
!    xtol an input double precision.  The method will terminate when parameters change by
!    less than xtol.
!
!    xrtol an input double precision.  The method will terminate if the relative change in
!    each of the parameters is less than xrtol.
!
!    ftol an input double precision.  The method will termiante if the Cost fails to decrease
!    by more than ftol for 3 consecutive iterations.
!
!    frtol an input double precision.  The method will terminate if the relative decrease in
!    Cost is less than frtol 3 consecutive iterations.
!
!    converged an output integer indicated the reason for termination:
!      converged = 1: artol
!      converged = 2: Cgoal
!      converged = 3: gtol
!      converged = 4: xtol
!      converged = 5: xrtol
!      converged = 6: ftol
!      converged = 7: frtol
!      converged = -1: maxiters exeeded
!      converged = -2: maxfev exceeded
!      converged = -3: maxjev exceeded
!      converged = -4: maxaev exceeded
!      converged = -10: user requested termination in callback via info
!      converged = -11: Either the initial function evalaution or subsequent jacobian
!                       evaluations produced Nans.
!
!    print_level an input integer specifying the amount of details to be printed.
!    acceptable values range from 0 to 5, with larger number printing more details.
!
!    print_unit an input integer specifying the unit number details should be written to.
!
!    imethod an input integer specifying the method for updating the LM parameter
!      imethod = 0: adjusted by fixed factors after accepted/rejected steps
!      imethod = 1: adjusted as described in Nielson
!      imethod = 2: adjusted according to an unpublished method due to Cyrus Umrigar and Peter Nightingal
!      imethod = 10: step size Delta adjusted by fixed factors after accepted/rejected steps
!      imethod = 11: step size adjusted as described in More'
!
!    initialfactor an input double precision for specifying either the initial LM parameter
!    of the initial step size.
!
!    factoraccept an input double precision (larger than 1.0) specifying the factor by which
!    either the LM parameter or the step size will be adjusted after an accepted step if
!    imethod = 0 or 10
!
!    factorreject an input double precision (larger than 1.0) specifying the factor by which
!    either the LM parameter of the step size will be adjusted after a rejected step if
!    imethod = 0 or 10
!
!    avmax an input double precision specifying the maximum norm of the geodesic acceleration
!    relative to the velocity vector.
!
!*****************************************************************

  IMPLICIT NONE
  !! Passed parameters
  EXTERNAL func, jacobian, Avv, callback

  REAL (KIND=8) x(n), fvec(m), fjac(m,n)
  INTEGER n, m
  LOGICAL analytic_jac, analytic_Avv, center_diff
  REAL (KIND=8) h1, h2
  REAL (KIND=8) dtd(n,n)
  INTEGER damp_mode, info
  INTEGER niters, nfev, njev, naev
  INTEGER maxiter, maxfev, maxaev, maxjev
  REAL (KIND=8) maxlam, minlam, artol, Cgoal, gtol, xtol, xrtol, ftol, frtol
  INTEGER converged
  INTEGER print_level, print_unit
  INTEGER iaccel, ibroyden, ibold, imethod
  REAL (KIND=8) avmax, initialfactor, factoraccept, factorreject

  !! Internal parameters

  REAL (KIND=8) acc(m), v(n), vold(n), a(n), lam, delta, cos_alpha, av
  REAL (KIND=8) fvec_new(m), fvec_best(m), C, Cnew, Cbest, Cold
  REAL (KIND=8) x_new(n), x_best(n)
  REAL (KIND=8) jtj(n,n), g(n,n)
  REAL (KIND=8) temp1, temp2, pred_red, dirder, actred, rho, a_param
  INTEGER i, j, istep, accepted, counter

  character(16) :: converged_info(-11:7)
  LOGICAL jac_uptodate, jac_force_update, valid_result

  ! strings for concluding print statement
  converged_info = '????????'
  converged_info(1) = 'artol reached'
  converged_info(2) = 'Cgoal reached'
  converged_info(3) = 'gtol reached'
  converged_info(4) = 'xtol reached'
  converged_info(5) = 'xrtol reached'
  converged_info(6) = 'ftol reached'
  converged_info(7) = 'frtol reached'
  converged_info(-1) = 'maxiters exeeded'
  converged_info(-2) = 'maxfev exceeded'
  converged_info(-3) = 'maxjev exceeded'
  converged_info(-4) = 'maxaev exceeded'
  converged_info(-10) = 'User Termination '
  converged_info(-11) = 'NaN Produced'

  IF(print_level .GE. 1) THEN
     WRITE(print_unit, *) "Optimizing with Geodesic-Levenberg-Marquardt algorithm, version 1.1"
     WRITE(print_unit, *) "Method Details:"
     WRITE(print_unit, *) "  Update method:   ", imethod
     WRITE(print_unit, *) "  acceleration:    ", iaccel
     WRITE(print_unit, *) "  Bold method:     ", ibold
     WRITE(print_unit, *) "  Broyden updates: ", ibroyden
     FLUSH(print_unit)
  ENDIF

  !! Initialize variables
  niters = 0
  nfev = 0
  naev = 0
  njev = 0
  converged = 0
  v(:) = 0.0d+0
  vold(:) = 0.0d+0
  a(:) = 0.0d+0
  cos_alpha = 1.0d+0
  av = 0.0d+0
  a_param = 0.5

  accepted = 0
  counter = 0
  CALL func(m,n,x,fvec)
  nfev = nfev + 1

  C = 0.5d+0*DOT_PRODUCT(fvec,fvec)
  IF(print_level .GE. 1) THEN
     WRITE(print_unit, *) "  Initial Cost:    ", C
     FLUSH(print_unit)
  ENDIF
  valid_result = .TRUE.
  !! Check for nans in fvec
  checkfvec: DO i = 1,m
     IF(fvec(i) /= fvec(i)) THEN
        valid_result = .FALSE.
        EXIT checkfvec
     END IF
  END DO checkfvec
  IF (.NOT. valid_result) THEN
     converged = -11
     maxiter = 0
  ENDIF
  Cbest = C
  fvec_best = fvec
  x_best = x
  IF(analytic_jac) THEN
     CALL jacobian(m,n,x,fjac)
     njev = njev + 1
  ELSE
     CALL fdjac(m,n,x,fvec,fjac,func,h1,center_diff)
     IF (center_diff) THEN
        nfev = nfev + 2*n
     ELSE
        nfev = nfev + n
     ENDIF
  ENDIF
  jac_uptodate = .TRUE.
  jac_force_update = .FALSE.
  jtj = MATMUL(TRANSPOSE(fjac), fjac)

  !! Check fjac for nans
  valid_result = .TRUE.
  checkfjac_initial: DO i = 1,m
     DO j = 1,n
        IF(fjac(i,j) /= fjac(i,j)) THEN
           valid_result = .FALSE.
           EXIT checkfjac_initial
        END IF
     END DO
  END DO checkfjac_initial
  IF( .NOT. valid_result) THEN
     converged = -11
     maxiter = 0
  ENDIF

  acc(:) = 0.0d+0
  a(:) = 0.0d+0

  !! Initialize scaling matrix
  IF(damp_mode.EQ.0) THEN
     dtd(:,:) = 0.0d+0
     DO i=1,n
        dtd(i,i) = 1.0d+0
     END DO
  ELSEIF(damp_mode.EQ.1) THEN
     DO i = 1,n
        dtd(i,i) = MAX(jtj(i,i),dtd(i,i))
     END DO
  ENDIF

  !! Initialize lambda
  IF(imethod .LT. 10) THEN
     lam = jtj(1,1)
     DO i = 2,n
        lam = MAX(jtj(i,i),lam)
     END DO
     lam = lam * initialfactor
  !! Initialize step bound if using trust region method
  ELSEIF(imethod .GE. 10) THEN
     delta = initialfactor*SQRT(DOT_PRODUCT(x,MATMUL(dtd,x)))
     lam = 1.0d+0
     IF(delta .EQ. 0.0d+0) delta = 100d+0
     IF( converged .EQ. 0) CALL TrustRegion(n,m,fvec,fjac,dtd,delta,lam) !! Do not call this if there were nans in either fvec or fjac
  ENDIF

  !! Main Loop
  main: DO istep=1, maxiter

     info = 0
     CALL callback(m,n,x,v,a,fvec,fjac,acc,lam,dtd,fvec_new,accepted,info)
     IF( info .NE. 0) THEN
        converged = -10
        exit main
     ENDIF
     !! Update Functions

     !! Full or partial Jacobian Update?
     IF (accepted .GT. 0 .AND. ibroyden .LE. 0) jac_force_update = .TRUE.
     IF (accepted + ibroyden .LE. 0 .AND. .NOT. jac_uptodate) jac_force_update = .TRUE.  !Force jac update after too many failed attempts

     IF (accepted .GT. 0 .AND. ibroyden .GT. 0 .AND. .NOT. jac_force_update) THEN !! Rank deficient update of jacobian matrix
        CALL UPDATEJAC(m,n,fjac, fvec, fvec_new, acc, v, a)
        jac_uptodate = .FALSE.
     ENDIF

     IF( accepted .GT. 0) THEN !! Accepted step
        fvec = fvec_new
        x = x_new
        vold = v
        C = Cnew
        IF( C .LE. Cbest) THEN
           x_best = x
           Cbest = C
           fvec_best = fvec
        ENDIF
     ENDIF

     IF( jac_force_update ) THEN !! Full rank update of jacobian
        IF(analytic_jac) THEN
           CALL jacobian(m,n,x,fjac)
           njev = njev + 1
        ELSE
           CALL fdjac(m,n,x,fvec,fjac,func,h1,center_diff)
           IF (center_diff) THEN
              nfev = nfev + 2*n
           ELSE
              nfev = nfev + n
           ENDIF
        ENDIF
        jac_uptodate = .TRUE.
        jac_force_update = .FALSE.
     ENDIF

     !! Check fjac for nans
     valid_result = .TRUE.
     checkfjac: DO i = 1,m
        DO j = 1,n
           IF(fjac(i,j) /= fjac(i,j)) THEN
              valid_result = .FALSE.
              EXIT checkfjac
           END IF
        END DO
     END DO checkfjac

     IF (valid_result) THEN !! If no nans in jacobian

        jtj = MATMUL(TRANSPOSE(fjac), fjac)

        !! Update Scaling/lam/TrustRegion
        IF(istep .GT. 1) THEN !! Only necessary after first step
           IF( damp_mode .EQ. 1) THEN
              DO i = 1,n
                 dtd(i,i) = MAX(jtj(i,i),dtd(i,i))
              END DO
           ENDIF
           !! Can add other lam-delta update methods
           SELECT CASE( imethod )
           CASE( 0 )
              !! Update lam directly by fixed factors
              CALL Updatelam_factor(lam, accepted, factoraccept, factorreject)
           CASE( 1 )
              !! Update lam directly based on Gain Factor rho (see Nielson reference)
              CALL Updatelam_nelson(lam, accepted, factoraccept, factorreject, rho)
           CASE( 2 )
              !! Update lam directly using method of Umrigar and Nightingale [unpublished]
              CALL Updatelam_Umrigar(m,n,lam,accepted, v, vold, fvec, fjac, dtd, a_param, Cold, Cnew)
           CASE(10)
              !! Update delta by fixed factors
              CALL UpdateDelta_factor(delta, accepted, factoraccept, factorreject)
              CALL TrustRegion(n,m,fvec, fjac, dtd, delta, lam)
           CASE(11)
              !! Update delta as described in More' reference
              CALL UpdateDelta_more(delta, lam, n, x, dtd, rho, C, Cnew, dirder, actred, av, avmax)
              CALL TrustRegion(n,m,fvec, fjac, dtd, delta, lam)
           END SELECT
        ENDIF

        !! Propose Step
        !! metric aray
        g = jtj + lam*dtd
        !! Cholesky decomposition
        CALL DPOTRF('U', n, g, n, info)
        !! CALL inv(n, g, info)
     ELSE !! If nans in jacobian
        converged = -11
        exit main
     ENDIF

     IF(info .EQ. 0) THEN  !! If matrix decomposition successful:
        !! v = -1.0d+0*MATMUL(g,MATMUL(fvec,fjac)) ! velocity
        v = -1.0d+0*MATMUL(fvec, fjac)
        CALL DPOTRS('U', n, 1, g, n, v, n, info)

        ! Calcualte the predicted reduction and the directional derivative -- useful for updating lam methods
        temp1 = 0.5d+0*DOT_PRODUCT(v,MATMUL(jtj, v))/C
        temp2 = 0.5d+0*lam*DOT_PRODUCT(v,MATMUL(dtd,v))/C
        pred_red = temp1 + 2.0d+0*temp2
        dirder = -1.0d+0*(temp1 + temp2)
        ! calculate cos_alpha -- cos of angle between step direction (in data space) and residual vector
        cos_alpha = ABS(DOT_PRODUCT(fvec, MATMUL(fjac, v)))
        cos_alpha = cos_alpha/SQRT( DOT_PRODUCT(fvec, fvec)*DOT_PRODUCT(MATMUL(fjac,v), MATMUL(fjac, v)) )
        IF ( imethod .LT. 10) delta = SQRT(DOT_PRODUCT(v, MATMUL(dtd, v)))  !! Update delta if not set directly
        !! update acceleration
        IF(iaccel .GT. 0) THEN
           IF( analytic_Avv ) THEN
              CALL Avv(m,n,x,v,acc)
              naev = naev + 1
           ELSE
              CALL FDAvv(m,n,x,v,fvec, fjac, func, acc, jac_uptodate, h2)
              IF(jac_uptodate) THEN
                 nfev = nfev + 1
              ELSE
                 nfev = nfev + 2 !! we don't use the jacobian if it is not up to date
              ENDIF
           ENDIF
           !! Check accel for nans
           valid_result = .TRUE.
           checkAccel: DO i = 1,m
              IF(acc(i) /= acc(i)) THEN
                 valid_result = .FALSE.
                 EXIT checkAccel
              END IF
           END DO checkAccel
           IF (valid_result ) THEN
              a = -1.0d+0*MATMUL(acc, fjac)
              CALL DPOTRS('U', n, 1, g, n, a, n, info)
              !!a = -1.0d+0*MATMUL(g,MATMUL(acc,fjac))
           ELSE
              a(:) = 0.0d+0 !! If nans in acc, we will ignore the acceleration term
           ENDIF
        ENDIF

        !! Evaluate at proposed step -- only necessary if av <= avmax
        av = SQRT(DOT_PRODUCT(a,MATMUL(dtd, a))/DOT_PRODUCT(v,MATMUL(dtd,v)))
        IF( av .LE. avmax) THEN
           x_new = x + v + 0.5d+0*a
           CALL func(m,n,x_new,fvec_new)
           nfev = nfev + 1
           Cnew = 0.5d+0*DOT_PRODUCT(fvec_new,fvec_new)
           Cold = C
           valid_result = .TRUE.
           !! Check for nans in fvec_new
           checkfvec_new: DO i = 1,m
              IF(fvec_new(i) /= fvec_new(i)) THEN
                 valid_result = .FALSE.
                 EXIT checkfvec_new
              END IF
           END DO checkfvec_new
           IF (valid_result) THEN  !! If no nans, proceed as normal
              ! update rho and actred
              actred = 1.0d+0 - Cnew/C
              rho = 0.0d+0
              IF(pred_red .NE. 0.0d+0) rho = (1.0d+0 - Cnew/C)/pred_red
              !! Accept or Reject proposed step
              CALL Acceptance(n,C, Cnew, Cbest, ibold, accepted, dtd, v, vold)
           ELSE !! If nans in fvec_new, reject step
              actred = 0.0d+0
              rho = 0.0d+0
              accepted = MIN(accepted - 1, -1)
           ENDIF
        ELSE !! If acceleration too large, then reject
           accepted = MIN(accepted - 1, -1)
        ENDIF
     ELSE !! If matrix factorization fails, reject the proposed step
        accepted = MIN(accepted -1, -1)
     ENDIF

     !! Check Convergence
     IF (converged .EQ. 0) THEN
        CALL convergence_check(m, n, converged, accepted, counter, &
             & C, Cnew, x, fvec, fjac, lam, x_new, &
             & nfev, maxfev, njev, maxjev, naev, maxaev, maxlam, minlam, &
             & artol, Cgoal, gtol, xtol, xrtol, ftol, frtol,cos_alpha)
        IF (converged .EQ. 1 .AND. .NOT. jac_uptodate) THEN
           !! If converged by artol with an out of date jacobian, update the jacoban to confirm true convergence
           converged = 0
           jac_force_update = .TRUE.
        END IF
     ENDIF

     !! Print status
     IF (print_level .EQ. 2 .AND. accepted .GT. 0) THEN
        WRITE(print_unit, *) "  istep, nfev, njev, naev, accepted", istep, nfev, njev, naev, accepted
        WRITE(print_unit, *) "  Cost, lam, delta", C, lam, delta
        WRITE(print_unit, *) "  av, cos alpha", av, cos_alpha
        FLUSH(print_unit)
     ELSEIF(print_level .EQ. 3) THEN
        WRITE(print_unit, *) "  istep, nfev, njev, naev, accepted", istep, nfev, njev, naev, accepted
        WRITE(print_unit, *) "  Cost, lam, delta", C, lam, delta
        WRITE(print_unit, *) "  av, cos alpha", av, cos_alpha
        FLUSH(print_unit)
     ENDIF
     IF (print_level .EQ. 4 .AND. accepted .GT. 0) THEN
        WRITE(print_unit, *) "  istep, nfev, njev, naev, accepted", istep, nfev, njev, naev, accepted
        WRITE(print_unit, *) "  Cost, lam, delta", C, lam, delta
        WRITE(print_unit, *) "  av, cos alpha", av, cos_alpha
        WRITE(print_unit, *) "  x = ", x
        WRITE(print_unit, *) "  v = ", v
        WRITE(print_unit, *) "  a = ", a
        FLUSH(print_unit)
     ELSEIF (print_level .EQ. 5) THEN
        WRITE(print_unit, *) "  istep, nfev, njev, naev, accepted", istep, nfev, njev, naev, accepted
        WRITE(print_unit, *) "  Cost, lam, delta", C, lam, delta
        WRITE(print_unit, *) "  av, cos alpha", av, cos_alpha
        WRITE(print_unit, *) "  x = ", x
        WRITE(print_unit, *) "  v = ", v
        WRITE(print_unit, *) "  a = ", a
        FLUSH(print_unit)
     ENDIF

     ! If converged -- return
     IF(converged .NE. 0) THEN
        exit main
     ENDIF

     IF (accepted .GE. 0) jac_uptodate = .FALSE. !jacobian is now out of date

  END DO main
  ! end main loop

  ! If not converged
  IF(converged .EQ. 0) converged = -1
  niters = istep

  ! Return best fit found
  ! If the method converged, but final x is different from x_best -- what to do?
  x = x_best
  fvec = fvec_best

  IF(print_level .GE. 1) THEN
     WRITE(print_unit,*) "Optimization finished"
     WRITE(print_unit,*) "Results:"
     WRITE(print_unit,*) "  Converged:    ", converged_info(converged), converged
     WRITE(print_unit,*) "  Final Cost: ", 0.5d+0*DOT_PRODUCT(fvec,fvec)
     WRITE(print_unit,*) "  Cost/DOF: ", 0.5d+0*DOT_PRODUCT(fvec,fvec)/(m-n)
     WRITE(print_unit,*) "  niters:     ", istep
     WRITE(print_unit,*) "  nfev:       ", nfev
     WRITE(print_unit,*) "  njev:       ", njev
     WRITE(print_unit,*) "  naev:       ", naev
     FLUSH(print_unit)
  ENDIF

END SUBROUTINE geodesiclm


text:
module ElementsMod
 implicit none
 private
 public :: elementnrQ,abundancesQ,weightsQ,fractionsQ,wfractionsQ,zQ ,&
           meanmwtQ,nelementsQ,elementnamesQ,wfractionQ,InitElements, &
           solarzQ,fractionQ,abundanceQ

  integer, parameter :: nelements=8
  character(len=4),dimension(nelements),parameter :: elementnames = &
             (/"H   ","He  ","C   ","N   ","O   ","Ne  ","Si  ","Fe  "/)
  integer, parameter :: H=1,He=2,metals(6)=(/ 3,4,5,6,7,8 /)

  real :: abundances(nelements)=(/0., 0., 0., 0., 0., 0., 0., 0. /)

!  real, parameter :: abundances(nelements)= &
!                      (/ 1.,.09,3.75e-5,8.7e-6,4.4e-5,2.6e-6,3.2e-6,3.2e-6 /)
!  real, parameter :: abundances(nelements)= &
!                      (/ 1.,.1,3.75e-4,8.7e-5,4.4e-4,2.6e-5,3.2e-5,3.2e-5 /)

  real, parameter :: weights(nelements) = &
                      (/ 1.0079,4.0026,12.011,14.007,15.999,20.18,28.086,55.845/)
  real, save :: fractions(nelements)=0., wfractions(nelements)=0., meanmwt=0.

  real, parameter :: solarz=0.01545

contains

  function elementnrQ(string) result(x)
   character(len=4), intent(in) :: string
   integer :: i,x
   x=0
   do i=1,nelements
    if(elementnames(i).EQ.string) x=i
   enddo
  end function

  subroutine InitElements(metallicity)
   character(len=*), optional :: metallicity

   if(any(abundances.NE.0)) return

   if(present(metallicity)) then
    select case (metallicity)
    case('2.5-solar')
      abundances=(/ 1.,.085,8.28e-4,2.1e-4,1.67e-3,3.e-4,8.88e-5,8.e-5 /)
    case('solar')
      abundances=(/ 1.,.085,3.31e-4,8.3e-5,6.76e-4,1.2e-4,3.55e-5,3.2e-5 /)
    case('0.2-solar')
      abundances=(/ 1.,.085,6.62e-5,1.66e-5,1.33e-4,2.4e-5,7.1e-6,6.4e-6 /)
    case('0.1-solar')
      abundances=(/ 1.,.085,3.31e-5,8.3e-6,6.76e-5,1.2e-5,3.55e-6,3.2e-6 /)
    case('0.02-solar')
      abundances=(/ 1.,.085,6.62e-6,1.66e-6,1.33e-5,2.4e-6,7.1e-7,6.4e-7 /)
    case('0.01-solar')
      abundances=(/ 1.,.065,3.31e-6,8.3e-7,6.76e-6,1.2e-6,3.55e-7,3.2e-7 /)
    case('FeSi-depl')
      abundances=(/ 1.,.085,3.31e-4,8.3e-5,6.76e-4,1.2e-4,1.18e-6,1.07e-6 /)
    case('C-depl')
      abundances=(/ 1.,.085,1.1e-5,8.3e-5,6.76e-4,1.2e-4,3.55e-5,3.2e-5 /)
    case default
      print*,'unknown metallicity (compilation error)'
      stop
    end select
   else
    abundances=(/ 1.,.085,3.31e-4,8.3e-5,6.76e-4,1.2e-4,3.55e-5,3.2e-5 /)
   endif

   fractions=abundances/sum(abundances)
   meanmwt=sum(fractions*weights)
   wfractions=weights*abundances/sum(weights*abundances)
  end subroutine

  function abundancesQ() result(x)
   real :: x(nelements)
   x=abundances
  end function
  function abundanceQ(i) result(x)
   integer :: i
   real :: x
   x=abundances(i)
  end function
  function weightsQ() result(x)
   real :: x(nelements)
   x=weights
  end function
  function fractionsQ() result(x)
   real :: x(nelements)
   x=fractions
  end function
  function wfractionsQ() result(x)
   real :: x(nelements)
   x=wfractions
  end function
  function meanmwtQ() result(x)
   real :: x
   x=meanmwt
  end function
  function nelementsQ() result(x)
   integer :: x
   x=nelements
  end function
  subroutine elementnamesQ(x)
   character(len=4),dimension(nelements) :: x
   x=elementnames
  end subroutine
  function wfractionQ(n) result(x)
   integer :: n
   real :: x
   x=wfractions(n)
  end function
  function fractionQ(n) result(x)
   integer :: n
   real :: x
   x=fractions(n)
  end function
  function zQ() result(x)
   real :: x
   x=sum(wfractions(metals))
  end function
  function solarzQ() result(x)
   real :: x
   x=solarz
  end function

end module ElementsMod

!program test
! use ElementsMod
! real weights(8)
! call InitElements

! print*,meanmwt/fractions(1),wfractions(1)

! weights=wfractionsQ()

! print*,weights
! print*,weights(1),weights(2),sum(weights(3:8))
! print*, zQ()
!end program

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE PCROSS (MXND, MXCORN, MLN, MXLOOP, MAXPRM, NUID,
     &   XN, YN, ZN, LXK, KXL, NXL, LXN, ANGLE, LNODES, BNSIZE, LINKPR,
     &   KPERIM, NODE, NODE1, NODE2, KKKOLD, LLLOLD, NNNOLD, IAVAIL,
     &   NAVAIL, DONE, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, DEV1, LLL,
     &   KKK, NNN, LCORN, NCORN, NLOOP, NEXTN1, KLOOP, GRAPH, VIDEO,
     &   KREG, NOROOM, ERR)
C***********************************************************************

C  SUBROUTINE PCROSS = CHECKS TO SEE IF ANY PERIMETERS CROSS AND HOOKS
C                      THEM TOGETHER IF THEY DO

C***********************************************************************

      COMMON /TIMING/ TIMEA, TIMEP, TIMEC, TIMEPC, TIMEAJ, TIMES

      DIMENSION XN (MXND), YN (MXND), ZN (MXND), NUID (MXND)
      DIMENSION LXK (4, MXND), KXL (2, 3*MXND)
      DIMENSION NXL (2, 3*MXND), LXN (4, MXND)
      DIMENSION ANGLE (MXND), LNODES (MLN, MXND), BNSIZE (2, MXND)
      DIMENSION LCORN (MXCORN)
      DIMENSION NLOOP (MXLOOP), NEXTN1 (MXLOOP), LINKPR (3, MAXPRM)

      CHARACTER*3 DEV1

      LOGICAL DONE,  ERR, NOROOM, DONE1
      LOGICAL GRAPH, BOK, LCROSS, LMATCH
      LOGICAL VIDEO, PMATCH

C  FIND THE FIRST OVERLAPPING LINE STARTING AT THE CURRENT NODE

      CALL GETIME (TIME1)
      ERR = .FALSE.
      PMATCH = .TRUE.

  100 CONTINUE
      if (node1 .eq. 0) return
      N1 = NODE1
      KOUNT = 0

  110 CONTINUE
      N0 = LNODES (2, N1)
      N2 = LNODES (3, N1)
      N3 = LNODES (3, N2)
      KOUNT = KOUNT + 1

C  CHECK FOR COMPLETION

      IF ((N1 .EQ. NODE2) .AND. (KOUNT .GT. 1)) THEN
         GOTO 140
      ELSEIF (KOUNT .GT. NLOOP (1) + 1) THEN
         CALL MESAGE ('** PROBLEMS WITH LOOP CLOSING IN PCROSS **')
         ERR = .TRUE.
         GOTO 140
      ENDIF

C  LOOP THROUGH ALL THE REMAINING PERIMETERS CHECKING FOR CROSSINGS

      IPERIM = LINKPR (2, KPERIM)
  120 CONTINUE

      IF (IPERIM .EQ. KPERIM) THEN
         N1 = N2
         GOTO 110
      ENDIF

      KOUNT2 = 0
      N1TEST = LINKPR (1, IPERIM)

  130 CONTINUE
      N0TEST = LNODES (2, N1TEST)
      N2TEST = LNODES (3, N1TEST)
      N3TEST = LNODES (3, N2TEST)

      CALL INTSCT (XN(N1), YN(N1), XN(N2), YN(N2), XN(N1TEST),
     &   YN(N1TEST), XN(N2TEST), YN(N2TEST), U, W, LCROSS)
      IF (.NOT. LCROSS) THEN
         N1TEST = N2TEST
         KOUNT2 = KOUNT2 + 1
         IF (N1TEST .EQ. LINKPR (1, IPERIM)) THEN
            IPERIM = LINKPR (2, IPERIM)
            GOTO 120
         ELSEIF (KOUNT2 .GT. LINKPR (3, IPERIM)) THEN
            CALL MESAGE ('** PROBLEMS IN PCROSS WITH UNCLOSED '//
     &         'PERIMETER **')
            ERR = .TRUE.
            GOTO 140
         ENDIF
         GOTO 130
      ENDIF

C  AN INTERSECTION HAS OCCURRED.
C  GET THE BEST SEAM FROM THIS INTERSECTION

      IF ((GRAPH) .OR. (VIDEO)) THEN
         CALL RPLOTL (MXND, XN, YN, ZN, NXL, XMIN, XMAX,
     &      YMIN, YMAX, ZMIN, ZMAX, LLL, DEV1, KREG)
         IF (VIDEO) CALL SNAPIT (2)
         IF (GRAPH) THEN
            CALL LCOLOR ('YELOW')
            CALL D2NODE (MXND, XN, YN, N1, N2)
            CALL D2NODE (MXND, XN, YN, N1TEST, N2TEST)
            CALL LCOLOR ('WHITE')
            CALL SFLUSH
         ENDIF
      ENDIF
      CALL MATCH2 (MXND, MLN, XN, YN, NXL, LXN, LNODES, ANGLE, N0, N1,
     &   N2, N3, N0TEST, N1TEST, N2TEST, N3TEST, I1, I2, J1, J2,
     &   KOUNTL, LMATCH, LINKPR (3, IPERIM), NODE, U, W, NLOOP (1),
     &   PMATCH, ERR)
      IF (ERR) GOTO 140
      IF (GRAPH) THEN
         CALL LCOLOR ('PINK ')
         CALL D2NODE (MXND, XN, YN, I1, I2)
         CALL D2NODE (MXND, XN, YN, J1, J2)
         CALL LCOLOR ('WHITE')
         CALL SFLUSH
      ENDIF
      IF (.NOT. LMATCH) THEN
         IF (N1TEST .EQ. LINKPR (1, IPERIM)) THEN
            IPERIM = LINKPR (2, IPERIM)
            GOTO 120
         ELSEIF (KOUNT2 .GT. LINKPR (3, IPERIM)) THEN
            CALL MESAGE ('** PROBLEMS IN PCROSS WITH UNCLOSED '//
     &         'PERIMETER **')
            ERR = .TRUE.
            GOTO 140
         ENDIF
         GOTO 130
      ENDIF

C  NOW CHECK TO SEE IF THE ATTACHMENT WOULD CAUSE
C  LINES ON THE BOUNDARY TO CROSS

      CALL BCROSS (MXND, MLN, XN, YN, ZN, LXK, KXL, NXL, LXN, LNODES,
     &   I1, I2, J1, J2, NLOOP(1), BOK, LLL, XMIN, XMAX, YMIN, YMAX,
     &   ZMIN, ZMAX, DEV1, KREG, ERR)
      IF (ERR) GOTO 140
      IF (.NOT. BOK) THEN
         N1TEST = N2TEST
         KOUNT2 = KOUNT2 + 1
         IF (N1TEST .EQ. LINKPR (1, IPERIM)) THEN
            IPERIM = LINKPR (2, IPERIM)
            GOTO 120
         ELSEIF (KOUNT2 .GT. LINKPR (3, IPERIM)) THEN
            CALL MESAGE ('** PROBLEMS IN PCROSS WITH UNCLOSED '//
     &         'PERIMETER **')
            ERR = .TRUE.
            GOTO 140
         ENDIF
         GOTO 130
      ENDIF

C  NOW THAT THE APPROPRIATE COLLAPSE HAS BEEN FOUND, THE TWO LINES
C  MUST BE JOINED AND THE PERIMETER LINKS RESTABLISHED

      CALL SEW2 (MXND, MLN, NUID, LXK, KXL, NXL, LXN, LNODES,
     &   IAVAIL, NAVAIL, LLL, KKK, NNN, I1, I2, J1, J2, NOROOM, ERR)
      IF ((NOROOM) .OR. (ERR)) GOTO 140

C  UPDATE THE CURRENT NODE

      IF (J1 .EQ. NODE) THEN
         NDUM = NODE
         NODE = I2
         IF (NODE1 .EQ. NDUM) NODE1 = I2
         IF (NODE2 .EQ. NDUM) NODE2 = I2
      ELSEIF (J2 .EQ. NODE) THEN
         NDUM = NODE
         NODE = I1
         IF (NODE1 .EQ. NDUM) NODE1 = I1
         IF (NODE2 .EQ. NDUM) NODE2 = I1
      ENDIF

      NLOOP (1) = NLOOP (1) + LINKPR (3, IPERIM) - 2
      LINKPR (3, KPERIM) = NLOOP (1)
      JPERIM = LINKPR (2, IPERIM)
      IF (JPERIM .EQ. KPERIM) KPERIM = IPERIM
      LINKPR (1, IPERIM) = LINKPR (1, JPERIM)
      LINKPR (2, IPERIM) = LINKPR (2, JPERIM)
      LINKPR (3, IPERIM) = LINKPR (3, JPERIM)
      IF (LINKPR (2, KPERIM) .EQ. KPERIM) LINKPR (2, KPERIM) = 0

C  NOW SMOOTH AND PLOT THE CURRENT MESH

      NNN2 = 1
      CALL GETIME (TIME2)
      TIMEPC = TIMEPC + TIME2 - TIME1
      CALL FILSMO  (MXND, MLN, XN, YN, ZN, LXK, KXL, NXL, LXN, LLL, NNN,
     &   NNN2, LNODES, BNSIZE, NLOOP (1), XMIN, XMAX, YMIN, YMAX, ZMIN,
     &   ZMAX, DEV1, KREG)
      CALL GETIME (TIME1)
      IF ((GRAPH) .OR. (VIDEO)) THEN
         CALL RPLOTL (MXND, XN, YN, ZN, NXL, XMIN, XMAX,
     &      YMIN, YMAX, ZMIN, ZMAX, LLL, DEV1, KREG)
         IF (VIDEO) CALL SNAPIT (2)
         IF (GRAPH) THEN
            CALL LCOLOR ('YELOW')
            CALL D2NODE (MXND, XN, YN, I1, I2)
            CALL LCOLOR ('WHITE')
            CALL SFLUSH
         ENDIF
      ENDIF

C  NOW TRY TO PINCH THE CONNECTION

      CALL LUPANG (MXND, MLN, XN, YN, ZN, LXK, KXL, NXL, LXN, NLOOP (1),
     &   ANGLE, LNODES, I2, LLL, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX,
     &   DEV1, KREG, ERR)
      IF (ERR) GOTO 140

      CALL GETIME (TIME2)
      TIMEPC = TIMEPC + TIME2 - TIME1
      CALL PINCH (MXND, MXCORN, MLN, NUID, XN, YN, ZN, LXK, KXL, NXL,
     &   LXN, ANGLE, LNODES, BNSIZE, NODE, NLOOP (1), KKKOLD, LLLOLD,
     &   NNNOLD, IAVAIL, NAVAIL, DONE1, XMIN, XMAX, YMIN, YMAX, ZMIN,
     &   ZMAX, DEV1, LLL, KKK, NNN, LCORN, NCORN, NODE1, NODE2, GRAPH,
     &   VIDEO, KREG, NOROOM, ERR)
      IF ((NOROOM) .OR. (ERR)) GOTO 140
      IF (LINKPR(2, KPERIM) .NE. 0) GO TO 100
      CALL GETIME (TIME1)

  140 CONTINUE
      CALL GETIME (TIME2)
      TIMEPC = TIMEPC + TIME2 - TIME1
      RETURN

      END

ext: f

max_stars_repo_name: utkarshayachit/seacas

max_issues_repo_name: utkarshayachit/seacas

max_forks_repo_name: utkarshayachit/seacas


text:
Hi! Im a librarian at the Physical Sciences and Engineering Library; I am the librarian for computer science, electrical engineering and physics. I can you help you find stuff, and I buy books in these areas as well. Feel free to come in or email me with reference questions or suggestions.

I work a lot on wiki:WikiPedia:user:phoebe another pretty wellknown wiki too.

contact
IM: phoebeayers
7529948

She also totally http://howwikipediaworks.com/ wrote a book on how wikipedia works!

libraries
There are many libraries in Davis.
 The Yolo Public Library is nicer than ever after their recent makeover and has some great places to study, as well as a lovely selection of pleasure reading.
 The main library on campus is Shields Library, where you can find most subjects. Protip: the recent magazines and newspapers are in the basement.
 the library I work in, the Physical Sciences and Engineering Library, is known as a quiet place to study as well as home of all of the engineering & physical science books
 the Carlson Health Sciences Library down by the stadium/genome building, is another quiet place to study

Wikipedia
There is a Wikipedia meetup group in the Bay Area/Berkeley. Should there be one in Davis?

Food
I use the daviswiki a lot for restaurant info. Some of my favorite cheap meals in Davis:

 Chicken shwarma with fries at Sams Mediterranean  as good as any Ive had anywhere, including the middle east.
 Village bakery pizza  nom. Especially the $5 lunch special of slice + salad
 California (potato) burrito at Taqueria Davis. Consistently outstanding. I love potato burritos.
 Crispy fish tacos at El Mariachi. Hard to eat but worth the challenge. Their veggie burritos are also yummy.
 Mac and cheese from the hot food bar at the coop. Yum. Also, the veggie sushi rolls from the coop.
 orange chicken from House of Chang. A particularly good rendition of this totally bad for you chineseamerican classic.
 Salt bagel and cream cheese from Posh Bagel. A good salt bagel is hard to find. I also often order an avocado & cream cheese bagel sandwich there, which I wish they would just go ahead and put on the menu.
 tofu burger from Burgers and brew  the place gets on my nerves sometimes but their smoked tofu is delicious, as are the fries.
 California crepe from Crepeville; I love the avocado + sundried tomato combo.
 the Malibu from that new sandwich place  its like a pound of pastrami plus avocado. Nontraditional but yummy, if one is in the mood for serious meat.
 Adams Special at Delta of Venus  almost the only thing I order at DoV because its so good.
 Lamb gyro from the gyro truck on campus, which is (fortunately? unfortunately?) mere steps away from my office.
 dumplings from the chinese dumpling formerly fish & chips place
 tom kha gai from Thai Nakorn (it has cabbage, which is bizarre, but its tasty); stirfrys from Taste of Thai.
 tamales from the davis farmers market  my numberone outoftowner meal pick
 popsicles from Jaymes Fat Face
 popcorn from the Davis Varsity  while not a meal per se I am picky about movie theater popcorn, and this is not only by far the best in town but one of the best in the area.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C######################################################################|
      module electric_fields
C----------------------------------------------------------------------|
C     This module is used to calculate characteristic electric fields.
C
C     Functions:
C       - E_ceff_ocer_E_c(.):       normalised effective critical field
C       - E_c(ne, T):               critical electric field
C       - E_D(ne, T):               Dreicer electric field
C
C     Function arguments:
C       - Z:    atomic number of each ion species
C       - Z0:   net charge of plasma for each ion species
C       - nj:   density of each ion species (m**-3)
C       - nSpe: number of ion species
C       - T:    electron temperature (eV)
C       - B:    magnetic field (T)
C       - ne:   free electron density (m**-3)
C
C----------------------------------------------------------------------|
      use double
      use physical_constants, only :
     >  alpha
      use Coulomb_logarithms, only :
     >  ln_Lambda_0, ln_Lambda_c
      use collision_frequencies, only :
     >  calc_nuBar_D01, calc_nuBar_S01

      implicit none

      private

      public ::
     >  E_ceff_over_E_c,
     >  E_c,
     >  E_D

      contains

C======================================================================|
      real(kind=dp) function E_ceff_over_E_c(Z, Z0, nj, nSpecies,
     >  T, B)
C----------------------------------------------------------------------|
C     The effective critical field is calculated as in
C     L. Hesslow et al., Plasma Phys. Control. Fusion 60, 074010 (2018)
C     in Eq. (23)-(24) for a plasma composition with any fully ionised
C     species in combination with partially ionized neon and/or argon.
C
C       Input:
C           Z   : atomic number of each ion species
C           Z0  : net charge of plasma for each ion species
C           nj  : density of each ion species (m**-3)
C           nSpe: number of ion species
C               Z, Z0 and nj must be vectors of length nSpecies
C           T   : temperature (eV)
C           B   : magnetic field (T)
C
C       Output:
C           Eceff/Ec, where Ec = ne e^3 lnLc/(4 pi epsilon_0^2 m_e c^2)
C
C----------------------------------------------------------------------|
        ! ----- Function arguments ------------------------------------|
      integer, intent(in) ::
     >  nSpecies, Z(nSpecies), Z0(nSpecies)

      real(kind=dp), intent(in) ::
     >  B, nj(nSpecies), T

        ! ----- Local variables ---------------------------------------|
      integer ::
     >  i

      real(kind=dp) ::
     >  a0, a1, a2, ne, nuBar_D0, nuBar_D1, nuBar_S0, nuBar_S1,
     >  tau_rad_inv

C----------------------------------------------------------------------|
C     Construct required quantities
C----------------------------------------------------------------------|
        ! ----- Free electron density ---------------------------------|
      ne = sum(Z0(:)*nj(:))

        ! ----- Generalised deflection frequency ----------------------|
      call calc_nuBar_D01(Z, Z0, nj, nSpecies, T, nuBar_D0, nuBar_D1)

        ! ----- Slowing down frequency --------------------------------|
      call calc_nuBar_S01(Z, Z0, nj, nSpecies, T, nuBar_S0, nuBar_S1)

        ! ----- Synchrotron radiation term ----------------------------|
      tau_rad_inv = B**2/(15.44_dp*ln_Lambda_c(ne, T))/(1.e-20_dp*ne)

C----------------------------------------------------------------------|
C     Calculate Eceff/Ec analytically
C----------------------------------------------------------------------|
C     With ideas from Vinodh Bandaru (IPP)

        ! ----- Coefficients of cubic equation ------------------------|
        ! a3*x**3 + a2*x**2 + a1*x + a0 = 0
        ! Note, that a3 = 1 in this case and is thus ignored.
      a0 = log(nuBar_D0/(2*nuBar_S1))

      a2 = -2._dp*(nuBar_S0 + nuBar_S1*(1+nuBar_D1/nuBar_D0)*a0)
      a1 = (.5_dp*a2)**2 - (.7_dp+.4_dp*a0)*alpha*nuBar_D0*nuBar_D1
     >      - nuBar_S1**2
      a0 = -2._dp*tau_rad_inv*nuBar_D0**2

        ! ----- Factors of the solution for x -------------------------|
        ! For the solution, the factors are
        ! a0 -> R = (9*a2*a1 - 27*a0 - 2*a2**3)/54
        ! a1 -> Q = (3*a1 - a2**2)/9
        ! a1 -> D = Q**3 + R**2
      a0 = (9*a2*a1 - 27*a0 - 2*a2**3)/54._dp
      a1 = (3*a1 - a2**2)/9._dp
      a1 = a1**3 + a0**2

        ! ----- Solution for x ----------------------------------------|
        ! The solution is given by:
        ! x = (R + sqrt(D))**(1/3) + (R - sqrt(D))**(1/3) - a2/3
        ! Since it may be that D < 0, complex numbers have to be used
        ! Note, that a1 -> D and a0 -> R
      E_ceff_over_E_c = real(
     >    (a0 + sqrt(cmplx(a1, 0._dp, kind=dp)))**(1._dp/3._dp)
     >  + (a0 - sqrt(cmplx(a1, 0._dp, kind=dp)))**(1._dp/3._dp)
     >  - a2/3._dp)

      return
      end function E_ceff_over_E_c
C======================================================================|

C======================================================================|
      real(kind=dp) function E_c(ne, T)
C----------------------------------------------------------------------|
C     Calculates the critical electric field E_c:
C       E_c = me*c/e/tauc = ne*e**3*ln Lambda_c/(4*pi*eps0**2*me*c**2).
C
C     Input:
C       ne  : free electron density (m**-3)
C       T   : electron temperature (eV)
C
C     Output:
C       E_c : critical electric field (V)
C
C----------------------------------------------------------------------|
        ! ----- Function arguments ------------------------------------|
      real(kind=dp), intent(in) ::
     >  ne, T

C----------------------------------------------------------------------|
      E_c = 5.09909908817339515e-23_dp*ne*ln_Lambda_c(ne, T)

      return
      end function E_c
C======================================================================|

C======================================================================|
      real(kind=dp) function E_D(ne, T)
C----------------------------------------------------------------------|
C     Calculates the Dreicer electric field E_D:
C       E_D = ne*e**3*ln Lambda/(4*pi*eps0**2*T).
C
C     Input:
C       ne  : free electron density (m**-3)
C       T   : electron temperature (eV)
C
C     Output:
C       E_D : Dreicer electric field (V)
C
C----------------------------------------------------------------------|
        ! ----- Function arguments ------------------------------------|
      real(kind=dp), intent(in) ::
     >  ne, T

C----------------------------------------------------------------------|
      E_D = 2.6056342609758366e-17_dp*ne*ln_Lambda_0(ne, T)/T

      return
      end function E_D
C======================================================================|

      end module electric_fields
C######################################################################|

max_stars_repo_name: o-linder/runawayelectrongeneration

max_issues_repo_name: o-linder/runawayelectrongeneration

max_forks_repo_name: o-linder/runawayelectrongeneration


text:
great_ne(ee1,kk1).
great_ne(cc1,v1).
great_ne(t1,ii1).
great_ne(jj1,kk1).
great_ne(y1,f1).
great_ne(bb1,ii1).
great_ne(z1,b1).
great_ne(u1,f1).
great_ne(ii1,ff1).
great_ne(p1,y1).
great_ne(t1,ff1).
great_ne(hh1,x1).
great_ne(ii1,a1).
great_ne(t1,x1).
great_ne(q1,t1).
great_ne(cc1,aa1).
great_ne(q1,dd1).
great_ne(y1,c1).
great_ne(k1,a1).
great_ne(aa1,bb1).
great_ne(jj1,x1).
great_ne(u1,jj1).
great_ne(q1,cc1).
great_ne(w1,aa1).
great_ne(bb1,ee1).
great_ne(h1,ff1).
great_ne(y1,p1).
great_ne(b1,kk1).
great_ne(i1,q1).
great_ne(aa1,cc1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
c$Header: /data/petsun4/data1/src_solaris/imglin/RCS/grad_opt.f,v 1.1 2007/09/05 04:05:45 avi Exp $
c$Log: grad_opt.f,v $
c Revision 1.1  2007/09/05  04:05:45  avi
c Initial revision
c
      subroutine grad_opt(imgt,nx,ny,nz,mmppix,center,param,x0,grad,xopt,ierr)
      real*4 imgt(nx,ny,nz)
      real*4 mmppix(3),center(3),x0(3),grad(3),xopt(3)
      real*4 param(16)
c     param(01:03)	normal to osculating plane (used by edge_trace)
c     param(04)		del arclength              (used by edge_trace)
c     param(05)		sigma eval
      real*4 coef(0:2),array(3,3)
      parameter (nk=3)
      parameter (npts=2*nk+1)
      real*4 u(-nk:nk),v(-nk-1:nk+1),dv(-nk:nk)
      real*4 x(3),unorm(3),e(3)
      logical*4 lfar,liter
      logical*4 ldebug/.false./

c     copy x0 into xopt	as inital guess
      do i=1,3
        xopt(i)=x0(i)
      enddo

      do 10 ncycle=1,5
      if(ldebug)write (*,"('x0 grad_opt ',3f10.4)")x0

      nnn=0
      niter=5
   83 liter=.false.
      lfar=.false.

c     compute gradient
      do i=1,3
        do k=1,3
          x(k)=xopt(k)
        enddo
        x(i)=xopt(i)+0.5*param(05)
        call imgvalx(imgt,nx,ny,nz,center,mmppix,x,v2,lslice)
        if(lslice.le.0)goto 99
        x(i)=xopt(i)-0.5*param(05)
        call imgvalx(imgt,nx,ny,nz,center,mmppix,x,v1,lslice)
        if(lslice.le.0)goto 99
        grad(i)=(v2-v1)/param(05)
      enddo

c     compute unit normal
      a=sqrt(dot(grad,grad))
      do i=1,3
        unorm(i)=-grad(i)/a
      enddo
c     write (*, "('unit normal',3f10.4)")unorm

      do 9 k=-nk-1,nk+1			! sample image along normal
      do i=1,3
        x(i)=xopt(i)-param(05)*unorm(i)*(float(k)/float(nk))
      enddo
      call imgvalx(imgt,nx,ny,nz,center,mmppix,x,v(k),lslice)
      if(lslice.le.0)goto 99
c     write (*, "('imgvalx x,v',3f7.2,f9.1)",x,v(k)
    9 continue
      do 8 k=-nk,nk			! differentiate
      u(k)=float(k)/float(nk)
    8 dv(k)=v(k+1)-v(k-1)
c     write (*,"('dv',10f9.1)")dv

      call polfit(u,dv,npts,coef,array,3,chisqr)
c     write (*,"('coef ',3g12.4)",coef
      t=-.5*coef(1)/coef(2)
      if(coef(2).gt.0..or.abs(t).gt.0.5)then
        t0=t
        t=sign(.5,coef(1))
        if(ldebug)write(*,"('t ',f10.6,' ->',f10.6)")t0,t
        lfar=.true.
      endif
      liter=liter.or.abs(t).gt.0.02
      do i=1,3
        xopt(i)=xopt(i)-t*param(05)*unorm(i)
      enddo
      nnn=nnn+1
      if(nnn.gt.20)then
        write(*,"('failed convergence')")
        ierr=1
        return
      endif
      if(.not.lfar)then
        niter=niter-1
      endif
      if(liter.and.niter.gt.0)goto 83

c     write (*,"('xopt       ',3f10.4)")xopt
      do i=1,3
        e(i)=x0(i)-xopt(i)
      enddo
      dist=dot(unorm,e)
      do i=1,3
        xopt(i)=x0(i)-dist*unorm(i)
      enddo
c     write (*,"('xopt new   ',3f10.4)")xopt
   10 continue
      ierr=0
      return

   99 write (*,"('grad_opt: image boundary reached')")
      ierr=-1
      return
      end

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
!##############################################################################
!# ****************************************************************************
!# <name> nonstationaryoptcsolver </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module realises a fully coupled time dependent solver for the
!# coupled Navier-Stokes optimal control problem.
!#
!# The following routines can be found here:
!#
!# 1.) cc_initParTimeDependence
!#     -> Initialise the parameters of the time dependent solver from DAT file
!#        parameters.
!#
!# 2.) cc_solveNonstationaryDirect
!#     -> Solves the space-time coupled system.
!# </purpose>
!##############################################################################

module nonstationaryoptcsolver

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

  use collection
  use convection

  use basicstructures
  use user_callback

  use spacepreconditioner
  use spacepreconditionerinit
  use stationaryoptcsolver
  use timeanalysis
  use spatialbc
  use spacediscretisation
  use postprocessing
  use nonlinearoneshotspacetimesolver
  use spacediscretisation

  implicit none

contains

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

!<subroutine>

  subroutine cc_initParTimeDependence (rproblem,ssection,rparams)

!<description>
  ! Initialises parameters in the problem structure according to whether the
  ! simulation is time dependent or not. Initialises the time stepping scheme,
  ! start proceduce, error bounds, etc.
  !
  ! Note: This will not allocate an memory but only initialise the parameters
  ! in rproblem according to the parameters in rparams from the DAT file.
!</description>

!<input>
  ! A parameter list from a DAT file containing parameters that configure the
  ! time dependence.
  type(t_parlist), intent(IN) :: rparams

  ! The name of the section in the parameter list containing the parameters
  ! for the time dependent simulation.
  character(LEN=*), intent(IN) :: ssection
!</input>

!<inputoutput>
  ! The problem structure to be initialised.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!</subroutine>

    ! Fetch the parameters. Initialise with standard settings if they don't
    ! exist.
    call parlst_getvalue_int (rparams,ssection,'itimedependence',  &
        rproblem%itimedependence, 0)
    call parlst_getvalue_int (rparams,ssection,'niterations',     &
        rproblem%rtimedependence%niterations, 1000)
    call parlst_getvalue_double (rparams,ssection,'dtimestart',   &
        rproblem%rtimedependence%dtimeInit, 0.0_DP)
    call parlst_getvalue_double (rparams,ssection,'dtimemax',     &
        rproblem%rtimedependence%dtimeMax, 20.0_DP)
    call parlst_getvalue_int (rparams, ssection,'ctimeStepScheme', &
          rproblem%rtimedependence%ctimeStepScheme, 0)
    call parlst_getvalue_double (rparams,ssection,'dtimeStepTheta',    &
        rproblem%rtimedependence%dtimeStepTheta, 1.0_DP)

  end subroutine

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

!<subroutine>

  subroutine cc_solveNonstationaryDirect (rproblem,rvector,ipropagateStartVector)

!<description>
  ! Solve the nonstationary optimal control problem. This allocates
  ! memory for the solution vector and solves the space-time coupled
  ! system.
  ! This is the 'single-grid' variant of the coupled solver, i.e.
  ! it calls the solver once without any multigrid schemes in space/time.
!</description>

!<input>
  ! OPTIONAL: A block vector that serves as stationary initial condition.
  type(t_vectorBlock), intent(IN), optional :: rvector

  ! OPTIONAL: If a start vector is present, this specifies if the start vector
  ! is propagated to all timesteps.
  integer, intent(in), optional :: ipropagateStartVector
!</input>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!</subroutine>

    type(t_ccoptSpaceTimeDiscretisation), dimension(:), pointer :: RspaceTimeDiscr
    type(t_spacetimeVector) :: rx,rd,rb
    integer :: i,ispacelevelcoupledtotimelevel,cspaceTimeSolverType
    integer :: nmaxSimulRefLevel,icurrentspace,icurrenttime
    real(DP) :: dspacetimeRefFactor,dcurrentfactor
    integer :: ctypePreconditioner
    integer(I32) :: TIMENLMIN,TIMENLMAX,istep
    integer, dimension(:,:), allocatable :: Ispacetimelevel

    call output_lbrk()

    ! Get the minimum and maximum time level from the parameter list
    call parlst_getvalue_int (rproblem%rparamList,'TIME-DISCRETISATION',&
                              'TIMENLMIN',TIMENLMIN,1)
    call parlst_getvalue_int (rproblem%rparamList,'TIME-DISCRETISATION',&
                              'TIMENLMAX',TIMENLMAX,1)

    ! Values <= 0 result in a relativ minimum level.
    if (TIMENLMIN .le. 0) TIMENLMIN = TIMENLMAX+TIMENLMIN

    ! Allocate memory for the space-time discretisation structures
    ! on all levels. Initialis the maximum level such that it represents
    ! the finest time level as well as the finest space level.
    allocate(RspaceTimeDiscr(1:TIMENLMAX))

    call output_line ('Initialising space-time discretisation of the maximum level...')

    call sptidis_initDiscretisation (rproblem,TIMENLMAX,&
        rproblem%NLMAX,RspaceTimeDiscr(TIMENLMAX))

    call sptivec_initVector (rx,RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr,&
        RspaceTimeDiscr(TIMENLMAX)%p_rlevelInfo%rdiscretisation)
    call sptivec_initVector (rb,RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr,&
        RspaceTimeDiscr(TIMENLMAX)%p_rlevelInfo%rdiscretisation)
    call sptivec_initVector (rd,RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr,&
        RspaceTimeDiscr(TIMENLMAX)%p_rlevelInfo%rdiscretisation)

    ! If a start vector is given, propagate it to all timesteps.
    if (present(rvector)) then
      ! Initialise the solution at timestep 1.
      call sptivec_setTimestepData(rx,1,rvector)

      ! Eventually propagate the initial solution to all timesteps.
      if (present(ipropagateStartVector)) then
        if (ipropagateStartVector .gt. 0) then
          call output_line ('Propagating start vector...')
          do istep = 2,rx%NEQtime
            call sptivec_setTimestepData(rx,istep,rvector)
          end do
        end if
      end if
    end if

    call output_line ('Reading target flow...')

    ! Read the target flow -- stationary or nonstationary
    call cc_initTargetFlow (rproblem,&
        RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr%dtimeInit,&
        RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr%dtimeMax,&
        RspaceTimeDiscr(TIMENLMAX)%rtimeDiscr%nintervals)

    ! Figure out which type of solver we should use to solve
    ! the problem.
    call parlst_getvalue_int (rproblem%rparamList,'TIME-SOLVER',&
                              'cspaceTimeSolverType',cspaceTimeSolverType,0)

    ! Get the preconditioner type which is to be used in the nonlinear
    ! solver.
    call parlst_getvalue_int (rproblem%rparamList,'TIME-SOLVER',&
                              'ctypePreconditioner',ctypePreconditioner,1)

    ! Call the solver for the space/time coupled system.
    select case (cspaceTimeSolverType)
    case (0)
      ! 1-level (in time) Gauss elimination solver.
      call output_line ('Invoking Gauss elimination solver...')
      call cc_solveSupersystemDirect (rproblem, &
          RspaceTimeDiscr(TIMENLMAX), rx, rb, rd)

    case (1)
      ! 1-level (in time) defect correction solver with a preconditioner
      ! in space.
      !
      ! Call the defect correction solver
      call output_line ('Invoking defect correction solver...')
      call cc_solveSupersystemDefCorr (rproblem, &
          RspaceTimeDiscr(TIMENLMAX), rx, rb, rd, ctypePreconditioner)

    case (2)
      ! Space-time multigrid solver in space and/or time.
      !
      ! Should we couple space and time coarsening/refinement?
      call parlst_getvalue_int (rproblem%rparamList,'TIME-MULTIGRID',&
          'ispacelevelcoupledtotimelevel',ispacelevelcoupledtotimelevel,1)

      ! Parameters of the refinement?
      call parlst_getvalue_int (rproblem%rparamList,'TIME-MULTIGRID',&
          'nmaxSimulRefLevel',nmaxSimulRefLevel,0)

      if (nmaxSimulRefLevel .le. 0) then
        nmaxSimulRefLevel = TIMENLMAX+nmaxSimulRefLevel
      end if
      nmaxSimulRefLevel = min(TIMENLMAX,max(TIMENLMIN,nmaxSimulRefLevel))

      call parlst_getvalue_double (rproblem%rparamList,'TIME-MULTIGRID',&
          'dspacetimeRefFactor',dspacetimeRefFactor,1.0_DP)

      call output_line ('Initialising space-time discretisation of the lower levels...')

      ! Initialise the supersystem on all levels below the topmost one
      select case (ispacelevelcoupledtotimelevel)
      case (0)
        ! Everywhere the same space level
        do i=TIMENLMIN,TIMENLMAX-1
          call sptidis_initDiscretisation (rproblem,i,&
              rproblem%NLMAX,RspaceTimeDiscr(i))
        end do

      case (1)
        ! Space level NLMAX-i = Time level TIMENLMAX-i.
        !
        ! But this is modified by nmaxSimulRefLevel and dspacetimeRefFactor!
        ! From level nmaxSimulRefLevel on, we have to use the maximum
        ! space level. Below of that, on every space refinement we have
        ! to do dspacetimeRefFactor time refinements (rounded).
        !
        ! To get a formula for that, one might get crazy. So we calculate
        ! the levels in advance. Allocate an array Ispacetimelevel
        ! and fill it: Ispacetimelevel(1,:) all the space levels,
        ! Ispacetimelevel(2,:) all the corresponding time levels.
        allocate(Ispacetimelevel(2,TIMENLMAX))

        icurrentspace = rproblem%NLMAX
        icurrenttime = TIMENLMAX
        dcurrentfactor = 0.0_DP
        do i=TIMENLMAX, TIMENLMIN,-1
          Ispacetimelevel(1,i) = icurrentspace
          Ispacetimelevel(2,i) = icurrenttime

          ! Reduce space and/or time level.
          if (dspacetimeRefFactor .ge. 1.0_DP) then
            ! Time coarsened as usual, space probably delayed
            icurrenttime = icurrenttime - 1
            if (i .le. nmaxSimulRefLevel) then
              ! Space coarsening allowed.
              dcurrentfactor = dcurrentfactor + 1.0_DP/dspacetimeRefFactor
              if (dcurrentfactor .ge. 1.0_DP) then
                ! Coarsening in space
                if (icurrentspace .gt. rproblem%NLMIN) &
                    icurrentspace = icurrentspace - 1
                dcurrentfactor = mod(dcurrentfactor,1.0_DP)
              end if
              ! Otherwise no coarsening in space.
            end if
          else
            ! Space coarsened as usual, time probably delayed
            if (icurrentspace .gt. rproblem%NLMIN) &
                icurrentspace = icurrentspace - 1
            if (i .le. nmaxSimulRefLevel) then
              ! Time coarsening allowed.
              dcurrentfactor = dcurrentfactor + dspacetimeRefFactor
              if (dcurrentfactor .ge. 1.0_DP) then
                ! Coarsening in space
                if (icurrenttime .gt. TIMENLMIN) icurrenttime = icurrenttime - 1
                dcurrentfactor = mod(dcurrentfactor,1.0_DP)
              end if
              ! Otherwise no coarsening in time.
            end if
          end if
        end do

        ! Now initialise the level hierarchy
        do i=TIMENLMIN,TIMENLMAX-1
          !CALL sptidis_initDiscretisation (rproblem,i,&
          !    MAX(rproblem%NLMIN,rproblem%NLMAX-(TIMENLMAX-i)),&
          !    RspaceTimeDiscr(i))
          call sptidis_initDiscretisation (rproblem,Ispacetimelevel(2,i),&
              Ispacetimelevel(1,i),RspaceTimeDiscr(i))
        end do

        deallocate (Ispacetimelevel)

      case (2)
        ! Space level NLMAX-i = Time level TIMENLMAX-i,
        ! but no restriction in time
        do i=TIMENLMIN,TIMENLMAX-1
          call sptidis_initDiscretisation (rproblem,TIMENLMAX,&
              max(rproblem%NLMIN,rproblem%NLMAX-(TIMENLMAX-i)),&
              RspaceTimeDiscr(i))
        end do
      end select

      ! Print some statistical information about the discretisation
      call output_lbrk ()
      call output_line ('Triangulation')
      call output_line ('-------------')
      do i=rproblem%NLMIN,rproblem%NLMAX
        call tria_infoStatistics (rproblem%RlevelInfo(i)%rtriangulation,&
          i .eq. rproblem%NLMIN,i)
      end do

      call output_lbrk ()
      call output_line ('Time discretisation')
      call output_line ('-------------------')
      do i=TIMENLMIN,TIMENLMAX
        call output_line ('Level: '//sys_siL(i,10))
        call sptidis_infoDiscretisation (RspaceTimeDiscr(i))
        call output_lbrk ()
      end do

      ! Call the multigrid solver to solve on all these levels
      call output_line ('Invoking nonlinear space-time solver...')
      call cc_solveSupersystemMultigrid (rproblem, &
          RspaceTimeDiscr(TIMENLMIN:TIMENLMAX), rx, rb, rd, ctypePreconditioner)
    end select

    call output_line ('Postprocessing...')

    ! POSTPROCESSING
    !
    call cc_spacetimepostproc (rproblem,RspaceTimeDiscr(TIMENLMAX),rx)

    call output_line ('Cleaning up...')

    call cc_doneTargetFlow (rproblem)

    do i=TIMENLMIN,TIMENLMAX-1
      call sptidis_doneDiscretisation (RspaceTimeDiscr(i))
    end do

    deallocate(RspaceTimeDiscr)

    call sptivec_releaseVector (rb)
    call sptivec_releaseVector (rd)
    call sptivec_releaseVector (rx)

  end subroutine

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

!<subroutine>

  subroutine cc_spacetimepostproc (rproblem,rspacetimediscr,rvector)

!<description>
  ! Performs postprocessing of the space time solution vector rvector.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT) :: rproblem
!</inputoutput>

!<input>
  ! The space-time solution vector on which postprocessing should be applied.
  type(t_spacetimevector), intent(IN) :: rvector

  ! The space-time discretisation structure belonging to that vector.
  type(t_ccoptSpaceTimeDiscretisation), intent(IN) :: rspaceTimeDiscr
!</input>

!</subroutine>

    ! local variables
    integer :: i
    character(SYS_STRLEN) :: stemp,sfileName
    type(t_vectorBlock) :: rvectorTmp
    real(dp) :: dtime

    ! Create a temp vector
    call lsysbl_createVecBlockByDiscr (&
        rspacetimediscr%p_rlevelInfo%rdiscretisation,rvectorTmp,.true.)

    ! Postprocessing of all solution vectors.
    do i = 0,tdiscr_igetNDofGlob(rvector%p_rtimeDiscretisation)-1

      dtime = &
          rvector%p_rtimeDiscretisation%dtimeInit + i * rvector%p_rtimeDiscretisation%dtstep

      ! Evaluate the space time function in rvector in the point
      ! in time dtime. Independent of the discretisation in time,
      ! this will give us a vector in space.
      !CALL sptivec_getTimestepData (rx, i, rvectorTmp)
      call tmevl_evaluate(rvector,dtime,rvectorTmp)

      call cc_postprocessingNonstat (rproblem,i,dtime,rvectorTmp)

    end do

    ! If there's a filename for the solution file given,
    ! write out the current solution.
    call parlst_getvalue_string (rproblem%rparamList, 'TIME-POSTPROCESSING', &
                                 'sfinalSolutionFileName', stemp, '')
    read (stemp,*) sfileName
    if (sfilename .ne. '') then
      stemp = '('''//trim(adjustl(sfilename))//'.'',I5.5)'
      call sptivec_saveToFileSequence(rvector,stemp,.true.)
    end if

    call lsysbl_releaseVector(rvectorTmp);

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


  MODULE pdcomp_module

  implicit none

  private
  public :: pidcomp

  CONTAINS

      subroutine pidcomp(dt,xh,rxh,arh1,arh2,uh,xf,rxf,arf1,arf2,uf,vh,vf,          &
                         gz,rgz,gzu,gzv,mh,rmh,mf,rmf,rds,rdsf,c1,c2,f2d,wprof,     &
                         pi0,th0,rth0,thv0,qv0,qc0,qi0,rho0,rr0,rf0,rrf0,u0,v0,     &
                         dum1,dum2,dum3,dum4,dum5,dum6,dum7,dum8,divx,              &
                         u3d,rru,uten,uten1,v3d,rrv,vten,vten1,w3d,rrw,wten,wten1,  &
                         rho,pp3d,th3d,q3d,udiag,vdiag,wdiag,pdiag,                 &
                         cfb,cfa,cfc,ad1,ad2,pdt,lgbth,lgbph,rhs,trans)

      use input
      use constants
      use adv_module
      use poiss_module

      implicit none

      real, intent(in) :: dt
      real, intent(in), dimension(ib:ie) :: xh,rxh,arh1,arh2,uh
      real, intent(in), dimension(ib:ie+1) :: xf,rxf,arf1,arf2,uf
      real, intent(in), dimension(jb:je) :: vh
      real, intent(in), dimension(jb:je+1) :: vf
      real, intent(in), dimension(itb:ite,jtb:jte) :: gz,rgz,gzu,gzv
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: mh,rmh
      real, intent(in), dimension(ib:ie,jb:je,kb:ke+1) :: mf,rmf
      real, intent(in), dimension(kb:ke) :: rds
      real, intent(in), dimension(kb:ke+1) :: rdsf
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: c1,c2
      real, intent(in), dimension(ib:ie,jb:je) :: f2d
      real, intent(in), dimension(kb:ke) :: wprof
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: pi0,th0,rth0,thv0,qv0,qc0,qi0,rho0,rr0,rf0,rrf0
      real, intent(in), dimension(ib:ie+1,jb:je,kb:ke) :: u0
      real, intent(in), dimension(ib:ie,jb:je+1,kb:ke) :: v0
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke) :: dum1,dum2,dum3,dum4,dum5,dum6,dum7,dum8,divx
      real, intent(in), dimension(ib:ie+1,jb:je,kb:ke) :: u3d
      real, intent(in), dimension(ib:ie,jb:je+1,kb:ke) :: v3d
      real, intent(in), dimension(ib:ie,jb:je,kb:ke+1) :: w3d
      real, intent(inout), dimension(ib:ie+1,jb:je,kb:ke) :: rru,uten,uten1
      real, intent(inout), dimension(ib:ie,jb:je+1,kb:ke) :: rrv,vten,vten1
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke+1) :: rrw,wten,wten1
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: rho,pp3d,th3d
      real, intent(in), dimension(ibm:iem,jbm:jem,kbm:kem,numq) :: q3d
      real, intent(inout) , dimension(ibdv:iedv,jbdv:jedv,kbdv:kedv,nudiag) :: udiag
      real, intent(inout) , dimension(ibdv:iedv,jbdv:jedv,kbdv:kedv,nvdiag) :: vdiag
      real, intent(inout) , dimension(ibdv:iedv,jbdv:jedv,kbdv:kedv,nwdiag) :: wdiag
      real, intent(inout) , dimension(ibdp:iedp,jbdp:jedp,kbdp:kedp,npdiag) :: pdiag
      real, intent(in), dimension(ipb:ipe,jpb:jpe,kpb:kpe) :: cfb
      real, intent(in), dimension(kpb:kpe) :: cfa,cfc,ad1,ad2
      complex, intent(inout), dimension(ipb:ipe,jpb:jpe,kpb:kpe) :: pdt,lgbth,lgbph
      complex, intent(inout), dimension(ipb:ipe,jpb:jpe) :: rhs,trans

      integer :: i,j,k,n
      real :: tem,dudz,dvdz
      double precision :: pavg1,pavg2

      !--------------------------------------------------------------!
      ! buoyancy pressure:

        IF( imoist.eq.1 )THEN
          ! buoyancy (with moisture terms):

          !$omp parallel do default(shared)  &
          !$omp private(i,j,k,n)
          do j=1,nj
            !---
            do k=1,nk
            do i=1,ni
              dum1(i,j,k) = 0.0
            enddo
            enddo
            !---
            do n=nql1,nql2
            do k=1,nk
            do i=1,ni
              dum1(i,j,k) = dum1(i,j,k)+q3d(i,j,k,n)
            enddo
            enddo
            enddo
            !---
            if( iice.eq.1 )then
            do n=nqs1,nqs2
            do k=1,nk
            do i=1,ni
              dum1(i,j,k) = dum1(i,j,k)+q3d(i,j,k,n)
            enddo
            enddo
            enddo
            endif
            !---
            do k=1,nk
            do i=1,ni
              dum6(i,j,k) = g*( th3d(i,j,k)*rth0(i,j,k) + repsm1*(q3d(i,j,k,nqv)-qv0(i,j,k)) - (dum1(i,j,k)-qc0(i,j,k)-qi0(i,j,k)) )
            enddo
            enddo
            !---
            do k=2,nk
            do i=1,ni
              wten(i,j,k)=(c1(i,j,k)*dum6(i,j,k-1)+c2(i,j,k)*dum6(i,j,k))
            enddo
            enddo
            !---
            do i=1,ni
              wten(i,j,1) = 0.0
              wten(i,j,nk+1) = 0.0
            enddo
            !---
          enddo

        ELSE
          ! buoyancy (dry):

          !$omp parallel do default(shared)  &
          !$omp private(i,j,k)
          do j=1,nj
            !---
            do k=1,nk
            do i=1,ni
              dum6(i,j,k) = g*th3d(i,j,k)*rth0(i,j,k)
            enddo
            enddo
            !---
            do k=2,nk
            do i=1,ni
              wten(i,j,k)=(c1(i,j,k)*dum6(i,j,k-1)+c2(i,j,k)*dum6(i,j,k))
            enddo
            enddo
            !---
            do i=1,ni
              wten(i,j,1) = 0.0
              wten(i,j,nk+1) = 0.0
            enddo
            !---
          enddo

        ENDIF

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj+1
        do i=1,ni+1
          uten(i,j,k) = 0.0
          vten(i,j,k) = 0.0
          divx(i,j,k) = 0.0
        enddo
        enddo
        enddo

        call   poiss(uh,vh,mh,rmh,mf,rmf,pi0,thv0,rho0,rf0,    &
                     dum2,dum3,dum4,dum5,divx,pdiag(ibdp,jbdp,kbdp,1),uten,vten,wten,  &
                     cfb,cfa,cfc,ad1,ad2,pdt,lgbth,lgbph,rhs,trans,dt)

        ! Adjust mean pressure ... set average value at upper-most level to
        ! mean of actual pp3d field:

        pavg1 = 0.0
        pavg2 = 0.0

        do j=1,nj
        do i=1,ni
          pavg1 = pavg1+pp3d(i,j,nk)
          pavg2 = pavg2+pdiag(i,j,nk,1)
        enddo
        enddo

        pavg1 = pavg1/dble(nx*ny)
        pavg2 = pavg2/dble(nx*ny)

        tem = sngl(pavg1-pavg2)

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          pdiag(i,j,k,1) = pdiag(i,j,k,1)+tem
        enddo
        enddo
        enddo

      !--------------------------------------------------------------!
      ! linear dynamic pressure:

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k,dudz,dvdz)
        do k=2,nk
          dudz = -0.5*rf0(1,1,k)*(u0(1,1,k)-u0(1,1,k-1))*rdz*mf(1,1,k)
          dvdz = -0.5*rf0(1,1,k)*(v0(1,1,k)-v0(1,1,k-1))*rdz*mf(1,1,k)
          do j=1,nj+1
          do i=1,ni+1
            dum1(i,j,k) = dudz*(w3d(i-1,j,k)+w3d(i,j,k))
            dum2(i,j,k) = dvdz*(w3d(i,j-1,k)+w3d(i,j,k))
          enddo
          enddo
        enddo

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do j=1,nj+1
        do i=1,ni+1
          dum1(i,j,1) = 0.0
          dum1(i,j,nk+1) = 0.0
          dum2(i,j,1) = 0.0
          dum2(i,j,nk+1) = 0.0
        enddo
        enddo

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k,tem)
        do k=1,nk
          tem = 0.5*rr0(1,1,k)
          do j=1,nj+1
          do i=1,ni+1
            uten(i,j,k) = tem*(dum1(i,j,k)+dum1(i,j,k+1))
            vten(i,j,k) = tem*(dum2(i,j,k)+dum2(i,j,k+1))
          enddo
          enddo
          do j=1,nj
          do i=1,ni
            wten(i,j,k) = 0.0
            divx(i,j,k) = 0.0
          enddo
          enddo
        enddo

        call   poiss(uh,vh,mh,rmh,mf,rmf,pi0,thv0,rho0,rf0,    &
                     dum2,dum3,dum4,dum5,divx,pdiag(ibdp,jbdp,kbdp,2),uten,vten,wten,  &
                     cfb,cfa,cfc,ad1,ad2,pdt,lgbth,lgbph,rhs,trans,dt)

        ! Adjust mean pressure ... set average value at upper-most level to zero:

        pavg2 = 0.0

        do j=1,nj
        do i=1,ni
          pavg2 = pavg2+pdiag(i,j,nk,2)
        enddo
        enddo

        pavg2 = pavg2/dble(nx*ny)

        tem = sngl(0.0d0-pavg2)

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          pdiag(i,j,k,2) = pdiag(i,j,k,2)+tem
        enddo
        enddo
        enddo

      !--------------------------------------------------------------!
      ! nonlinear dynamic pressure:

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        DO k=1,nk
          ! store perturbation u,v in uten1,vten1 arrays:
          do j=jb,je
          do i=ib,ie+1
            rru(i,j,k)=rho0(1,1,k)*u3d(i,j,k)
            uten1(i,j,k)=u3d(i,j,k)-u0(i,j,k)
          enddo
          enddo
          do j=jb,je+1
          do i=ib,ie
            rrv(i,j,k)=rho0(1,1,k)*v3d(i,j,k)
            vten1(i,j,k)=v3d(i,j,k)-v0(i,j,k)
          enddo
          enddo
          IF(k.eq.1)THEN
            do j=0,nj+1
            do i=0,ni+1
              rrw(i,j,   1) = 0.0
              rrw(i,j,nk+1) = 0.0
            enddo
            enddo
          ELSE
            do j=0,nj+1
            do i=0,ni+1
              rrw(i,j,k)=rf0(1,1,k)*w3d(i,j,k)
            enddo
            enddo
          ENDIF
        ENDDO

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+1
        do i=0,ni+1
          divx(i,j,k)=( (rru(i+1,j,k)-rru(i,j,k))*rdx*uh(i)        &
                       +(rrv(i,j+1,k)-rrv(i,j,k))*rdy*vh(j) )      &
                       +(rrw(i,j,k+1)-rrw(i,j,k))*rdz*mh(1,1,k)
          if(abs(divx(i,j,k)).lt.smeps) divx(i,j,k)=0.0
          uten(i,j,k) = 0.0
          vten(i,j,k) = 0.0
          wten(i,j,k) = 0.0
        enddo
        enddo
        enddo

        call advu(nrkmax,arh1,arh2,uh,xf,rxf,arf1,arf2,uf,vh,gz,rgz,gzu,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,dum5,dum6,dum7,divx, &
                   rru,uten1,uten,rrv,rrw,rdsf,c1,c2,rho, dt  ,.false.,udiag,wprof)
        call advv(nrkmax,xh,rxh,arh1,arh2,uh,xf,vh,vf,gz,rgz,gzv,mh,rho0,rr0,rf0,rrf0,dum1,dum2,dum3,dum4,dum5,dum6,dum7,divx, &
                   rru,rrv,vten1,vten,rrw,rdsf,c1,c2,rho, dt  ,.false.,vdiag,wprof)
        call   advw(nrkmax,xh,rxh,arh1,arh2,uh,xf,vh,gz,rgz,mh,mf,rho0,rr0,rf0,rrf0,  &
                    dum1,dum2,dum3,dum4,dum5,dum6,dum7,divx,                       &
                    rru,rrv,rrw,w3d  ,wten,rds,rdsf,c1,c2,rho, dt  ,               &
                    .false.,wdiag,hadvordrv,vadvordrv,advwenov)

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          divx(i,j,k) = 0.0
        enddo
        enddo
        enddo

        call   poiss(uh,vh,mh,rmh,mf,rmf,pi0,thv0,rho0,rf0,    &
                     dum2,dum3,dum4,dum5,divx,pdiag(ibdp,jbdp,kbdp,3),uten,vten,wten,  &
                     cfb,cfa,cfc,ad1,ad2,pdt,lgbth,lgbph,rhs,trans,dt)

        ! Adjust mean pressure ... set average value at upper-most level to zero:

        pavg2 = 0.0

        do j=1,nj
        do i=1,ni
          pavg2 = pavg2+pdiag(i,j,nk,3)
        enddo
        enddo

        pavg2 = pavg2/dble(nx*ny)

        tem = sngl(0.0d0-pavg2)

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          pdiag(i,j,k,3) = pdiag(i,j,k,3)+tem
        enddo
        enddo
        enddo

      !--------------------------------------------------------------!
      ! Coriolis pressure:

      IF( icor.eq.1 )THEN

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj+1
        do i=1,ni+1
          uten(i,j,k) =  0.125*(f2d(i,j)+f2d(i-1,j))           &
                              *( (v3d(i  ,j,k)+v3d(i  ,j+1,k)) &
                                +(v3d(i-1,j,k)+v3d(i-1,j+1,k)) )
          vten(i,j,k) = -0.125*(f2d(i,j)+f2d(i,j-1))           &
                              *( (u3d(i,j  ,k)+u3d(i+1,j  ,k)) &
                                +(u3d(i,j-1,k)+u3d(i+1,j-1,k)) )
          wten(i,j,k) = 0.0
          divx(i,j,k) = 0.0
        enddo
        enddo
        enddo

        call   poiss(uh,vh,mh,rmh,mf,rmf,pi0,thv0,rho0,rf0,    &
                     dum2,dum3,dum4,dum5,divx,pdiag(ibdp,jbdp,kbdp,4),uten,vten,wten,  &
                     cfb,cfa,cfc,ad1,ad2,pdt,lgbth,lgbph,rhs,trans,dt)

        ! Adjust mean pressure ... set average value at upper-most level to zero:

        pavg2 = 0.0

        do j=1,nj
        do i=1,ni
          pavg2 = pavg2+pdiag(i,j,nk,4)
        enddo
        enddo

        pavg2 = pavg2/dble(nx*ny)

        tem = sngl(0.0d0-pavg2)

        !$omp parallel do default(shared)  &
        !$omp private(i,j,k)
        do k=1,nk
        do j=1,nj
        do i=1,ni
          pdiag(i,j,k,4) = pdiag(i,j,k,4)+tem
        enddo
        enddo
        enddo

      ENDIF

      !--------------------------------------------------------------!

      end subroutine pidcomp

  END MODULE pdcomp_module


!     This script was initially intended for
!     making a simple source_sink.in at all elem., and also .th
!
!     ifort -O2 -mcmodel=medium -CB -Bstatic -o gen_elev gen_elev.f90
      implicit real*8(a-h,o-z)
      allocatable :: x(:),y(:),dp(:),area(:),vso(:),r_rough(:),nlayers(:)
      integer, allocatable :: i34(:),elnode(:,:),i_rain(:),isource(:),iDB(:)
      real(8) :: slope(4)

!     Rain rate
      rain_rate= 0.001 !m/hour
      h0=1e-6

      open(10,file='include.gr3',status='old')
      read(10,*); read(10,*)

      open(14,file='hgrid.gr3')
      read(14,*)
      read(14,*)ne,np
      allocate(x(np),y(np),dp(np),r_rough(np),nlayers(np),iDB(np),i_rain(np),isource(ne),i34(ne),elnode(4,ne),area(ne),vso(ne))
      do i=1,np
        read(14,*)j,x(i),y(i),dp(i)
        read(10,*)dummy,dummy,dummy,iDB(i)
      enddo !i
      close(10)
      vol0=0 !inti volume
      do i=1,ne
        read(14,*)j,i34(i),elnode(1:i34(i),i)
        n1=elnode(1,i)
        n2=elnode(2,i)
        n3=elnode(3,i)
        area(i)=signa(x(n1),x(n2),x(n3),y(n1),y(n2),y(n3))
        if(area(i)<=0) then
          print*, 'area<=0 (0):',i,area(i),i34(i)
          stop
        endif

        if(i34(i)==3) then
        else if(i34(i)==4) then
          n4=elnode(4,i)
          area(i)=area(i)+signa(x(n1),x(n3),x(n4),y(n1),y(n3),y(n4))
          if(area(i)<=0) then
            print*, 'area<=0:',i,area(i)
            stop
          endif
        else
          stop 'Unknown elem'
        endif

        tmp=minval(dp(elnode(1:i34(i),i)))
        if(tmp>0) vol0=vol0+area(i)*sum(dp(elnode(1:i34(i),i)))/i34(i)
      enddo !i
      close(14)

      open(15,file='vgrid.in')
      read(15,*) ivcor
      read(15,*) nvrt
      if (ivcor==2) then !s-layer
        nlayers=nvrt-1
      elseif (ivcor==1) then
        do i=1,np
          read(15,*) j,kbp
          nlayers(i)=nvrt-kbp
        enddo
      endif
      close(15)

      total_area=sum(area)
      print*, 'Total area=',total_area !sum(area)
      print*, 'Init volume=',vol0

      open(8,file='rough.gr3',status='replace')
      open(9,file='elev.ic',status='replace')
      write(8,*); write(8,*)ne,np
      write(9,*); write(9,*)ne,np

      r_rough_background=2e-4
      r_rough_land=0.01
      slope_thres=999999.; depth2=-3.0; depth1=-1.0
      !init
      r_rough = r_rough_background

      !based on bathymetry
      do i=1,np
        if (iDB(i) ==0 ) cycle

        if (y(i) > 4452200 ) then !upstream river, y value in utm espg:26918
          r_rough(i)=r_rough_land
        else !DB
          if(dp(i)<depth2) then !upland
            r_rough(i)=r_rough_land
          elseif(dp(i)<=depth1) then !transitional zone
            r_rough(i)=(  (r_rough_land-r_rough_background) * abs(depth1-dp(i))/(depth1-depth2) + r_rough_background )
          endif
        endif ! DE R. or DB
        if (r_rough(i) < -1000) then
          print *, i,dp(i),r_rough(i)
        endif
      enddo

      !multi-layers use small r_rough
      do i=1,ne
        if (maxval(nlayers(elnode(1:i34(i),i))) >1 ) then
          do j=1,i34(i)
            n1=elnode(j,i)
            r_rough(n1) = r_rough_background
          enddo
        endif
      enddo

      !!based on slope and number of layers
      !do i=1,ne
      !  !low land uses small r_rough
      !  if (maxval(dp(elnode(1:i34(i),i))) >= depth_thres) cycle
      !  !multi-layers use small r_rough
      !  if (maxval(nlayers(elnode(1:i34(i),i))) >1 ) cycle

      !  slope = 0.0
      !  do j=1,i34(i)
      !    n1=elnode(j,i);
      !    if ((j+1) <= i34(i)) then
      !      n2=elnode(j+1,i)
      !    else
      !      n2=elnode(1,i)
      !    endif
      !    distance = sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2)
      !    slope(j)=abs((dp(n1)-dp(n2))/distance)
      !  enddo
      !  if (maxval(slope) > slope_thres) then
      !    do j=1,i34(i)-1
      !      n1=elnode(j,i)
      !      r_rough(n1)=-r_CD_manning !negative as drag
      !    enddo
      !  else
      !    do j=1,i34(i)-1
      !      n1=elnode(j,i)
      !      r_rough(n1)=-r_CD_manning !negative as drag
      !      r_rough(n1)=-(  (-r_rough(n1)-drag_min) * maxval(slope)/slope_thres + drag_min ) !linearly varying, slope=0: drag_min; slope=slope_thres: largest drag
      !    enddo
      !  endif
      !enddo

      do i=1,np
        write(8,'(i8,3(1x,f15.6))')i,x(i),y(i),r_rough(i)
        if (iDB(i) ==0 ) then
          write(9,'(i8,3(1x,f15.6))')i,x(i),y(i),max(0.d0,-dp(i))
        else
          write(9,'(i8,3(1x,f15.6))')i,x(i),y(i),max(0.d0,-dp(i)-1e-6)
        endif
      enddo
      do i=1,ne
        write(8,*)i,i34(i),elnode(1:i34(i),i)
        write(9,*)i,i34(i),elnode(1:i34(i),i)
      enddo
      close(8)
      close(9)

      nsource = 0
      do i=1,ne
         !only in selected regions
         if ( i==64871 .or. i==497175 )  then
           nsource = nsource +1
         endif
      enddo !i

      open(10,file='source_sink.in',status='replace')
      write(10,*)nsource
      vso=0.0
      isource =0
      nsource = 0
      do i=1,ne
         !only in selected regions
         if ( i==64871 .or. i==497175 )  then
           write(10,*)i
           nsource = nsource +1
           isource(nsource)=i
           if (i==497175) then
             vso(i)=1500.0
             print *, 'setting Susquehanna'
           elseif (i==64871) then
             vso(i)=400.0
             print *, 'setting Delaware'
           endif
         endif
      enddo !i
      write(10,*)
      write(10,*)0 !# of sinks
      close(10)

      open(10,file='vsource.th',status='replace')
      write(10,'(9000000(1x,e16.5))')0.,vso(isource(1:nsource))
      write(10,'(9000000(1x,e16.5))')8640000.,vso(isource(1:nsource))
      close(10)

      open(12,file='msource.th',status='replace')
      write(12,'(9000000(1x,f3.0))')0.,(24.,j=1,nsource),(0.,j=1,nsource)
      write(12,'(f14.1,1x,9000000(1x,f3.0))')8640000.,(4.,j=1,nsource),(0.,j=1,nsource)
      close(12)

!     Analytical change in volume from init volume (m^3)
!     When comparing to model, remember to remove the init volume in model
      dt=100 !sec
      do it=1,5.*86400/dt
        time=it*dt/86400 !days
        if(time<=0.5) then
          dV=total_area*rain_rate*24*time
        else if(time<=1) then
          dV=total_area*rain_rate*24*(0.75-(time-1)**2)
        else
          dV=total_area*rain_rate*24*0.75
        endif
        write(7,*)time,dV
      enddo !it

      stop
      end

function signa(x1,x2,x3,y1,y2,y3)
!-------------------------------------------------------------------------------
! Compute signed area formed by pts 1,2,3
!-------------------------------------------------------------------------------
  implicit real*8(a-h,o-z)
  real*8,intent(in) :: x1,x2,x3,y1,y2,y3

  signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2d0

end function signa