text:
      SUBROUTINE WBTED( numwbt )
!                  PROCESS WBT  (WATER BOTTOM TIME DEFINITION)
!                  ------- ---
!
!  WRITEUP DATE:  29 April 1993
!
!     PROCESS WBT ASSOCIATES A WATER BOTTOM TIME WITH EVERY TRACE.  THIS TIME IS
!  PLACED IN THE TRACE HEADER AND THEN MAY BE USED BY OTHER PROCESSES SUCH
!  AS MUTE, NMO, OR DECON FOR 'HANGING' WINDOWS FROM. 'HANGING' WINDOWS MEANS
!  THAT THE WINDOW TIMES ARE ADDED TO THE WATER BOTTOM TIME, THUS THE
!  WINDOW TIMES MAY BE THE SAME DISTANCE FROM THE WATER BOTTOM ON ALL TRACES.
!    THE WATER BOTTOM TIMES MAY BE ENTERED IN ONE OF SEVERAL WAYS: 1) BY RP-TIME
!  PAIRS;   2) BY GMT TIME-TIME PAIRS;  3)  CONVERTING DEPTHS TO TIME.
!     THE RP-TIME PAIRS METHOD OF ENTERING THE WATER BOTTOM TIMES CONSIST OF
!  MEARLY ENTERING A LIST OF PAIRS OF RP NUMBERS FOLLOWED BY THE TWO WAY TRAVEL
!  TIME FOR THE RP.  THE RP NUMBERS MUST STRICTLY INCREASE.  RP numbers
!  may be assigned using SIOSEIS processes GEOM or HEADER.
!  THE WATER BOTTOM TIME SPECIFIED WILL BE ASSOCIATED WITH ALL TRACES WITH
!  THE SAME RP NUMBER (A SINGLE CDP HAS THE SAME RP NUMBER!).  WATER BOTTOM
!  TIMES NOT SPECIFIED ARE CALCULATED BY INTERPOLATION OR HELD CONSTANT.
!     THE GMT TIME METHOD CONSISTS OF SPECIFING THE JULIAN DAY VIA THE PARAMETER
!  'DAY' AND THEN A LIST OF GMT TIME - WATER BOTTOM TIME PAIRS.  THE PARAMETER
!  'DAY' MUST BE GIVEN PRIOR TO THE GMT TIME AND MUST BE GIVEN ON THE FIRST
!  PAIR.  FURTHERMORE, THE DAY MUST BE GIVEN ON DAY CHANGES.  THE GMT TIME IS
!  GIVEN IN TERMS OF DECIMAL MINUTES OF THE 24 HOUR CLOCK, THUS 1532.75
!  REPRESENTS 1532 AND 45 SECONDS.  WATER BOTTOM TIMES NOT
!  SPECIFIED ARE CALCULATED BY INTERPOLATION OR ARE HELD CONSTANT.
!  EXAMPLE:
!      WBT
!         41 1.197
!         43 1.199
!      END
!  RPS LESS THAN 41 WILL RECEIVE A WATER BOTTOM TIME OF 1.197, RP 41 A TIME
!  OF 1.197, 42 A TIME OF 1.198, 43 A TIME OF 1.199, AND RPS LARGER THAN 43 A
!  TIME OF 1.199.
!  EXAMPLE 2):
!      WBT
!       DAY 265 1045 .21
!               1100 .25
!               2200 2.0
!       DAY 266 0300 3.5
!      END
!
!    THE DEPTH CONVERSION METHOD CONVERTS THE WATER BOTTOM DEPTH IN WORD 107
!  OF THE TRACE HEADER INTO A TIME USING A VELOCITY OF 1500 METERS PER SECOND.
!
!  THE PARAMETER DICTIONARY
!  --- --------- ----------
!
!  DAY    - The Julian day of the GMT-time pairs that follow.
!  VEL    - The velocity to use in converting the SeaBeam depth to seismic time.
!           Preset = 0.
!  OFFSET - The number of traces between the SeaBeam depth (center of
!            the ship) and the reflection point.
!           Preset = 3
!  NSCAN  - The number of traces to scan for the "closest" SeaBeam depth.
!           only.   NSCAN/2 traces will be scanned in both forward and aft
!           directions.
!           Preset = 10
!
!    *** NOTE ***
!    1)  THE ONLY 'END' NEEDED IS THE 'END' TO TERMINATE THE WBT PARAMETERS.
!
!  COPYRIGHT (C), PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, MARCH 1980
!  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
!  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
!
!
!    mod April 1987 to add the SeaBeam picker.
!    mod 24 July 1990 to make the day/gmt a REAL
!    mod 29 Apr 93 - change the documentation for depth option
!    mod 21 Sep 95 - Add SEL, SES, SOLRAT - an auto picker.
!                  - Remove TYPE
!                  - add PRESTK
!    mod 7 Feb 96 - Adding SEl, SES etc made lists need extra end
!    mod 6 June 97 - Add threshold picking
!    mod Aug 97 - Correctly preset thres
!    mod 22 Oct 98 - Add PEAK
!    mod May 99 - Add SEPP
!    mod June 00 - Add INDEX
!    mod June 01 - Add numwbts and make it re-enterable
!    mod July 01 - Add an edit check for solrat when ses or sel given
!    mod 3 Aug 01 - Add parameter TRACK
!    mod 21 Nov 02 - preset the common to zero FIRST
!                  - Add some error checks
!    mod 9 Sep 05 - Change nscan preset from 11 to 0
!    mod 22 May 07 - Add error if PRESTK is not YES or 1
!    mod 5 Oct 09 - Add bpass
!    mod 10 Jul 14 - Add guided & seg
!
      PARAMETER (npars = 18)
      PARAMETER (MAXWBTS=3)
      REAL ldummy, lrp, lastrp
      DIMENSION ldummy(2)
      EQUIVALENCE (lrp,ldummy(1)), (time,ldummy(2))
      CHARACTER*80 TOKEN
      CHARACTER *6 names(npars)
      CHARACTER *1 types(npars)
      DIMENSION length(npars)
      DIMENSION vals(npars),lvals(npars)
      EQUIVALENCE (vals(1),lvals(1))
      INTEGER day, offset, ns, prestk, prestack
      REAL solrat, selong, seshort, sepeak
      COMMON /EDITS/ IERROR,IWARN,IRUN,NOW,ICOMPT
      COMMON /WBOT/ IWUNIT(MAXWBTS), iday(MAXWBTS), maxwd(MAXWBTS),
     &  lprin(MAXWBTS), loff(MAXWBTS), vela(MAXWBTS), nsca(MAXWBTS),
     &  sol(MAXWBTS), sel(MAXWBTS,2), ses(MAXWBTS,2), prestack(MAXWBTS),
     &  threshold(MAXWBTS), peakp(MAXWBTS), sep(MAXWBTS),
     & sepp(MAXWBTS,2), index1(MAXWBTS), trac(MAXWBTS), pass(MAXWBTS,2),
     &  guide(MAXWBTS), seg(MAXWBTS,2)
      EQUIVALENCE (day,lvals(1)),
     2            (lprint,lvals(3)),
     3            (offset,lvals(4)),
     4            (vel,vals(5)),
     5            (nscan,lvals(6)),
     7            (maxwrd,lvals(7)),
     8            (solrat,vals(8)),
     9            (selong,vals(9)),
     *            (seshort,vals(10)),
     1            (thres,vals(11)),
     2            (peak,vals(12)),
     3            (sepeak,vals(13)),
     4            (index,vals(14)),
     5            (track,vals(15)),
!     6            (ipass,vals(16))
     7            (guided,vals(17)),
     8            (seguide,vals(18))
      DATA names/'DAY   ','PRESTK','LPRINT','OFFSET','VEL   ',
     *           'NSCAN ','MAXWRD','SOLRAT','SEL   ','SES   ',
     *           'THRES ','PEAK  ','SEPP  ','INDEX ','TRACK ',
     *           'PASS  ','GUIDED','SEG   '/
      DATA length/3,6,6,6,3,5,6,6,3,3,5,4,4,5,5,4,6,3/
      DATA types/'L','L','L','L','F','L','L',4*'F','A','F','L','F',
     &           'F','F','F'/
!****
!****       Preset the parameters and GET A FILE FOR THE PARAMETERS
!****
      LASTRP=-32768
      lprint = 0
      day = 0
      lrp = 0
      iend = 0
      prestk = 0
      offset = 3
      vel = 0.
      nscan = 0
      maxwrd = 11000
      solrat = 0.
      selong = 0.
      seshort = 0.
      thres = 0.
      peak = 0.
      sepeak = 0.
      track = 99.
      guided = 0.
      CALL GETFIL(1,IWUNIT(numwbt),TOKEN,ISTAT)
      ns = 0
      ntokes = 1
      index = 50
      iday(numwbt) = 0
      maxwd(numwbt) = 0
      lprin(numwbt) = 0
      loff(numwbt) = 0
      vela(numwbt) = 0.
      nsca(numwbt) = 0
      sol(numwbt) = 0.
      sel(numwbt,1) = 0.
      sel(numwbt,2) = 0.
      ses(numwbt,1) = 0.
      ses(numwbt,2) = 0.
      prestack(numwbt) = 0.
      threshold(numwbt) = 0.
      peakp(numwbt) = 0.
      sep(numwbt) = 0.
      sepp(numwbt,1) = 0.
      sepp(numwbt,2) = 0.
      index1(numwbt) = 0
      trac(numwbt) = 0.
      pass(numwbt,1) = 0.
      pass(numwbt,2) = 0.
      guide(numwbt) = 0.
      seg(numwbt,1) = 0.
      seg(numwbt,2) = 0.
!****
!****      GET A TOKEN
!****
  100 CONTINUE
      CALL GETOKE(TOKEN,NCHARS)
      IF( NCHARS .EQ. 0 ) THEN
          IF( now .EQ. 1 ) PRINT *,' <  ENTER PARAMETERS  >'
          CALL RDLINE
          NTOKES=0
          GO TO 100
      ENDIF
      CALL UPCASE(TOKEN,NCHARS)
  150 ntokes = ntokes + 1
      DO 160 i = 1, npars
         len = length(i)
         nparam = i
         IF( token(1:nchars) .EQ. names(i)(1:len) ) GOTO 200
  160 CONTINUE
      IF( token(1:nchars) .EQ. 'END' ) THEN
          iend = iend + 1
          GOTO 300
      ENDIF
      CALL dcode( token, nchars, value, istat )
      IF( istat .NE. 2 ) THEN
          PRINT *,' ***  ERROR  ***  WBT does not have a parameter ',
     *    'named ',token(1:nchars)
          ierror = ierror + 1
          GOTO 100
      ENDIF
      IF( ns .EQ. 0 ) THEN
          lrp = value
          ns = 1
          GOTO 100
      ELSE
          time = value
          ns = 0
      ENDIF
      IF( day .NE. 0 ) THEN
          IF( day .LT. 0 .OR. day .GT. 366 ) THEN
              PRINT *,' ***  ERROR  ***  Incorrect julian day of ',day
              ierror = ierror + 1
          ENDIF
          IF( lrp .LT. 0 .OR. lrp .GT. 2400 ) THEN
              PRINT *,' ***  ERROR  ***  Incorrect GMT of ',lrp
              ierror = ierror + 1
          ENDIF
          lrp = day * 10000. + lrp
          IF( lrp .LT. lastrp ) THEN
              PRINT *,' ***  ERROR  ***  GMT and DAY must increase.'
              ierror = ierror + 1
          ENDIF
      ELSE
          IF( lrp .LT. lastrp ) THEN
              PRINT *,' ***  ERROR  ***  RP numbers and GMT must ',
     *             'increase.'
              ierror = ierror + 1
          ENDIF
      ENDIF
      CALL wrdisc( iwunit(numwbt), ldummy, 2 )
      IF( IAND(lprint,1) .NE. 0 ) PRINT *,' lrp=',lrp,' time=',time
      GOTO 100
  200 CONTINUE
      ns = 0
  210 CALL getoke( token, nchars )
      IF( nchars .EQ. 0 ) THEN
          IF( now .EQ. 1 ) PRINT *,' <  Enter parameters  >'
          CALL rdline
          ntokes = 0
          GOTO 200
      ENDIF
      CALL upcase( token, nchars )
      ntokes = ntokes + 1
      IF( names(nparam) .EQ. 'PRESTK' ) THEN
          IF( token(1:1) .NE. 'Y' .AND. token(1:1) .NE. '1' ) THEN
              PRINT *,' ***  ERROR  ***  PRESTK must be YES or 1'
              ierror = ierror + 1
          ENDIF
          IF( token(1:1) .EQ. 'Y' ) prestk = 1
          IF( token(1:1) .EQ. '1' ) prestk = 1
          GOTO 100
      ENDIF
      IF( names(nparam) .EQ. 'PEAK' ) THEN
          IF( token(1:3) .EQ. 'POS' ) peak = 1.
          IF( token(1:3) .EQ. 'NEG' ) peak = 2.
          IF( token(1:3) .EQ. 'ABS' ) peak = 3.
          GOTO 100
      ENDIF
      IF( names(nparam) .EQ. 'GUIDED' ) THEN
          IF( token(1:3) .EQ. 'POS' ) guided = 1.
          IF( token(1:3) .EQ. 'NEG' ) guided = 2.
          IF( token(1:3) .EQ. 'ABS' ) guided = 3.
          GOTO 100
      ENDIF
      CALL dcode( token, nchars, value, istat )
      IF( istat .NE. 2 ) THEN
          PRINT *,' ***   ERROR  ***  WBT expected a number instead of',
     *      token(1:nchars)
         ierror = ierror + 1
         GOTO 100
      ENDIF
      IF( types(nparam) .EQ. 'F' ) THEN
          vals(nparam) = value
          IF( names(nparam) .EQ. 'SEL' ) THEN
              ns = ns + 1
              sel(numwbt,ns) = value
              IF( ns .EQ. 1 ) GOTO 210
              ns = 0
              GOTO 100
          ENDIF
          IF( names(nparam) .EQ. 'SES' ) THEN
              ns = ns + 1
              ses(numwbt,ns) = value
              IF( ns .EQ. 1 ) GOTO 210
              ns = 0
              GOTO 100
          ENDIF
          IF( names(nparam) .EQ. 'SEPP' ) THEN
              ns = ns + 1
              sepp(numwbt,ns) = value
              IF( ns .EQ. 1 ) GOTO 210
              ns = 0
              GOTO 100
          ENDIF
          IF( names(nparam) .EQ. 'PASS' ) THEN
              ns = ns + 1
              pass(numwbt,ns) = value
              IF( ns .EQ. 1 ) GOTO 210
              ns = 0
              GOTO 100
          ENDIF
          IF( names(nparam) .EQ. 'SEG' ) THEN
              ns = ns + 1
              seg(numwbt,ns) = value
              IF( ns .EQ. 1 ) GOTO 210
              ns = 0
              GOTO 100
          ENDIF
      ENDIF
      IF( types(nparam) .EQ. 'L' ) lvals(nparam) = value
      GOTO 100
!****
!****     Check the parameters for validity
!****
  300 CONTINUE
      IF( sel(numwbt,2) + ses(numwbt,2).NE. 0 .AND. solrat .EQ. 0 ) THEN
          PRINT *,' ***  ERROR  ***  SOLRAT must be given with ses/sel.'
          ierror = ierror + 1
      ENDIF
      IF( sel(numwbt,1) .GT. sel(numwbt,2) ) THEN
          PRINT *,' ***  ERROR  ***  Bad SELL.'
          ierror = ierror + 1
      ENDIF
      IF( ses(numwbt,1) .GT. ses(numwbt,2) ) THEN
          PRINT *,' ***  ERROR  ***  Bad SES.'
          ierror = ierror + 1
      ENDIF
      IF( sepp(numwbt,1) .GT. sepp(numwbt,2) ) THEN
          PRINT *,' ***  ERROR  ***  Bad SEPP.',
     &        sepp(numwbt,1), sepp(numwbt,2)
          ierror = ierror + 1
      ENDIF
      IF( track .NE. 99. .AND. track .GT. 1 ) THEN
          PRINT *,' ***  WARNING  ***  Unusually large TRACK.'
          PRINT *,' TRACK is in units of time, not depth.'
          iwarn = iwarn + 1
      ENDIF
      IF( guided .NE. 0 .AND. seg(numwbt,1) .EQ. 0 .AND.
     &    seg(numwbt,2) .EQ. 0 ) THEN
      PRINT *,' ***  ERROR  ***  SEG must be given when GUIDED is used'
          ierror = ierror + 1
      ENDIF
      iday(numwbt) = day
      nsca(numwbt) = nscan
      vela(numwbt) = vel
      loff(numwbt) = offset
      maxwd(numwbt) = maxwrd
      lprin(numwbt) = lprint
      sol(numwbt) = solrat
      prestack(numwbt) = prestk
      threshold(numwbt) = thres
      peakp(numwbt) = peak
      sep(numwbt) = sepeak
      index1(numwbt) = index
      trac(numwbt) = track
      guide(numwbt) = guided
      IF( IAND(lprint,1) .NE. 0 ) THEN
          PRINT *,' numwbt=',numwbt,' off=',loff,' vel=',vela,
     &            ' nscan=',nscan
          PRINT *,' solrat=',solrat,' sel=',sel,' ses=',ses
          PRINT *,' prestack=',prestack,' thres=',threshold
          PRINT *,' peak=',peakp,' index=',index1,' track=',track
          PRINT *,' guided=',guide,' seg=',seg
      ENDIF
      IF( ns .EQ. 1 ) THEN
          PRINT *,' ***  ERROR  ***  WBT pairs pairs are not pairs!'
          ierror = ierror + 1
      ENDIF
 2000 IF( lrp .NE. 0 ) RETURN
! 2000 IF( solrat+selong+seshort+thres+peak .EQ. 0. ) RETURN
      CALL getoke( token, nchars )
      IF( nchars .LE. 0 ) THEN
          IF( now .EQ. 1 ) PRINT *,' <  ENTER PARAMETERS  >'
          CALL rdline
           ! get a new line of parameters
          ntokes = 0
          GOTO 2000
      ENDIF
      CALL upcase( token, nchars )
      ntokes = ntokes + 1
      IF(  token(1:nchars) .NE. 'END' .OR. nchars .NE. 3 ) THEN
          PRINT *,' ***  ERROR  ***  WBT allows only one fno/lno list.'
          ierror = ierror + 1
      ENDIF
      RETURN

      END


text:
module lattice_mod
  use mod_d2q9_type, only: d2q9_lattice_t
  use mod_kinds, only: ik, rk
  implicit none

  class(d2q9_lattice_t), allocatable :: lattice

contains

  subroutine initialize()
    allocate(d2q9_lattice_t :: lattice)
    call lattice%initialize()
  end subroutine initialize

  subroutine test_basic()
    !< Test basic settings like weights, dimension, etc.

    if(lattice%n_weights /= 9) then
      if(this_image() == 1) write(*, '(a,i0)') "Error: lattice%n_weights should be 9, but is: ", lattice%n_weights
      error stop
    endif

    if(lattice%dimension /= 2) then
      if(this_image() == 1) write(*, '(a,i0)') "Error: lattice%dimension should be 2, but is: ", lattice%dimension
      error stop
    endif
  end subroutine test_basic

  subroutine test_neighbors()
    !< Test the lattice neighbor values, so that the proper directional velocities are used in streaming

    integer :: i, j

    do i = 1, lattice%n_weights - 1
      associate(neighbor=>lattice%neighbor_indices(i))
        associate(node=>int(lattice%v_hat(i, :)), &
                  other_node=>int(lattice%v_hat(neighbor, :)))

          do j = 1, 2
            if(node(j) /= 0 .or. other_node(j) /= 0) then
              if(node(j) /= -other_node(j)) then
                if(this_image() == 1) then
                  write(*, *) "Error: lattice%neighbor_indices is wrong: [", node, "] -> neighbor: [", other_node, "]"
                end if
                error stop
              end if
            end if
          end do

        end associate
      end associate
    end do

  end subroutine test_neighbors

end module lattice_mod

program test_lattices

  use lattice_mod
  implicit none

  call initialize()
  call test_basic()
  call test_neighbors()

end program test_lattices

max_stars_repo_name: smillerc/lbm

max_issues_repo_name: smillerc/lbm

max_forks_repo_name: smillerc/lbm


text:
      DOUBLE PRECISION FUNCTION slGMSA (DATE, UT)
*+
*     - - - - - -
*      G M S A
*     - - - - - -
*
*  Conversion from Universal Time to Greenwich mean sidereal time,
*  with rounding errors minimized.
*
*  double precision
*
*  Given:
*    DATE    d      UT1 date (MJD: integer part of JD-2400000.5))
*    UT      d      UT1 time (fraction of a day)
*
*  The result is the Greenwich mean sidereal time (double precision,
*  radians, in the range 0 to 2pi).
*
*  There is no restriction on how the UT is apportioned between the
*  DATE and UT arguments.  Either of the two arguments could, for
*  example, be zero and the entire date+time supplied in the other.
*  However, the routine is designed to deliver maximum accuracy when
*  the DATE argument is a whole number and the UT lies in the range
*  0 to 1 (or vice versa).
*
*  The algorithm is based on the IAU 1982 expression (see page S15 of
*  the 1984 Astronomical Almanac).  This is always described as giving
*  the GMST at 0 hours UT1.  In fact, it gives the difference between
*  the GMST and the UT, the steady 4-minutes-per-day drawing-ahead of
*  ST with respect to UT.  When whole days are ignored, the expression
*  happens to equal the GMST at 0 hours UT1 each day.  Note that the
*  factor 1.0027379... does not appear explicitly but in the form of
*  the coefficient 8640184.812866, which is 86400x36525x0.0027379...
*
*  In this routine, the entire UT1 (the sum of the two arguments DATE
*  and UT) is used directly as the argument for the standard formula.
*  The UT1 is then added, but omitting whole days to conserve accuracy.
*
*  See also the routine slGMST, which accepts the UT as a single
*  argument.  Compared with slGMST, the extra numerical precision
*  delivered by the present routine is unlikely to be important in
*  an absolute sense, but may be useful when critically comparing
*  algorithms and in applications where two sidereal times close
*  together are differenced.
*
*  Called:  slDA2P
*
*  P.T.Wallace   Starlink   14 October 2001
*
*  Copyright (C) 2001 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION DATE,UT

*  Seconds of time to radians
      DOUBLE PRECISION S2R
      PARAMETER (S2R=7.272205216643039903848712D-5)

      DOUBLE PRECISION D1,D2,T
      DOUBLE PRECISION slDA2P

*  Julian centuries since J2000.
      IF (DATE.LT.UT) THEN
         D1=DATE
         D2=UT
      ELSE
         D1=UT
         D2=DATE
      END IF
      T=(D1+(D2-51544.5D0))/36525D0

*  GMST at this UT1.
      slGMSA=slDA2P(S2R*(24110.54841D0+
     :                         (8640184.812866D0+
     :                         (0.093104D0
     :                         -6.2D-6*T)*T)*T
     :                         +86400D0*(MOD(D1,1D0)+MOD(D2,1D0))))

      END

ext: f

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


!=============================================================================!
        subroutine impRK()
!
!  Driver for implicit R-K time advancement of the compressible Navier-Stokes
!  solver.  This routine could also be used as the backbone of a MF-GMRES
!  driver.  However, the RK scheme doesn't really work well enough to justify
!  further work right now.  The R-K scheme is the scheme proposed by Jameson
!  that is supposed to be quite dissipative.  The coefficients come from
!  Charles Fenno's thesis from NC State [1993].
!
!  Author:   Scott Collis
!
!  Written:  6-5-95
!
!=============================================================================!
        use stuff
        use local
        implicit none
!=============================================================================!
        real, allocatable :: v(:), vold(:), v2old(:), r(:)
        real, allocatable :: vint(:)
#ifdef __GFORTRAN__
  external :: fstat
#endif
        real :: rtime, cpu, dtau
        integer :: i, j, idof
        integer :: ier, mem = 0
        integer, external :: igetver
        real, external :: second
!=============================================================================!

!.... allocate the storage area for the disturbance field

        mem = mem + ny*nx*ndof
        allocate (v(ny*nx*ndof), vold(ny*nx*ndof), STAT=ier)
        if (ier .ne. 0) call error('imprk$','Insufficient Memory for v$')

!.... generate the initial condition

        call genini(v, vold)

!.... allocate memory for the RHS

        mem = mem + ny*nx*ndof
        allocate (r(ny*nx*ndof), STAT=ier)
        if (ier .ne. 0) call error('impRK$','Insufficient Memory for r$')

!.... allocate memory for the temporal v

        mem = mem + ny*nx*ndof
        allocate (vint(ny*nx*ndof), STAT=ier)
        if (ier .ne. 0) call error('impRK$','Insufficient Memory for vint$')

!.... allocate memory for the previous, previous disturbance field

        if (impl.eq.2) then
          mem = mem + ny*nx*ndof
          allocate (v2old(ny*nx*ndof), STAT=ier)
          if (ier .ne. 0) call error('impRK$', &
                                     'Insufficient Memory for v2old$')
        end if

        write(*,"(' ImpRK allocated ===> ',1pe13.6,' words')") float(mem)
!=============================================================================!

!.... dtau is the time step in pseudotime.  I don't have a rational way
!.... setting it right now.

        dtau = 0.5

!.... open the history file

        iver = igetver( base, 'h'//char(0) ) + 1
        call addver(base, 'h'//char(0), iver, filen, lfile)
        open(unit=ihist, file=filen, form='formatted', status='unknown')

!.... compute the CFL

        call dtcfl(vm)

!.... loop through the time steps

        do istep = 1, nstep

!.... -----------------------> predictor phase <-----------------------

          if (lstep .gt. 0) v2old = vold
          vold  = v
          time  = time + Delt
          call itrbc(v,vm)

!.... --------------------> multi-corrector phase <--------------------

        do iter = 1, niter

!.... Two-step backward

            if ( lstep .gt. 0 ) then

              call lrhs(v, r, vm, x, y)
              r = v - onept33 * vold + pt33 * v2old +           &
                  pt66 * delt * r
              call rkbc(r)
              call resstat(r)

              vint = v - 0.25 * dtau * r
              call itrbc(vint,vm)
              call lrhs(vint, r, vm, x, y)
              r = vint - onept33 * vold + pt33 * v2old +        &
                  pt66 * delt * r

              vint = v - 0.50 * dtau * r
              call itrbc(vint,vm)
              call lrhs(vint, r, vm, x, y)
              r = vint - onept33 * vold + pt33 * v2old +        &
                  pt66 * delt * r

              vint = v - 0.55 * dtau * r
              call itrbc(vint,vm)
              call lrhs(vint, r, vm, x, y)
              r = vint - onept33 * vold + pt33 * v2old +        &
                  pt66 * delt * r

            else

!.... Backward Euler

              call lrhs( v, r, vm, x, y)
              r = v - vold + delt * r
              call rkbc(r)
              call resstat(r)

              vint = v - 0.25 * dtau * r
              call itrbc(vint,vm)
              call lrhs( vint, r, vm, x, y)
              r = vint - vold + delt * r

              vint = v - 0.50 * dtau * r
              call itrbc(vint,vm)
              call lrhs( vint, r, vm, x, y)
              r = vint - vold + delt * r

              vint = v - 0.55 * dtau * r
              call itrbc(vint,vm)
              call lrhs( vint, r, vm, x, y)
              r = vint - vold + delt * r

            end if

!.... update v

            v = v - dtau * r

!.... satisfy boundary conditions on v

            call itrbc(v,vm)

          end do                ! loop on iter

          lstep = lstep + 1

!.... quit if less than 10 seconds remain

          call tremain(rtime)
          if (rtime .le. 10) goto 5000

!.... compute temporal statistics

          if(.true. .and. xper .and. istep .ne. nstep .and. &
             mod(istep-1,itout) .eq. 0 ) then
            call fstat(.false., v, vm)          ! don't print profiles
          end if

!.... compute spatial statistics

          if (.false. .and. linear.eq.1) then
            call sstat(v, vm, x, y)
          end if

!.... write out restart and q files

          if ( mod(lstep,ntout) .eq. 0 .and. istep .ne. nstep) then

            iver = igetver( base, 'R'//char(0) ) + 1
            call addver(base, 'R'//char(0), iver, filen, lfile)
            open(unit=ifile, file=filen, form='unformatted', &
                  status='unknown')
            write(ifile) lstep, time, nx, ny, nz, ndof
            write(ifile) v
            if (impl.eq.2) write(ifile) vold
            close(ifile)

            write(*,10) filen(1:index(filen,char(0))-1), lstep, time
            write(ihist,10) filen(1:index(filen,char(0))-1), lstep, time

          end if

        end do          ! loop on istep

!.... ---------------------->  Post Processing  <----------------------

5000    continue

!.... compute final statistics

        if (.true. .and. xper) then
          call fstat(.true., v, vm)     ! print profiles
        end if

!.... write the restart file

        iver = igetver( base, 'R'//char(0) ) + 1
        call addver(base, 'R'//char(0), iver, filen, lfile)
        open(unit=ifile, file=filen, form='unformatted', &
             status='unknown')
        write(ifile) lstep, time, nx, ny, nz, ndof
        write(ifile) v
        if (impl.eq.2) write(ifile) vold
        close(ifile)

        write(*,10) filen(1:index(filen,char(0))-1), lstep, time
        write(ihist,10) filen(1:index(filen,char(0))-1), lstep, time

!.... close the history file

        close(ihist)

        return

  10    format(' Wrote ',a,' at step = ',i7,' time = ',1pe13.6)

        end


text:
	PROGRAM TESTTB
C************************************************************************
C* TESTTB								*
C*									*
C* This program tests the TABLE library subroutines			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 5/84						*
C* M. desJardins/GSFC	10/88	Rewritten				*
C* S. Jacobs/EAI	 7/92	Added FL_TOPN				*
C* K. Brill/NMC		 8/93	stnid*4 -> stnid*8, ISPRI		*
C* D. Keiser/GSC	12/95	Added FL_TBOP and removed FL_TOPN	*
C* G. Krueger/EAI	 5/96	Crnr pts, default proj; Check bad table	*
C* L. Sager/NCEP         6/96   Changed calling sequence for TB_RSTN    *
C*                              to add additional parameter string      *
C* L. Sager/NCEP	 6/96   Added TB_ASTN				*
C* S. Schotz/NCEP	 7/96	Added option to write out all stations  *
C*				for TB_ASTN, added FL_CLOS		*
C* S. Jacobs/NCEP	11/96	Added station name to TB_ASTN; Clean up	*
C* S. Maxwell/GSC	 2/97	Added TB_PARM				*
C* S. Jacobs/NCEP	12/97	Added TB_NIDS				*
C* S. Jacobs/NCEP	 7/98	Added TB_FONT				*
C* A. Hardy/GSC		 8/99   Added TB_IDST				*
C* S. Jacobs/NCEP	12/99	Changed size of tbchr and tbchrs 14->20	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER 	file*72, geog*20, cdproj*30, stnid*8, stnnam*32,
     +			state*2, contry*2, type*72, tbchrs*20,
     +			id (LLSTFL)*8, tbchr (LLSTFL)*20,
     +			st (LLSTFL)*2, con (LLSTFL)*2,
     +			stnam (LLSTFL)*32, prmlst*72, colors*72,
     +			parms*72, prdnam*8, units*8, desc*20, font*12
	REAL      	rla (LLSTFL), rlo (LLSTFL), hgh (LLSTFL)
	INTEGER		inum (LLSTFL), ispi (LLSTFL)
C------------------------------------------------------------------------
	CALL IN_BDTA ( iret )
	iostat = 0
	DO WHILE ( iostat .eq. 0 )
	    WRITE(6,20)
20	    FORMAT
     +  ('  1 = TB_FGEO    2 = TB_RSTN    3 = TB_ASTN '/
     +   '  4 = TB_PARM    5 = TB_NIDS    6 = TB_FONT '/
     +   '  7 = TB_IDST '/
     +   ' 20 = FL_SWOP   21 = FL_APND   22 = FL_TBOP '/
     +   ' 23 = FL_CLOS   24 = FL_REWD'/ )
C
	    CALL TM_INT ( 'Select a subroutine number', .false.,
     +      .false., 1, numsub, n, ier )
	IF ( ier .eq. 2 ) THEN
	   iostat = -1
           numsub = -1
 	END IF
2	    FORMAT (A)
C------------------------------------------------------------------------
	    IF (numsub .eq. 1) THEN
		WRITE (6,*) 'Enter GEOG'
		READ  (5,2)  geog
		CALL  TB_FGEO  ( geog, rlatll, rlonll, rlatur, rlonur,
     +				 cdproj, cenlat, cenlon, iret )
		WRITE (6,*) 'RLATLL,RLONLL,RLATUR,RLONUR:',
     +			    rlatll,rlonll,rlatur,rlonur
		WRITE (6,*) 'CDPROJ:', cdproj
		WRITE (6,4) 'CENLAT,CENLON,IRET:', cenlat,cenlon,iret
4	    	FORMAT (A, 2F8.2, I4)
	        IF ( (iret .ne. 0) .and.
     +		     (iret .ne. 1) .and. (iret .ne. -10) ) THEN
		    CALL ER_WMSG ( 'TB', iret, geog, ier )
		END IF
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 2) THEN
	        WRITE (6,*) 'ENTER LUN'
	        READ  (5,*)  lun
		CALL TB_RSTN  ( lun, stnid, stnnam, isnum, state,
     +				contry, rlat, rlon, hght, ispri,
     +				tbchrs, iret )
		WRITE (6,*) 'IRET = ', IRET
	        WRITE (6,1000)  stnid, isnum, stnnam, state, contry,
     +				rlat, rlon, hght, ispri, tbchrs
1000		FORMAT ( A, 1X, I6, 1X, A, 1X, A, 1X, A, 1X,
     +                   F9.2, 1X, F9.2, 1X, F9.2, 1X, I2, 1X, A )
	        IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
              ELSE IF (numsub .eq. 3) THEN
                WRITE (6,*) 'ENTER LUN'
                READ  (5,*)  lun
                WRITE (6,*) 'ENTER MAXSTN, the number of stns to read'
                READ  (5,*)  maxstn
                CALL TB_ASTN  ( lun, maxstn, nstn, id, stnam, inum, st,
     +				con, rla, rlo, hgh, ispi, tbchr, iret )
		WRITE (6,*) 'Number of stations read, NSTN = ', nstn
                WRITE (6,*) 'IRET = ', iret
		IF (nstn .gt. 0 ) THEN
    		    WRITE (6,*) ' To list all the stations read,',
     +				' enter 1:'
		    READ  (5,*) lststn
		    IF ( lststn .eq. 1 ) THEN
		        DO i = 1, nstn
		           WRITE (6,1000) id (i), inum (i), stnam (i),
     +					  st (i), con (i), rla (i),
     +					  rlo (i), hgh (i), ispi (i),
     +					  tbchr (i)
	      	        END DO
		    END IF
		END IF
                IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 4)  THEN
		WRITE (6,*) ' ENTER ALIAS NAME:'
		READ  (5,2) prmlst
		CALL TB_PARM ( prmlst, parms, colors, iret )
		WRITE (6,*) 'PARMS = ',parms
		WRITE (6,*) 'COLORS = ',colors
		WRITE (6,*) 'IRET = ', iret
                IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 5)  THEN
		WRITE (6,*) 'Enter source ID number:'
		READ  (5,*) isrc
		WRITE (6,*) 'Enter product ID number:'
		READ  (5,*) iprod
		CALL TB_NIDS ( isrc, iprod, prdnam, nlev, units,
     +			       res, desc, iret )
		IF  ( iret .eq. 0 )  THEN
		    CALL ST_RNUL(prdnam, prdnam, lens, ier )
		    CALL ST_RNUL(units, units, lens, ier )
		    CALL ST_RNUL(desc, desc, lens, ier )
		    WRITE (6,*) 'PRDNAM = ', prdnam
		    WRITE (6,*) 'NLEV   = ', nlev
		    WRITE (6,*) 'UNITS  = ', units
		    WRITE (6,*) 'RES    = ', res
		    WRITE (6,*) 'DESC   = ', desc
		END IF
		WRITE (6,*) 'IRET   = ', iret
                IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 6)  THEN
		WRITE (6,*) ' Enter font size name:'
		READ  (5,2) font
		CALL TB_FONT ( font, fsize, iret )
		WRITE (6,*) 'Font size = ', fsize
		WRITE (6,*) 'IRET   = ', iret
                IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 7) THEN
		WRITE (6,*) 'Enter the 2 character state ID'
                READ(5,2)state
		CALL TB_IDST  ( state, stnnam, iret )
		WRITE (6,*) 'IRET = ', IRET
		WRITE (6,*) 'State is : ',stnnam
	        IF (iret .ne. 0) CALL ER_WMSG ( 'TB', iret, ' ', ier )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 20) THEN
		WRITE (6,*) 'ENTER FILE'
		READ  (5,2)  file
		CALL FL_SWOP  ( file, lun, iret )
		WRITE (6,*) 'LUN,IRET: ', lun, iret
	        IF (iret .ne.. 0) CALL ER_WMSG ('FL', iret, ' ', ier)
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 21 )  THEN
		WRITE (6,*) 'ENTER LUN'
		READ  (5,*)  lun
		CALL FL_APND  ( lun, iret )
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 22)  THEN
		WRITE (6,*)' ENTER FILE'
		READ  (5,2) file
		WRITE (6,*)' ENTER TYPE'
		READ  (5,2) type
		CALL FL_TBOP ( file, type, lun, iret )
		WRITE (6,*) 'lun, iret ', lun, iret
		IF (iret .ne. 0) CALL ER_WMSG ( 'FL',  iret,
     +                                          file, ierr)
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 23)  THEN
		WRITE (6,*) ' ENTER LUN'
		READ (5,*) lun
		CALL FL_CLOS (lun, iret)
C------------------------------------------------------------------------
	      ELSE IF (numsub .eq. 24)  THEN
		WRITE (6,*) ' ENTER LUN'
		READ (5,*) lun
		CALL FL_REWD (lun, iret)
C------------------------------------------------------------------------
	    END IF
	END DO
C
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2015 NCAR/RAL
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.

module check_icond_module
USE nrtype

! access missing values
USE globalData,only:integerMissing  ! missing integer
USE globalData,only:realMissing     ! missing double precision number

! define modeling decisions
USE mDecisions_module,only:  &
 moisture,                   & ! moisture-based form of Richards' equation
 mixdform                      ! mixed form of Richards' equation

implicit none
private
public::check_icond
contains

 ! ************************************************************************************************
 ! public subroutine check_icond: read model initial conditions
 ! ************************************************************************************************
 subroutine check_icond(nGRU,                          & ! number of GRUs and HRUs
                        progData,                      & ! model prognostic (state) variables
                        mparData,                      & ! model parameters
                        indxData,                      & ! layer index data
                        err,message)                     ! error control
 ! --------------------------------------------------------------------------------------------------------
 ! modules
 USE nrtype
 USE var_lookup,only:iLookParam                          ! variable lookup structure
 USE var_lookup,only:iLookProg                           ! variable lookup structure
 USE var_lookup,only:iLookIndex                          ! variable lookup structure
 USE globalData,only:gru_struc                           ! gru-hru mapping structures
 USE data_types,only:gru_hru_doubleVec                   ! actual data
 USE data_types,only:gru_hru_intVec                      ! actual data
 USE globaldata,only:iname_soil,iname_snow               ! named variables to describe the type of layer
 USE multiconst,only:&
                       LH_fus,    &                      ! latent heat of fusion                (J kg-1)
                       iden_ice,  &                      ! intrinsic density of ice             (kg m-3)
                       iden_water,&                      ! intrinsic density of liquid water    (kg m-3)
                       gravity,   &                      ! gravitational acceleration           (m s-2)
                       Tfreeze                           ! freezing point of pure water         (K)
 USE snow_utils_module,only:fracliquid                   ! compute volumetric fraction of liquid water in snow based on temperature
 USE updatState_module,only:updateSnow                   ! update snow states
 USE updatState_module,only:updateSoil                   ! update soil states
 implicit none

 ! --------------------------------------------------------------------------------------------------------
 ! variable declarations
 ! dummies
 integer(i4b)           ,intent(in)    :: nGRU           ! number of grouped response units
 type(gru_hru_doubleVec),intent(inout) :: progData       ! prognostic vars
 type(gru_hru_doubleVec),intent(in)    :: mparData       ! parameters
 type(gru_hru_intVec)   ,intent(in)    :: indxData       ! layer indexes
 integer(i4b)           ,intent(out)   :: err            ! error code
 character(*)           ,intent(out)   :: message        ! returned error message

 ! locals
 character(len=256)                     :: cmessage      ! downstream error message
 integer(i4b)                           :: iGRU          ! loop index
 integer(i4b)                           :: iHRU          ! loop index

 ! temporary variables for realism checks
 integer(i4b)                      :: iLayer             ! index of model layer
 integer(i4b)                      :: iSoil              ! index of soil layer
 real(dp)                          :: fLiq               ! fraction of liquid water on the vegetation canopy (-)
 real(dp)                          :: vGn_m              ! van Genutchen "m" parameter (-)
 real(dp)                          :: tWat               ! total water on the vegetation canopy (kg m-2)
 real(dp)                          :: scalarTheta        ! liquid water equivalent of total water [liquid water + ice] (-)
 real(dp)                          :: h1,h2              ! used to check depth and height are consistent
 integer(i4b)                      :: nLayers            ! total number of layers
 real(dp)                          :: kappa              ! constant in the freezing curve function (m K-1)
 integer(i4b)                      :: nSnow              ! number of snow layers

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

 ! Start procedure here
 err=0; message="check_icond/"

 ! --------------------------------------------------------------------------------------------------------
 ! Check that the initial conditions do not conflict with parameters, structure, etc.
 ! --------------------------------------------------------------------------------------------------------
 do iGRU = 1,nGRU
  do iHRU = 1,gru_struc(iGRU)%hruCount
   ! ensure the spectral average albedo is realistic
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMax)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMax)%dat(1)
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinWinter)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinWinter)%dat(1)
   ! ensure the visible albedo is realistic
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxVisible)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxVisible)%dat(1)
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinVisible)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinVisible)%dat(1)
   ! ensure the nearIR albedo is realistic
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxNearIR)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxNearIR)%dat(1)
   if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinNearIR)%dat(1)) &
      progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinNearIR)%dat(1)
  end do
 end do

 ! ensure the initial conditions are consistent with the constitutive functions
 do iGRU = 1,nGRU
  do iHRU = 1,gru_struc(iGRU)%hruCount

   ! associate local variables with variables in the data structures
   associate(&
   ! state variables in the vegetation canopy
   scalarCanopyTemp  => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyTemp)%dat(1)   , & ! canopy temperature
   scalarCanopyIce   => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyIce)%dat(1)    , & ! mass of ice on the vegetation canopy (kg m-2)
   scalarCanopyLiq   => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyLiq)%dat(1)    , & ! mass of liquid water on the vegetation canopy (kg m-2)
   ! state variables in the snow+soil domain
   mLayerTemp        => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerTemp)%dat            , & ! temperature (K)
   mLayerVolFracLiq  => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracLiq)%dat      , & ! volumetric fraction of liquid water in each snow layer (-)
   mLayerVolFracIce  => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracIce)%dat      , & ! volumetric fraction of ice in each snow layer (-)
   mLayerMatricHead  => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerMatricHead)%dat      , & ! matric head (m)
   mLayerLayerType   => indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%layerType)%dat            , & ! type of layer (ix_soil or ix_snow)
   ! depth varying soil properties
   vGn_alpha         => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%vGn_alpha)%dat            , & ! van Genutchen "alpha" parameter (m-1)
   vGn_n             => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%vGn_n)%dat                , & ! van Genutchen "n" parameter (-)
   theta_sat         => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%theta_sat)%dat            , & ! soil porosity (-)
   theta_res         => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%theta_res)%dat            , & ! soil residual volumetric water content (-)
   ! snow parameters
   snowfrz_scale     => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%snowfrz_scale)%dat(1)     , & ! scaling parameter for the snow freezing curve (K-1)
   FCapil            => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%FCapil)%dat(1)              & ! fraction of pore space in tension storage (-)
   )  ! (associate local variables with model parameters)

   ! compute the constant in the freezing curve function (m K-1)
   kappa  = (iden_ice/iden_water)*(LH_fus/(gravity*Tfreeze))  ! NOTE: J = kg m2 s-2

   ! modify the liquid water and ice in the canopy
   if(scalarCanopyIce > 0._dp .and. scalarCanopyTemp > Tfreeze)then
    message=trim(message)//'canopy ice > 0 when canopy temperature > Tfreeze'
    err=20; return
   end if
   fLiq = fracliquid(scalarCanopyTemp,snowfrz_scale)  ! fraction of liquid water (-)
   tWat = scalarCanopyLiq + scalarCanopyIce           ! total water (kg m-2)
   scalarCanopyLiq = fLiq*tWat                        ! mass of liquid water on the canopy (kg m-2)
   scalarCanopyIce = (1._dp - fLiq)*tWat              ! mass of ice on the canopy (kg m-2)

   ! number of layers
   nLayers = gru_struc(iGRU)%hruInfo(iHRU)%nSnow + gru_struc(iGRU)%hruInfo(iHRU)%nSoil
   nSnow   = gru_struc(iGRU)%hruInfo(iHRU)%nSnow

   ! loop through all layers
   do iLayer=1,nLayers

    ! compute liquid water equivalent of total water (liquid plus ice)
    if (iLayer>nSnow) then ! soil layer = no volume expansion
     iSoil       = iLayer - nSnow
     vGn_m       = 1._dp - 1._dp/vGn_n(iSoil)
     scalarTheta = mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)
    else ! snow layer = volume expansion allowed
     iSoil       = integerMissing
     vGn_m       = realMissing
     scalarTheta = mLayerVolFracIce(iLayer)*(iden_ice/iden_water) + mLayerVolFracLiq(iLayer)
    end if

    ! *****
    ! * check that the initial volumetric fraction of liquid water and ice is reasonable...
    ! *************************************************************************************
    select case(mlayerLayerType(iLayer))

     ! ***** snow
     case(iname_snow)
      ! (check liquid water)
      if(mLayerVolFracLiq(iLayer) < 0._dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < 0: layer = ',iLayer; err=20; return; end if
      if(mLayerVolFracLiq(iLayer) > 1._dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > 1: layer = ',iLayer; err=20; return; end if
      ! (check ice)
      if(mLayerVolFracIce(iLayer) > 0.80_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > 0.80: layer = ',iLayer; err=20; return; end if
      if(mLayerVolFracIce(iLayer) < 0.05_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0.05: layer = ',iLayer; err=20; return; end if
      ! check total water
      if(scalarTheta > 0.80_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > 0.80: layer = ',iLayer; err=20; return; end if
      if(scalarTheta < 0.05_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < 0.05: layer = ',iLayer; err=20; return; end if

     ! ***** soil
     case(iname_soil)
      ! (check liquid water)
      if(mLayerVolFracLiq(iLayer) < theta_res(iSoil) )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < theta_res: layer = ',iLayer; err=20; return; end if
      if(mLayerVolFracLiq(iLayer) > theta_sat(iSoil) )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > theta_sat: layer = ',iLayer; err=20; return; end if
      ! (check ice)
      if(mLayerVolFracIce(iLayer) < 0._dp            )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0: layer = '        ,iLayer; err=20; return; end if
      if(mLayerVolFracIce(iLayer) > theta_sat(iSoil) )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > theta_sat: layer = ',iLayer; err=20; return; end if
      ! check total water
      if(scalarTheta < theta_res(iSoil) )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < theta_res: layer = ',iLayer; err=20; return; end if
      if(scalarTheta > theta_sat(iSoil) )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > theta_sat: layer = ',iLayer; err=20; return; end if

     case default
      err=20; message=trim(message)//'cannot identify layer type'; return
    end select

    ! *****
    ! * check that the initial conditions are consistent with the constitutive functions...
    ! *************************************************************************************
    select case(mLayerLayerType(iLayer))

     ! ** snow
     case(iname_snow)

      ! check that snow temperature is less than freezing
      if(mLayerTemp(iLayer) > Tfreeze)then
       message=trim(message)//'initial snow temperature is greater than freezing'
       err=20; return
      end if

      ! ensure consistency among state variables
      call updateSnow(&
                      ! input
                      mLayerTemp(iLayer),                        & ! intent(in): temperature (K)
                      scalarTheta,                               & ! intent(in): mass fraction of total water (-)
                      snowfrz_scale,                             & ! intent(in): scaling parameter for the snow freezing curve (K-1)
                      ! output
                      mLayerVolFracLiq(iLayer),                  & ! intent(out): volumetric fraction of liquid water (-)
                      mLayerVolFracIce(iLayer),                  & ! intent(out): volumetric fraction of ice (-)
                      fLiq,                                      & ! intent(out): fraction of liquid water (-)
                      err,cmessage)                                ! intent(out): error control
      if(err/=0)then; message=trim(message)//trim(cmessage); return; end if  ! (check for errors)

     ! ** soil
     case(iname_soil)

      ! ensure consistency among state variables
      call updateSoil(&
                      ! input
                      mLayerTemp(iLayer),                                                    & ! intent(in): layer temperature (K)
                      mLayerMatricHead(iLayer-nSnow),                                        & ! intent(in): matric head (m)
                      vGn_alpha(iSoil),vGn_n(iSoil),theta_sat(iSoil),theta_res(iSoil),vGn_m, & ! intent(in): van Genutchen soil parameters
                      ! output
                      scalarTheta,                                                           & ! intent(out): volumetric fraction of total water (-)
                      mLayerVolFracLiq(iLayer),                                              & ! intent(out): volumetric fraction of liquid water (-)
                      mLayerVolFracIce(iLayer),                                              & ! intent(out): volumetric fraction of ice (-)
                      err,cmessage)                                                            ! intent(out): error control
      if(err/=0)then; message=trim(message)//trim(cmessage); return; end if  ! (check for errors)

     case default; err=10; message=trim(message)//'unknown case for model layer'; return
    end select

   end do  ! (looping through layers)

   ! end association to variables in the data structures
   end associate

   ! if snow layers exist, compute snow depth and SWE
   if(nSnow > 0)then
    progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSWE)%dat(1)       = sum( (progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
                                                                               progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice)        * &
                                                                               progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
   end if  ! if snow layers exist

   ! check that the layering is consistent
   do iLayer=1,nLayers
    h1 = sum(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:iLayer)) ! sum of the depths up to the current layer
    h2 = progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(iLayer) - progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(0)  ! difference between snow-atm interface and bottom of layer
    if(abs(h1 - h2) > 1.e-12_dp)then
     write(message,'(a,1x,i0)') trim(message)//'mis-match between layer depth and layer height [suggest round numbers in initial conditions file]; layer = ', iLayer
     err=20; return
    end if
   end do

  end do ! iHRU
 end do ! iGRU

 end subroutine check_icond

end module check_icond_module

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


!*==slaebz.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download SLAEBZ + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
!                          RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
!                          NAB, WORK, IWORK, INFO )
!
!       .. Scalar Arguments ..
!       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
!       REAL               ABSTOL, PIVMIN, RELTOL
!       ..
!       .. Array Arguments ..
!       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
!       REAL               AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
!      $                   WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> SLAEBZ contains the iteration loops which compute and use the
!> function N(w), which is the count of eigenvalues of a symmetric
!> tridiagonal matrix T less than or equal to its argument  w.  It
!> performs a choice of two types of loops:
!>
!> IJOB=1, followed by
!> IJOB=2: It takes as input a list of intervals and returns a list of
!>         sufficiently small intervals whose union contains the same
!>         eigenvalues as the union of the original intervals.
!>         The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
!>         The output interval (AB(j,1),AB(j,2)] will contain
!>         eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
!>
!> IJOB=3: It performs a binary search in each input interval
!>         (AB(j,1),AB(j,2)] for a point  w(j)  such that
!>         N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
!>         the search.  If such a w(j) is found, then on output
!>         AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
!>         (AB(j,1),AB(j,2)] will be a small interval containing the
!>         point where N(w) jumps through NVAL(j), unless that point
!>         lies outside the initial interval.
!>
!> Note that the intervals are in all cases half-open intervals,
!> i.e., of the form  (a,b] , which includes  b  but not  a .
!>
!> To avoid underflow, the matrix should be scaled so that its largest
!> element is no greater than  overflow**(1/2) * underflow**(1/4)
!> in absolute value.  To assure the most accurate computation
!> of small eigenvalues, the matrix should be scaled to be
!> not much smaller than that, either.
!>
!> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
!> Matrix", Report CS41, Computer Science Dept., Stanford
!> University, July 21, 1966
!>
!> Note: the arguments are, in general, *not* checked for unreasonable
!> values.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] IJOB
!> \verbatim
!>          IJOB is INTEGER
!>          Specifies what is to be done:
!>          = 1:  Compute NAB for the initial intervals.
!>          = 2:  Perform bisection iteration to find eigenvalues of T.
!>          = 3:  Perform bisection iteration to invert N(w), i.e.,
!>                to find a point which has a specified number of
!>                eigenvalues of T to its left.
!>          Other values will cause SLAEBZ to return with INFO=-1.
!> \endverbatim
!>
!> \param[in] NITMAX
!> \verbatim
!>          NITMAX is INTEGER
!>          The maximum number of "levels" of bisection to be
!>          performed, i.e., an interval of width W will not be made
!>          smaller than 2^(-NITMAX) * W.  If not all intervals
!>          have converged after NITMAX iterations, then INFO is set
!>          to the number of non-converged intervals.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The dimension n of the tridiagonal matrix T.  It must be at
!>          least 1.
!> \endverbatim
!>
!> \param[in] MMAX
!> \verbatim
!>          MMAX is INTEGER
!>          The maximum number of intervals.  If more than MMAX intervals
!>          are generated, then SLAEBZ will quit with INFO=MMAX+1.
!> \endverbatim
!>
!> \param[in] MINP
!> \verbatim
!>          MINP is INTEGER
!>          The initial number of intervals.  It may not be greater than
!>          MMAX.
!> \endverbatim
!>
!> \param[in] NBMIN
!> \verbatim
!>          NBMIN is INTEGER
!>          The smallest number of intervals that should be processed
!>          using a vector loop.  If zero, then only the scalar loop
!>          will be used.
!> \endverbatim
!>
!> \param[in] ABSTOL
!> \verbatim
!>          ABSTOL is REAL
!>          The minimum (absolute) width of an interval.  When an
!>          interval is narrower than ABSTOL, or than RELTOL times the
!>          larger (in magnitude) endpoint, then it is considered to be
!>          sufficiently small, i.e., converged.  This must be at least
!>          zero.
!> \endverbatim
!>
!> \param[in] RELTOL
!> \verbatim
!>          RELTOL is REAL
!>          The minimum relative width of an interval.  When an interval
!>          is narrower than ABSTOL, or than RELTOL times the larger (in
!>          magnitude) endpoint, then it is considered to be
!>          sufficiently small, i.e., converged.  Note: this should
!>          always be at least radix*machine epsilon.
!> \endverbatim
!>
!> \param[in] PIVMIN
!> \verbatim
!>          PIVMIN is REAL
!>          The minimum absolute value of a "pivot" in the Sturm
!>          sequence loop.
!>          This must be at least  max |e(j)**2|*safe_min  and at
!>          least safe_min, where safe_min is at least
!>          the smallest number that can divide one without overflow.
!> \endverbatim
!>
!> \param[in] D
!> \verbatim
!>          D is REAL array, dimension (N)
!>          The diagonal elements of the tridiagonal matrix T.
!> \endverbatim
!>
!> \param[in] E
!> \verbatim
!>          E is REAL array, dimension (N)
!>          The offdiagonal elements of the tridiagonal matrix T in
!>          positions 1 through N-1.  E(N) is arbitrary.
!> \endverbatim
!>
!> \param[in] E2
!> \verbatim
!>          E2 is REAL array, dimension (N)
!>          The squares of the offdiagonal elements of the tridiagonal
!>          matrix T.  E2(N) is ignored.
!> \endverbatim
!>
!> \param[in,out] NVAL
!> \verbatim
!>          NVAL is INTEGER array, dimension (MINP)
!>          If IJOB=1 or 2, not referenced.
!>          If IJOB=3, the desired values of N(w).  The elements of NVAL
!>          will be reordered to correspond with the intervals in AB.
!>          Thus, NVAL(j) on output will not, in general be the same as
!>          NVAL(j) on input, but it will correspond with the interval
!>          (AB(j,1),AB(j,2)] on output.
!> \endverbatim
!>
!> \param[in,out] AB
!> \verbatim
!>          AB is REAL array, dimension (MMAX,2)
!>          The endpoints of the intervals.  AB(j,1) is  a(j), the left
!>          endpoint of the j-th interval, and AB(j,2) is b(j), the
!>          right endpoint of the j-th interval.  The input intervals
!>          will, in general, be modified, split, and reordered by the
!>          calculation.
!> \endverbatim
!>
!> \param[in,out] C
!> \verbatim
!>          C is REAL array, dimension (MMAX)
!>          If IJOB=1, ignored.
!>          If IJOB=2, workspace.
!>          If IJOB=3, then on input C(j) should be initialized to the
!>          first search point in the binary search.
!> \endverbatim
!>
!> \param[out] MOUT
!> \verbatim
!>          MOUT is INTEGER
!>          If IJOB=1, the number of eigenvalues in the intervals.
!>          If IJOB=2 or 3, the number of intervals output.
!>          If IJOB=3, MOUT will equal MINP.
!> \endverbatim
!>
!> \param[in,out] NAB
!> \verbatim
!>          NAB is INTEGER array, dimension (MMAX,2)
!>          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
!>          If IJOB=2, then on input, NAB(i,j) should be set.  It must
!>             satisfy the condition:
!>             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
!>             which means that in interval i only eigenvalues
!>             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
!>             NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with
!>             IJOB=1.
!>             On output, NAB(i,j) will contain
!>             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
!>             the input interval that the output interval
!>             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
!>             the input values of NAB(k,1) and NAB(k,2).
!>          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
!>             unless N(w) > NVAL(i) for all search points  w , in which
!>             case NAB(i,1) will not be modified, i.e., the output
!>             value will be the same as the input value (modulo
!>             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
!>             for all search points  w , in which case NAB(i,2) will
!>             not be modified.  Normally, NAB should be set to some
!>             distinctive value(s) before SLAEBZ is called.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is REAL array, dimension (MMAX)
!>          Workspace.
!> \endverbatim
!>
!> \param[out] IWORK
!> \verbatim
!>          IWORK is INTEGER array, dimension (MMAX)
!>          Workspace.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:       All intervals converged.
!>          = 1--MMAX: The last INFO intervals did not converge.
!>          = MMAX+1:  More than MMAX intervals were generated.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup OTHERauxiliary
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>      This routine is intended to be called only by other LAPACK
!>  routines, thus the interface is less user-friendly.  It is intended
!>  for two purposes:
!>
!>  (a) finding eigenvalues.  In this case, SLAEBZ should have one or
!>      more initial intervals set up in AB, and SLAEBZ should be called
!>      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
!>      Intervals with no eigenvalues would usually be thrown out at
!>      this point.  Also, if not all the eigenvalues in an interval i
!>      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
!>      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
!>      eigenvalue.  SLAEBZ is then called with IJOB=2 and MMAX
!>      no smaller than the value of MOUT returned by the call with
!>      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
!>      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
!>      tolerance specified by ABSTOL and RELTOL.
!>
!>  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
!>      In this case, start with a Gershgorin interval  (a,b).  Set up
!>      AB to contain 2 search intervals, both initially (a,b).  One
!>      NVAL element should contain  f-1  and the other should contain  l
!>      , while C should contain a and b, resp.  NAB(i,1) should be -1
!>      and NAB(i,2) should be N+1, to flag an error if the desired
!>      interval does not lie in (a,b).  SLAEBZ is then called with
!>      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
!>      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
!>      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
!>      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
!>      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
!>      w(l-r)=...=w(l+k) are handled similarly.
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE SLAEBZ(Ijob,Nitmax,N,Mmax,Minp,Nbmin,Abstol,Reltol,    &
     &                  Pivmin,D,E,E2,Nval,Ab,C,Mout,Nab,Work,Iwork,    &
     &                  Info)
      IMPLICIT NONE
!*--SLAEBZ323
!
!  -- LAPACK auxiliary routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      INTEGER Ijob , Info , Minp , Mmax , Mout , N , Nbmin , Nitmax
      REAL Abstol , Pivmin , Reltol
!     ..
!     .. Array Arguments ..
      INTEGER Iwork(*) , Nab(Mmax,*) , Nval(*)
      REAL Ab(Mmax,*) , C(*) , D(*) , E(*) , E2(*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ZERO , TWO , HALF
      PARAMETER (ZERO=0.0E0,TWO=2.0E0,HALF=1.0E0/TWO)
!     ..
!     .. Local Scalars ..
      INTEGER itmp1 , itmp2 , j , ji , jit , jp , kf , kfnew , kl ,     &
     &        klnew
      REAL tmp1 , tmp2
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , MAX , MIN
!     ..
!     .. Executable Statements ..
!
!     Check for Errors
!
      Info = 0
      IF ( Ijob<1 .OR. Ijob>3 ) THEN
         Info = -1
         RETURN
      ENDIF
!
!     Initialize NAB
!
      IF ( Ijob==1 ) THEN
!
!        Compute the number of eigenvalues in the initial intervals.
!
         Mout = 0
         DO ji = 1 , Minp
            DO jp = 1 , 2
               tmp1 = D(1) - Ab(ji,jp)
               IF ( ABS(tmp1)<Pivmin ) tmp1 = -Pivmin
               Nab(ji,jp) = 0
               IF ( tmp1<=ZERO ) Nab(ji,jp) = 1
!
               DO j = 2 , N
                  tmp1 = D(j) - E2(j-1)/tmp1 - Ab(ji,jp)
                  IF ( ABS(tmp1)<Pivmin ) tmp1 = -Pivmin
                  IF ( tmp1<=ZERO ) Nab(ji,jp) = Nab(ji,jp) + 1
               ENDDO
            ENDDO
            Mout = Mout + Nab(ji,2) - Nab(ji,1)
         ENDDO
         RETURN
      ENDIF
!
!     Initialize for loop
!
!     KF and KL have the following meaning:
!        Intervals 1,...,KF-1 have converged.
!        Intervals KF,...,KL  still need to be refined.
!
      kf = 1
      kl = Minp
!
!     If IJOB=2, initialize C.
!     If IJOB=3, use the user-supplied starting point.
!
      IF ( Ijob==2 ) THEN
         DO ji = 1 , Minp
            C(ji) = HALF*(Ab(ji,1)+Ab(ji,2))
         ENDDO
      ENDIF
!
!     Iteration loop
!
      DO jit = 1 , Nitmax
!
!        Loop over intervals
!
         IF ( kl-kf+1>=Nbmin .AND. Nbmin>0 ) THEN
!
!           Begin of Parallel Version of the loop
!
            DO ji = kf , kl
!
!              Compute N(c), the number of eigenvalues less than c
!
               Work(ji) = D(1) - C(ji)
               Iwork(ji) = 0
               IF ( Work(ji)<=Pivmin ) THEN
                  Iwork(ji) = 1
                  Work(ji) = MIN(Work(ji),-Pivmin)
               ENDIF
!
               DO j = 2 , N
                  Work(ji) = D(j) - E2(j-1)/Work(ji) - C(ji)
                  IF ( Work(ji)<=Pivmin ) THEN
                     Iwork(ji) = Iwork(ji) + 1
                     Work(ji) = MIN(Work(ji),-Pivmin)
                  ENDIF
               ENDDO
            ENDDO
!
            IF ( Ijob<=2 ) THEN
!
!              IJOB=2: Choose all intervals containing eigenvalues.
!
               klnew = kl
               DO ji = kf , kl
!
!                 Insure that N(w) is monotone
!
                  Iwork(ji) = MIN(Nab(ji,2),MAX(Nab(ji,1),Iwork(ji)))
!
!                 Update the Queue -- add intervals if both halves
!                 contain eigenvalues.
!
                  IF ( Iwork(ji)==Nab(ji,2) ) THEN
!
!                    No eigenvalue in the upper interval:
!                    just use the lower interval.
!
                     Ab(ji,2) = C(ji)
!
                  ELSEIF ( Iwork(ji)==Nab(ji,1) ) THEN
!
!                    No eigenvalue in the lower interval:
!                    just use the upper interval.
!
                     Ab(ji,1) = C(ji)
                  ELSE
                     klnew = klnew + 1
                     IF ( klnew<=Mmax ) THEN
!
!                       Eigenvalue in both intervals -- add upper to
!                       queue.
!
                        Ab(klnew,2) = Ab(ji,2)
                        Nab(klnew,2) = Nab(ji,2)
                        Ab(klnew,1) = C(ji)
                        Nab(klnew,1) = Iwork(ji)
                        Ab(ji,2) = C(ji)
                        Nab(ji,2) = Iwork(ji)
                     ELSE
                        Info = Mmax + 1
                     ENDIF
                  ENDIF
               ENDDO
               IF ( Info/=0 ) RETURN
               kl = klnew
            ELSE
!
!              IJOB=3: Binary search.  Keep only the interval containing
!                      w   s.t. N(w) = NVAL
!
               DO ji = kf , kl
                  IF ( Iwork(ji)<=Nval(ji) ) THEN
                     Ab(ji,1) = C(ji)
                     Nab(ji,1) = Iwork(ji)
                  ENDIF
                  IF ( Iwork(ji)>=Nval(ji) ) THEN
                     Ab(ji,2) = C(ji)
                     Nab(ji,2) = Iwork(ji)
                  ENDIF
               ENDDO
            ENDIF
!
         ELSE
!
!           End of Parallel Version of the loop
!
!           Begin of Serial Version of the loop
!
            klnew = kl
            DO ji = kf , kl
!
!              Compute N(w), the number of eigenvalues less than w
!
               tmp1 = C(ji)
               tmp2 = D(1) - tmp1
               itmp1 = 0
               IF ( tmp2<=Pivmin ) THEN
                  itmp1 = 1
                  tmp2 = MIN(tmp2,-Pivmin)
               ENDIF
!
               DO j = 2 , N
                  tmp2 = D(j) - E2(j-1)/tmp2 - tmp1
                  IF ( tmp2<=Pivmin ) THEN
                     itmp1 = itmp1 + 1
                     tmp2 = MIN(tmp2,-Pivmin)
                  ENDIF
               ENDDO
!
               IF ( Ijob<=2 ) THEN
!
!                 IJOB=2: Choose all intervals containing eigenvalues.
!
!                 Insure that N(w) is monotone
!
                  itmp1 = MIN(Nab(ji,2),MAX(Nab(ji,1),itmp1))
!
!                 Update the Queue -- add intervals if both halves
!                 contain eigenvalues.
!
                  IF ( itmp1==Nab(ji,2) ) THEN
!
!                    No eigenvalue in the upper interval:
!                    just use the lower interval.
!
                     Ab(ji,2) = tmp1
!
                  ELSEIF ( itmp1==Nab(ji,1) ) THEN
!
!                    No eigenvalue in the lower interval:
!                    just use the upper interval.
!
                     Ab(ji,1) = tmp1
                  ELSEIF ( klnew<Mmax ) THEN
!
!                    Eigenvalue in both intervals -- add upper to queue.
!
                     klnew = klnew + 1
                     Ab(klnew,2) = Ab(ji,2)
                     Nab(klnew,2) = Nab(ji,2)
                     Ab(klnew,1) = tmp1
                     Nab(klnew,1) = itmp1
                     Ab(ji,2) = tmp1
                     Nab(ji,2) = itmp1
                  ELSE
                     Info = Mmax + 1
                     RETURN
                  ENDIF
               ELSE
!
!                 IJOB=3: Binary search.  Keep only the interval
!                         containing  w  s.t. N(w) = NVAL
!
                  IF ( itmp1<=Nval(ji) ) THEN
                     Ab(ji,1) = tmp1
                     Nab(ji,1) = itmp1
                  ENDIF
                  IF ( itmp1>=Nval(ji) ) THEN
                     Ab(ji,2) = tmp1
                     Nab(ji,2) = itmp1
                  ENDIF
               ENDIF
            ENDDO
            kl = klnew
!
         ENDIF
!
!        Check for convergence
!
         kfnew = kf
         DO ji = kf , kl
            tmp1 = ABS(Ab(ji,2)-Ab(ji,1))
            tmp2 = MAX(ABS(Ab(ji,2)),ABS(Ab(ji,1)))
            IF ( tmp1<MAX(Abstol,Pivmin,Reltol*tmp2) .OR. Nab(ji,1)     &
     &           >=Nab(ji,2) ) THEN
!
!              Converged -- Swap with position KFNEW,
!                           then increment KFNEW
!
               IF ( ji>kfnew ) THEN
                  tmp1 = Ab(ji,1)
                  tmp2 = Ab(ji,2)
                  itmp1 = Nab(ji,1)
                  itmp2 = Nab(ji,2)
                  Ab(ji,1) = Ab(kfnew,1)
                  Ab(ji,2) = Ab(kfnew,2)
                  Nab(ji,1) = Nab(kfnew,1)
                  Nab(ji,2) = Nab(kfnew,2)
                  Ab(kfnew,1) = tmp1
                  Ab(kfnew,2) = tmp2
                  Nab(kfnew,1) = itmp1
                  Nab(kfnew,2) = itmp2
                  IF ( Ijob==3 ) THEN
                     itmp1 = Nval(ji)
                     Nval(ji) = Nval(kfnew)
                     Nval(kfnew) = itmp1
                  ENDIF
               ENDIF
               kfnew = kfnew + 1
            ENDIF
         ENDDO
         kf = kfnew
!
!        Choose Midpoints
!
         DO ji = kf , kl
            C(ji) = HALF*(Ab(ji,1)+Ab(ji,2))
         ENDDO
!
!        If no more intervals to refine, quit.
!
         IF ( kf>kl ) EXIT
      ENDDO
!
!     Converged
!
      Info = MAX(kl+1-kf,0)
      Mout = kl
!
!
!     End of SLAEBZ
!
      END SUBROUTINE SLAEBZ


text:
!----------------------------------------------------------------------
! Module: module_apply_convol.f90
!> Purpose: convolution routines
!>

!> Author: Benjamin Menetrier
!>

!> Licensing: this code is distributed under the CeCILL-C license
!>

!> Copyright © 2017 METEO-FRANCE
!----------------------------------------------------------------------
module module_apply_convol

use module_apply_com, only: alpha_com_AB,alpha_com_BA,alpha_com_CA,alpha_copy_AC,alpha_copy_BC
use module_namelist, only: nproc
use omp_lib
use tools_missing, only: msr
use type_fields, only: alphatype
use type_linop, only: apply_linop_sym
use type_mpl, only: mpl
use type_sdata, only: sdatatype,sdatampitype
implicit none

private
public :: convol

contains

!----------------------------------------------------------------------
! Subroutine: convol
!> Purpose: convolution
!----------------------------------------------------------------------
subroutine convol(sdata,alpha)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                 !< Sampling data
type(alphatype),intent(inout) :: alpha(sdata%nproc) !< Subgrid variable

! Local variables
integer :: iproc

if (nproc==0) then
   ! Apply global convolution
   call convol_global(sdata,alpha(1))
elseif (nproc>0) then
   ! Communication
   if (sdata%mpicom==1) then
      ! Copy zone B into zone C
      do iproc=1,sdata%nproc
         call alpha_copy_BC(sdata%mpi(iproc),alpha(iproc))
      end do
   elseif (sdata%mpicom==2) then
      ! Halo reduction from zone B to zone A
      call alpha_com_BA(sdata,alpha)

      ! Copy zone A into zone C
      do iproc=1,sdata%nproc
         call alpha_copy_AC(sdata%mpi(iproc),alpha(iproc))
      end do
   end if

   ! Apply local convolution
   do iproc=1,sdata%nproc
      call convol_local(sdata%mpi(iproc),alpha(iproc))
   end do
end if

! Communication
if (nproc>0) then
   ! Halo reduction from zone C to zone A
   call alpha_com_CA(sdata,alpha)

   ! Halo extension from zone A to zone B
   call alpha_com_AB(sdata,alpha)
end if

end subroutine convol

!----------------------------------------------------------------------
! Subroutine: convol_global
!> Purpose: convolution, global
!----------------------------------------------------------------------
subroutine convol_global(sdata,alpha)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata    !< Sampling data
type(alphatype),intent(inout) :: alpha !< Subgrid variable

! Apply linear operator, symmetric
call apply_linop_sym(sdata%c,alpha%val)

end subroutine convol_global

!----------------------------------------------------------------------
! Subroutine: convol_local
!> Purpose: convolution, local
!----------------------------------------------------------------------
subroutine convol_local(sdatampi,alpha)

implicit none

! Passed variables
type(sdatampitype),intent(in) :: sdatampi !< Sampling data
type(alphatype),intent(inout) :: alpha    !< Subgrid variable

! Apply linear operator, symmetric
call apply_linop_sym(sdatampi%c,alpha%valc)

end subroutine convol_local

end module module_apply_convol

max_stars_repo_name: thomasauligne/NICAS

max_issues_repo_name: thomasauligne/NICAS

max_forks_repo_name: thomasauligne/NICAS