text:
subroutine cont_weath(t)
!************************
use co2_interpolation, only: climate_interpolation
use dynsoil, only: regolith_time_integration
use dynsoil_steady_state, only: steady_state_surface_regolith
use constante, only: PI_n_CO2_atm, Rgas, TSS_density, rsw_litho, CaMg_rock, P_rock, P2C_ker, P2C_carb, OC_rock, cp_cont, &
                     Sulf_rock, BASALT_LITHO_NUM
!use dynsoil_lithium, only: lithium
implicit none
double precision :: loc_granwth
include 'combine_foam.inc'
! WARNING: carbonate must be the last lithological class, whatever the chosen
! number of classes

! current atmospheric CO2 level (in PAL) and climate interpolation (from look-up table)
! -------------------------------------------------------------------------------------

p=var(12,nbasin)/PI_n_CO2_atm
call climate_interpolation(p)

! Total DISCHARGE of water to the ocean
!--------------------------------------

discharge=0.

do j0=1,ncontpxl
    j = list_cont_pixel(j0)
    discharge = discharge + 1e10*runclim(j)*areaclimber(j)
end do !                    ^^^^^
! convert area [1e6km2 -> m2] and runoff [cm/y -> m/y]

!========================================================!
!========================================================!
!==               WEATHERING COMPUTATION               ==!
!========================================================!
!========================================================!

call FACCO2(t,pco2,fco2,fe)

!=========================================================================!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! VARIABLES DEPENDENT OF DYNSOIL USAGE !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if (.not. coupling_dynsoil) then
! ######################################
! # if not coupled with dynsoil module #
! ######################################

    ! initialising:
    fsilw=0.
    fbasw=0.
    fkerw=0.
    total_cont_POC_export=0.
    weighted_rsw=0.

    ! Riverine sediment discharge => parameteric power law of water discharge
    tss = 2.9d6 * discharge**0.5

    do j0=1,ncontpxl
        j = list_cont_pixel(j0)

        ! =====================
        ! Basaltic weathering
        ! =====================
        ! "basalt" lithology => number BASALT_LITHO_NUM
        facbas=(-42300./Rgas)*(1./(Tclim(j)+273.15)-1./288.15)
        wth_litho(BASALT_LITHO_NUM,j) = (0.0483d+10)*dexp(facbas)*runclim(j) !0.044086200331d+10  !483685972.093689  !172421904
        wth_litho_wgh(BASALT_LITHO_NUM,j) = wth_litho(BASALT_LITHO_NUM,j) * litho_frac(BASALT_LITHO_NUM,j)
        fbas(j) = wth_litho_wgh(BASALT_LITHO_NUM,j) * areaclimber(j)
        fbasw=fbasw+fbas(j)*clo

        ! =====================
        ! Granitic weathering
        ! =====================
        ! skip basalts (BASALT_LITHO_NUM) and carbonates (last litho)
        facgra=(-48200./Rgas)*(1./(Tclim(j)+273.15)-1./288.15)
        loc_granwth = (0.11379d+10)*dexp(facgra)*runclim(j) * fco2*fe !0.101357275d+10
        fsil(j) = 0
        do k = 1,BASALT_LITHO_NUM-1 ! skip basalts
            wth_litho(k,j) = loc_granwth
            wth_litho_wgh(k,j) = wth_litho(k,j) * litho_frac(k,j)
            fsil(j) = fsil(j) + wth_litho_wgh(k,j)
            weighted_rsw = weighted_rsw + rsw_litho(k)*wth_litho_wgh(k,j)*1e12*areaclimber(j)
        end do
        do k = BASALT_LITHO_NUM+1,nlitho-1 ! skip carbonates
            wth_litho(k,j) = loc_granwth
            wth_litho_wgh(k,j) = wth_litho(k,j) * litho_frac(k,j)
            fsil(j) = fsil(j) + wth_litho_wgh(k,j)
            weighted_rsw = weighted_rsw + rsw_litho(k)*wth_litho_wgh(k,j)*1e12*areaclimber(j)
        end do
        fsil(j) = fsil(j) * areaclimber(j)
        fsilw=fsilw+fsil(j)*clo

        ! conversion /1e6km2 => /m2 for ALL silicate weathering rates
        wth_litho(1:nlitho-1,j)     = 1e-12*wth_litho(1:nlitho-1,j)
        wth_litho_wgh(1:nlitho-1,j) = 1e-12*wth_litho_wgh(1:nlitho-1,j)

        wth_allsil(j) = sum(wth_litho_wgh(1:nlitho-1,j))

        ! =====================
        ! Kerogen weathering
        ! =====================
        ! Old formulation:
        !! facker=(-43000./Rgas)*(1./(Tclim(j)+273.15)-1./288.15)
        !fker(j)=4.4184d+8*runclim(j)*areaclimber(j)   !*var(11,nbasin)/38.d+18 !& !2.4 !1.432
        !                                              ! *dexp(facker)  !0.672d+10
        !fkerw=fkerw+fker(j)*clo
        fker(j) = 0.5 * sum(litho_frac(:,j)*OC_rock)*(tss/TSS_density/areatot)
        !         ^^^ : oxidation efficiency
        if (.not. lock_oxygen_cycle)  fkerw = fkerw + fker(j)*areaclimber(j)*clo
        fker(j) = 1e-12*fker(j) ! conversion /1e6km2 => /m2

    end do

    ! ======================
    ! Continental POC export
    ! ======================

    ! erosion export (t/km2/an)
    reg_eros_galy_unit=(tss/areatot)*1e-6*1e-3  !moving from kg/1e6km2/yr to t/km2/yr
    do j0=1,ncontpxl
        j = list_cont_pixel(j0)
        !scaling:
        POC_export_rate(j)=0.081*(reg_eros_galy_unit)**0.56  !in t/km2/yr
        POC_export_rate(j)=POC_export_rate(j)/12.         !mol/m2/yr
        total_cont_POC_export=total_cont_POC_export+POC_export_rate(j)*areaclimber(j)*1.e+12  !mol/yr
    end do

    ! Lithium riverine flux^M
    FrivLi=10.d+9

    ! lithology-averaged Strontium isotopic ratio
    weighted_rsw = weighted_rsw/fsilw

    ! ===============================================================================
    ! Pyrite weathering and distribution in carb dissol, sil dissol and H2SO4 release
    ! ===============================================================================

    ! molar flux of SO4^2- from pyrite oxidative weathering:
    fpyrw = 0.025*tss
    ! molar flux of Ca released to the ocean from carbonate dissolution by sulfuric acid
    fcarbsulfw  =  0.724*fpyrw
    ! molar flux of Ca released to the ocean from silicate dissolution by sulfuric acid
    fsilsulfw   =  0.276*fpyrw
    ! molar flux of H2SO4 from pyrite oxidation directly reaching the ocean
    fH2SO4sulfw =  0*fpyrw

else
! ######################################
! #   if coupled with dynsoil module   #
! ######################################

    ! increase time counter of dynsoil (=> asynchronous coupling) !
    icount_DS_int = icount_DS_int + 1
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

    if ( icount_DS_int >= ijump_DS_integration ) then ! if integration time has been reached

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !! integrate dynsoil module !!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ! counter reinitialization:
        icount_DS_int = 0

        ! =============================
        ! DYNSOIL (silicate weathering)
        ! =============================
        !
        if (use_dynsoil_steady_state) then
            call steady_state_surface_regolith(reg_thick, reg_x_surf, reg_tau_surf, reg_prod, reg_eros, reg_P_diss, &
                                               Tclim, runclim, slope, DS_timestep, veget_factor, veget_eros_factor, &
                                               list_cont_pixel, ncontpxl                                            )
        else
            call regolith_time_integration(xlevs, reg_thick, reg_x_surf, reg_tau_surf, reg_P_vol, reg_ktop,                 &
                                           reg_z_prof, reg_tau_prof,                                                        &
                                           reg_prod,reg_eros,reg_P_diss,reg_P_eros, reg_x_surf_eros, Tclim, runclim, slope, &
                                           DS_timestep, DS_varout_missval(9), DS_varout_missval(10),                        &
                                           veget_factor, veget_eros_factor, list_cont_pixel, ncontpxl                       )
        end if
        !call lithium(reg_P_diss, Tclim, runclim, reg_thick, reg_Li_Friv, reg_Li_Fsp, reg_Li_driv, list_cont_pixel, ncontpxl)

        ! Continental sum:
        ! ----------------

        ! initialising:
        fsilw=0.
        fbasw=0.
        FrivLi=0.
        dLiriv=0.
        tss=0.
        fkerw=0.
        total_cont_POC_export=0.
        weighted_rsw=0.
        fcarbsulfw=0
        fsilsulfw=0
        fH2SO4sulfw=0
        faddsulfw=0

        do j0=1,ncontpxl
            j = list_cont_pixel(j0)

            wth_allsil(j) = 0
            fsil(j) = 0
            do k = 1,nlitho-1 ! as k=nlitho is carbonate
                !
                wth_litho(k,j) = CaMg_rock(k) * reg_P_diss(k,j)
                wth_litho_wgh(k,j) = wth_litho(k,j) * litho_frac(k,j)
                wth_allsil(j) = wth_allsil(j) + wth_litho_wgh(k,j)
                fsil(j) = fsil(j)  +  wth_litho_wgh(k,j) * 1e12*areaclimber(j)
                !       WARNING: conversion: 1e6km2 -> m2: ^^^^
                weighted_rsw = weighted_rsw + rsw_litho(k)*wth_litho_wgh(k,j)*1e12*areaclimber(j)
                FrivLi = FrivLi + reg_Li_Friv(k,j)*litho_frac(k,j)*1e12*areaclimber(j)   !mol/yr
                dLiriv = dLiriv + reg_Li_driv(k,j)*(reg_Li_Friv(k,j)*litho_frac(k,j)*1e12*areaclimber(j))
            end do
            !
            ! Impose erosion of carbonates = erosion of sediments
            reg_eros(nlitho,j) = reg_eros(nlitho-1,j)
            !
            ! riverine sediment discharge
            tss = tss + TSS_density*sum(reg_eros(:,j)*litho_frac(:,j))*1e12*areaclimber(j)
            !
            ! Split silicates weathering in 2 fluxes: "basalt" weathering (mafic rocks) and "granite" weathering (rest of silicates)
            fbas(j) =  wth_litho_wgh(BASALT_LITHO_NUM,j) * 1e12*areaclimber(j)
            fsil(j) = fsil(j) - fbas(j)
            fsilw = fsilw + fsil(j)*clo
            fbasw = fbasw + fbas(j)*clo
            !
            ! same for strontium isotopic ratio
            weighted_rsw = weighted_rsw - rsw_litho(BASALT_LITHO_NUM)*wth_litho_wgh(BASALT_LITHO_NUM,j)*1e12*areaclimber(j)

            ! =====================
            ! Kerogen weathering
            ! =====================
            ! Old formulation:
            ! facker=(-43000./Rgas)*(1./(Tclim(j)+273.15)-1./288.15)
            ! fker(j)=4.4184d+8*runclim(j)*areaclimber(j)   !*var(11,nbasin)/38.d+18 !& !2.4 !1.432
                                                            ! *dexp(facker)  !0.672d+10
            fker(j) = 0.5 * sum(litho_frac(:,j)*OC_rock*reg_eros(:,j))
            !         ^^^ : oxidation efficiency
            ! *************************************
            ! Kerogen weathering perturbation: +50%
            !fker(j) = 1.5 * fker(j)
            ! +10.32% (eq. in C source as + 50% of sulfuric weathering)
            !fker(j) = 1.1032 * fker(j)
            ! *************************************
            if (.not. lock_oxygen_cycle)  fkerw = fkerw + fker(j)*1e12*areaclimber(j)*clo

            ! ======================
            ! Continental POC export
            ! ======================

            ! erosion export (t/km2/an)
            reg_eros_lithmean(j) = sum(litho_frac(:,j)*reg_eros(:,j)) / sum(litho_frac(:,j))
            !scaling:
            reg_eros_galy_unit=reg_eros_lithmean(j)*TSS_density*1.e6*1.e-3  !moving from m/yr to t/km2/yr
            POC_export_rate(j)=0.081*(reg_eros_galy_unit)**0.56  !in t/km2/yr
            POC_export_rate(j)=POC_export_rate(j)/12.         !mol/m2/yr
            !! add O2 feedback
            !POC_export_rate(j) = POC_export_rate(j) / sqrt(var(11,nbasin)/38.002d18)
            !! add stronger O2 feedback
            !POC_export_rate(j) = POC_export_rate(j) / (var(11,nbasin)/38.002d18)
            total_cont_POC_export=total_cont_POC_export+POC_export_rate(j)*areaclimber(j)*1.e+12  !mol/yr

            ! ===============================================================================
            ! Pyrite weathering and distribution in carb dissol, sil dissol and H2SO4 release
            ! ===============================================================================

            ! molar flux of SO4^2- from pyrite oxidative weathering:
            fpyrw = sum(Sulf_rock(:) * reg_eros(:,j) * litho_frac(:,j))
            ! molar flux of Ca released to the ocean from carbonate dissolution by sulfuric acid
            fcarbsulfw  = fcarbsulfw    +  0.724*fpyrw * 1e12*areaclimber(j)
            ! molar flux of Ca released to the ocean from silicate dissolution by sulfuric acid
            fsilsulfw   = fsilsulfw     +  0.276*fpyrw * 1e12*areaclimber(j)
            ! molar flux of H2SO4 from pyrite oxidation directly reaching the ocean
            fH2SO4sulfw = fH2SO4sulfw   +  0*fpyrw  *1e12*areaclimber(j)

            !*********************************************************************!
            ! SULFURIC WEATHERING PERTURBATION - PART 1/3
            ! -------------------------------------------
            ! PERTURBATION: +50% of pyrite weathering flux (with kerw => +10.32%):
            ! abrupt perturbation
            !faddsulfw = faddsulfw + 0.5*fpyrw*1e12*areaclimber(j)
            !faddsulfw = faddsulfw + 0.1032*fpyrw*1e12*areaclimber(j)
            ! progressive pertubation over 40 Ma
            !faddsulfw = faddsulfw + 0.5*min(1d0, t/40d6)*fpyrw*1e12*areaclimber(j)
            !faddsulfw = faddsulfw + 0.1032*min(1d0, t/40d6)*fpyrw*1e12*areaclimber(j)
            !*********************************************************************!

        end do

        !*********************************************************************!
        ! SULFURIC WEATHERING PERTURBATION - PART 2/3
        ! -------------------------------------------
        !
        !fcarbsulfw  = fcarbsulfw  + faddsulfw
        !fsilsulfw   = fsilsulfw   + faddsulfw
        !fH2SO4sulfw = fH2SO4sulfw + faddsulfw
        !
        ! Trade-off between sulfuric weathering and carbonic weathering:
        ! => deduce additional sulfuric weathering from carbonic weathering
        ! (Silicate sulfuric weathering only)
        !fsilw  = fsilw  - faddsulfw
        !*********************************************************************!

        ! Lithium riverine isotopic delta
        !dLiriv = dLiriv/FrivLi

        ! lithology-averaged Strontium isotopic ratio
        weighted_rsw = weighted_rsw/fsilw

    end if

end if

!=========================================================================!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! VARIABLES INDEPENDENT OF DYNSOIL USAGE !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! initialising:
fcarbw=0.
fpw=0.
fpyrw=0.
ftrapw=0.

do j0=1,ncontpxl
    j = list_cont_pixel(j0)

    ! =====================
    ! Carbonate weathering
    ! =====================
    !fcarb(j)=0.18577793d+11*dsqrt(runclim(j))*areaclimber(j)*fco2*fe !old
    call carbo(p,Tclim(j),runclim(j),carb_weath_conc,j)
    ! Carbonate is assumed to be the last lithological class (k = nlitho)
    wth_litho(nlitho,j) = 3.589d10*runclim(j)*carb_weath_conc ! to fit 12.3 Tmol/yr of Carb wth
    !wth_litho(nlitho,j) = 16.*0.483483601042728d+10*runclim(j)*carb_weath_conc*0.978903 !full newton, old tuning
    wth_litho_wgh(nlitho,j) = litho_frac(nlitho,j) * wth_litho(nlitho,j)
    fcarb(j) = areaclimber(j) * wth_litho_wgh(nlitho,j)
    ! conversion because areaclimber is in 1e6km2 instead of m2:
    wth_litho(nlitho,j) = 1e-12 * wth_litho(nlitho,j)
    wth_litho_wgh(nlitho,j) = 1e-12 * wth_litho_wgh(nlitho,j)

    fcarbw=fcarbw+fcarb(j)*clo

    ! =====================
    ! Phosphorus weathering
    ! =====================
    fp(j) = sum((P_rock(1:nlitho-1)/CaMg_rock(1:nlitho-1)) * wth_litho_wgh(1:nlitho-1,j)) & ! silicate part
            +  P2C_carb * wth_litho_wgh(nlitho,j) & ! carbonate part
            +  P2C_ker * fker(j) ! kerogen part
            ! *************************************
            ! Do not consider "extra" kerogen weathering (ie, perturbation) for P weathering
            !+  P2C_ker * fker(j) / 1.5    ! kerogen part
            !+  P2C_ker * fker(j) / 1.1032 ! kerogen part
            ! *************************************
    fpw = fpw + fp(j)*1d12*areaclimber(j)*clo*phosss

end do

! Substract Phosphorus that is in exported biospheric organic matter
fpw = fpw - total_cont_POC_export/cp_cont

! + + + + + + + + + !
! Lock sulfur cycle !
! + + + + + + + + + !
if (lock_sulfur_cycle) then
    ! => keep sil. sulf. wth (because it eventually consumes CO2), set H2SO4 release to 0 and
    !    use carb. sulf. wth as adjustment variable to balance sulfate-reduction
    fH2SO4sulfw = 0
    fcarbsulfw = sum( fSulfRed(1:nbasin-1) ) - fsilsulfw
end if

! + + + + + + + + + !
! Lock oxygen cycle !
! + + + + + + + + + !
if (lock_oxygen_cycle)  fkerw = sum( fO2_odc(1:nbasin-1) ) - (15./8.)*(fcarbsulfw + fsilsulfw + fH2SO4sulfw)

!*********************************************************************!
! SULFURIC WEATHERING PERTURBATION - PART 3/3
! -------------------------------------------
! Trade-off between sulfuric weathering and carbonic weathering:
! => deduce additional sulfuric weathering from carbonic weathering
! (Carbonate sulfuric weathering only)
!fcarbw = fcarbw - faddsulfw
!*********************************************************************!

return
end


text:
      SUBROUTINE slF45Z (R1950,D1950,BEPOCH,R2000,D2000)
*+
*     - - - - - -
*      F 4 5 Z
*     - - - - - -
*
*  Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
*  proper motion in the FK5 frame (double precision)
*
*  This routine converts stars from the old, Bessel-Newcomb, FK4
*  system to the new, IAU 1976, FK5, Fricke system, in such a
*  way that the FK5 proper motion is zero.  Because such a star
*  has, in general, a non-zero proper motion in the FK4 system,
*  the routine requires the epoch at which the position in the
*  FK4 system was determined.
*
*  The method is from Appendix 2 of Ref 1, but using the constants
*  of Ref 4.
*
*  Given:
*     R1950,D1950     dp    B1950.0 FK4 RA,Dec at epoch (rad)
*     BEPOCH          dp    Besselian epoch (e.g. 1979.3D0)
*
*  Returned:
*     R2000,D2000     dp    J2000.0 FK5 RA,Dec (rad)
*
*  Notes:
*
*  1)  The epoch BEPOCH is strictly speaking Besselian, but
*      if a Julian epoch is supplied the result will be
*      affected only to a negligible extent.
*
*  2)  Conversion from Besselian epoch 1950.0 to Julian epoch
*      2000.0 only is provided for.  Conversions involving other
*      epochs will require use of the appropriate precession,
*      proper motion, and E-terms routines before and/or
*      after FK45Z is called.
*
*  3)  In the FK4 catalogue the proper motions of stars within
*      10 degrees of the poles do not embody the differential
*      E-term effect and should, strictly speaking, be handled
*      in a different manner from stars outside these regions.
*      However, given the general lack of homogeneity of the star
*      data available for routine astrometry, the difficulties of
*      handling positions that may have been determined from
*      astrometric fields spanning the polar and non-polar regions,
*      the likelihood that the differential E-terms effect was not
*      taken into account when allowing for proper motion in past
*      astrometry, and the undesirability of a discontinuity in
*      the algorithm, the decision has been made in this routine to
*      include the effect of differential E-terms on the proper
*      motions for all stars, whether polar or not.  At epoch 2000,
*      and measuring on the sky rather than in terms of dRA, the
*      errors resulting from this simplification are less than
*      1 milliarcsecond in position and 1 milliarcsecond per
*      century in proper motion.
*
*  References:
*
*     1  Aoki,S., et al, 1983.  Astron.Astrophys., 128, 263.
*
*     2  Smith, C.A. et al, 1989.  "The transformation of astrometric
*        catalog systems to the equinox J2000.0".  Astron.J. 97, 265.
*
*     3  Yallop, B.D. et al, 1989.  "Transformation of mean star places
*        from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
*        Astron.J. 97, 274.
*
*     4  Seidelmann, P.K. (ed), 1992.  "Explanatory Supplement to
*        the Astronomical Almanac", ISBN 0-935702-68-7.
*
*  Called:  slDS2C, slEPJ, slEB2D, slDC2S, slDA2P
*
*  P.T.Wallace   Starlink   21 September 1998
*
*  Copyright (C) 1998 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 R1950,D1950,BEPOCH,R2000,D2000

      DOUBLE PRECISION D2PI
      PARAMETER (D2PI=6.283185307179586476925287D0)

      DOUBLE PRECISION W
      INTEGER I,J

*  Position and position+velocity vectors
      DOUBLE PRECISION R0(3),A1(3),V1(3),V2(6)

*  Radians per year to arcsec per century
      DOUBLE PRECISION PMF
      PARAMETER (PMF=100D0*60D0*60D0*360D0/D2PI)

*  Functions
      DOUBLE PRECISION slEPJ,slEB2D,slDA2P

*
*  CANONICAL CONSTANTS  (see references)
*

*  Vectors A and Adot, and matrix M (only half of which is needed here)
      DOUBLE PRECISION A(3),AD(3),EM(6,3)
      DATA A,AD/ -1.62557D-6,  -0.31919D-6, -0.13843D-6,
     :           +1.245D-3,    -1.580D-3,   -0.659D-3/

      DATA (EM(I,1),I=1,6) / +0.9999256782D0,
     :                       +0.0111820610D0,
     :                       +0.0048579479D0,
     :                       -0.000551D0,
     :                       +0.238514D0,
     :                       -0.435623D0 /

      DATA (EM(I,2),I=1,6) / -0.0111820611D0,
     :                       +0.9999374784D0,
     :                       -0.0000271474D0,
     :                       -0.238565D0,
     :                       -0.002667D0,
     :                       +0.012254D0 /

      DATA (EM(I,3),I=1,6) / -0.0048579477D0,
     :                       -0.0000271765D0,
     :                       +0.9999881997D0,
     :                       +0.435739D0,
     :                       -0.008541D0,
     :                       +0.002117D0 /

*  Spherical to Cartesian
      CALL slDS2C(R1950,D1950,R0)

*  Adjust vector A to give zero proper motion in FK5
      W=(BEPOCH-1950D0)/PMF
      DO I=1,3
         A1(I)=A(I)+W*AD(I)
      END DO

*  Remove e-terms
      W=R0(1)*A1(1)+R0(2)*A1(2)+R0(3)*A1(3)
      DO I=1,3
         V1(I)=R0(I)-A1(I)+W*R0(I)
      END DO

*  Convert position vector to Fricke system
      DO I=1,6
         W=0D0
         DO J=1,3
            W=W+EM(I,J)*V1(J)
         END DO
         V2(I)=W
      END DO

*  Allow for fictitious proper motion in FK4
      W=(slEPJ(slEB2D(BEPOCH))-2000D0)/PMF
      DO I=1,3
         V2(I)=V2(I)+W*V2(I+3)
      END DO

*  Revert to spherical coordinates
      CALL slDC2S(V2,W,D2000)
      R2000=slDA2P(W)

      END

ext: f

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text: Katie Johnson was the 20052006 Academic Affairs Commission chair.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!========================================================================
!
!                          T o m o f a s t - x
!                        -----------------------
!
!              Main authors: Vitaliy Ogarko, Roland Martin,
!                   Jeremie Giraud, Dimitri Komatitsch.
!
! (c) CNRS, France, and University of Western Australia. January 2018
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

!======================================================================================================
! In the Finite Volume formulation one has the nodes of the unknowns located at the center of a cell,
! and the nodes of the fluxes located in the middle of each edge of the cell.
! Permittivities are defined at the same location at a corner of a cell.
! It is necessary to do that because one can not impose two different permittivities
! at a sharp interface of the tube (at the electrodes and at i = i1).
!======================================================================================================
module geometry

  use global_typedefs
  use parameters_ect

  implicit none

  private

  public :: get_refined_ntheta
  public :: set_geometry
  public :: compute_grid_coordinates
  public :: compute_shifted_grid_coordinates
  public :: refine_geometry
  public :: calculate_cell_volumes

  private :: calc_electrode_coordinates

contains

!========================================================
! Returns the number of points per electrode.
!========================================================
function get_npelec(ntheta, nel)
  integer, intent(in) :: ntheta, nel
  integer :: get_npelec

  get_npelec = 2*ntheta/nel/3
  return
end function get_npelec

!========================================================
! Returns the number of points per gaps.
!========================================================
function get_npgaps(ntheta, nel)
  integer, intent(in) :: ntheta, nel
  integer :: get_npgaps

  get_npgaps = ntheta/nel - get_npelec(ntheta, nel)
  return
end function get_npgaps

!========================================================
! Returns the ntheta for the refined mesh.
!========================================================
function get_refined_ntheta(ntheta, nel)
  integer, intent(in) :: ntheta, nel
  integer :: get_refined_ntheta

  get_refined_ntheta = ntheta + get_npelec(ntheta,nel)-1 + 2*(get_npgaps(ntheta,nel)+1)
  return
end function get_refined_ntheta

!=============================================================================================
! Calculate the bottom and top (z-direction) grid coordinates of the electrodes.
!=============================================================================================
subroutine calc_electrode_coordinates(nrings, nel, kguards, elecs)
  integer, intent(in) :: nrings, nel, kguards
  type(t_electrode), intent(inout) :: elecs(:)

  integer :: ielectrode
  integer :: height, iring

  ! TODO: Add (vertical) gaps between the rings.

  if (mod(2 * kguards, nrings) /= 0) then
    print *, "2 * kguards is not a multiple of nrings; exiting..."
    stop
  endif

  if (mod(nel, nrings) /= 0) then
    print *, "nel is not a multiple of nrings; exiting..."
    stop
  endif

  height = 2 * kguards / nrings

  do ielectrode = 1, nel
    ! The ring number to which this electrode belongs, e.g.,
    ! 1-12 electrodes belong to the 1st rings, 13-24 electrodes to the 2nd ring etc.
    iring = (ielectrode - 1) / (nel / nrings) + 1

    elecs(ielectrode)%izmin = kguards + 1 + (iring - 1) * height
    elecs(ielectrode)%izmax = elecs(ielectrode)%izmin + height - 1

    !print *, '---->', ielectrode, iring, elecs(ielectrode)%izmin, elecs(ielectrode)%izmax
  enddo

end subroutine calc_electrode_coordinates

!=============================================================================================
! Sets the system geometry.
!=============================================================================================
subroutine set_geometry(ilevel, par, nr, ntheta, nz, r, dr, theta, dz, z, izspace, kguards, &
                        i1, i2, elecs, guards)
  ! Parameters.
  type(t_parameters_ect), intent(in) :: par
  ! Dimensions (locally).
  integer, intent(in) :: nr, ntheta, nz, kguards
  integer, intent(in) :: ilevel
  ! theta, r and z angle, radius and size of each point, step along r and z
  real(kind=CUSTOM_REAL), intent(out) :: r(0:)
  real(kind=CUSTOM_REAL), intent(out) :: dr(0:)
  real(kind=CUSTOM_REAL), intent(out) :: theta(0:)
  real(kind=CUSTOM_REAL), intent(out) :: dz(0:)
  real(kind=CUSTOM_REAL), intent(out) :: z(0:)
  integer, intent(out) :: izspace
  integer, intent(out) :: i1, i2
  type(t_electrode), intent(out) :: elecs(:)
  type(t_guards), intent(out) :: guards

  ! Local variables --------------------------------------
  real(kind=CUSTOM_REAL) :: dtheta_elec, dtheta_gap, heiguard
  real(kind=CUSTOM_REAL) :: dr1, dr2, dr3, dz1, dz2, dz3, dz4, dz5
  integer :: i, j, k, i1_read, i2_read
  integer :: npelec, npgaps, ielec, ielectrode, jelectrode
  integer :: nr1, nr2
  integer :: igaps_on_theta,igaps_on_z
  integer :: jleft
  integer :: nelr

  !=======================================================
  ! Compute radii.
  !=======================================================

  if (par%ifixed_elecgeo == 1) then
    nr1 = int(par%sens%radiusin * par%dims%nr / par%sens%radiusoutout)
    nr2 = int((par%sens%radiusout - par%sens%radiusin) * par%dims%nr / par%sens%radiusoutout)
    i1_read = nr1-mod(nr1,8)
    i2_read = i1_read+nr2-mod(nr2,8)
    nr1 = i1_read
    nr2 = i2_read

    i1 = i1_read/2**(ilevel-1)
    i2 = i2_read/2**(ilevel-1)
    nr1 = i1
    nr2 = i2
  else
    i2 = 8*nr/9
    i1 = 3*i2/4
    nr1 = i1
    nr2 = i2
  endif

! VO VO RM RM 6 Jan 2015 fixed the mapping according to Mapping #3 in the document Mappings.pdf
! VO VO RM RM cell index i --> radius index 2i-1,
! VO VO RM RM with ghost points i=0 and i=nr+1 separated by 2dr from regular points i=1 and i=nr.
! VO VO RM RM To be consistent with matrix_init and r(2i1-1)=R1, r(2i2-1)=R2 and r(2(nr+1)-1)=R3
  dr1 = (par%sens%radiusin) / dble(2*nr1)
  dr2 = (par%sens%radiusout - par%sens%radiusin) / dble(2*nr2-2*nr1)
  dr3 = (par%sens%radiusoutout - par%sens%radiusout) / dble(2*nr-2*nr2+2)

  do i=0,2*i1-1
    r(i) = dr1+dble(i)*dr1
  enddo

  do i=2*i1,2*i2-1
    r(i) = par%sens%radiusin+dble(i-2*nr1+1)*dr2
  enddo

  do i=2*i2,2*nr+1
    r(i) = par%sens%radiusout+dble(i-2*nr2+1)*dr3
  enddo
! VO VO RM RM 6 Jan 2015 End of the fix.

   !Set equal dr everywhere.
!  do i=0,2*nr+1
!    r(i) = dble(i)*sens%radiusoutout/dble(2*nr+1)
!  enddo

! VO VO 9 Jan fixed dr(0)=r(0)-r(-1) according to radius index Mapping #3.
  dr(0) = r(0)
  do i=1,2*nr+1
    dr(i) = r(i)-r(i-1)
  enddo

  !=======================================================
  ! Compute theta angles.
  !=======================================================

  ! Number of electrodes per ring.
  nelr = par%nel / par%nrings

  igaps_on_theta = nint(par%sens%space_electrodes)
  !if (myrank == 0) print *,'ilevel, igaps_on_theta=', ilevel, igaps_on_theta

! VO VO 8 Jan 2015 Fixed to have the same angle at the tips of electrodes at different level.
! VO VO 8 Jan 2015 Added unit tests for this check marked with (***)
! VO VO 8 Jan 2015 Made the calculation robust, i.e., avoiding iterative summation in the loop.
! VO VO 14 Jan 2015 Made sure that there is at least one point on the gap, i.e., a point without Dirichlet BCs.
! VO VO So for ntheta=36 electrode points will be 1-2, 4-5, 7-8, etc.
  if (igaps_on_theta /= 0) then
    ! The number of points per electrode.
    npelec = get_npelec(ntheta, nelr)
    ! The number of points per gap, i.e., points without Dirichlet BCs.
    npgaps = get_npgaps(ntheta, nelr)
    ! Compute angles.
    dtheta_elec = (2._CUSTOM_REAL*PI/nelr - par%sens%space_electrodes*PI/180._CUSTOM_REAL) / dble(npelec-1)
    dtheta_gap = par%sens%space_electrodes*PI/180._CUSTOM_REAL / dble(npgaps+1)

    theta(0) = dtheta_gap
    do j=1,ntheta+1
      ielec = (j-1) / (npelec+npgaps)
      jleft = (j-1) - ielec*(npelec+npgaps)

      if (jleft <= npelec-1) then
        theta(j) = -2._CUSTOM_REAL*PI/nelr*ielec - dtheta_elec*jleft
      else
        theta(j) = -2._CUSTOM_REAL*PI/nelr*ielec - dtheta_elec*(npelec-1) - dtheta_gap*(jleft-npelec+1)
      endif
    enddo

    ! Calculate the first and the last indices in the theta-array for the electrodes.
    do ielectrode=1,nelr
      elecs(ielectrode)%ithetamin = (ielectrode-1)*(npelec+npgaps)+1
      elecs(ielectrode)%ithetamax = (ielectrode-1)*(npelec+npgaps)+npelec
    enddo
! VO VO 14 Jan 2015 End of the fix.

  else
    ! The number of points per electrode.
    npelec = ntheta/nelr
    ! Compute angle.
    dtheta_elec = (2._CUSTOM_REAL*PI/nelr)/dble(npelec)

    do j=0,ntheta+1
      theta(j) = dtheta_elec*dble(1-j)
    enddo

    ! Calculate the first and the last indices in the theta-array for the electrodes.
    do ielectrode=1,nelr
      elecs(ielectrode)%ithetamin = (ielectrode-1)*(npelec)+1
      elecs(ielectrode)%ithetamax = (ielectrode-1)*(npelec)+npelec
    enddo
  endif
! VO VO 8 Jan 2015 End of the fix.

  if (par%nrings > 1) then
  ! Copy the theta angles from the first to the other electrode rings.
    do ielectrode = nelr + 1, par%nel
      jelectrode = mod(ielectrode, nelr)
      if (jelectrode == 0) jelectrode = nelr

      elecs(ielectrode)%ithetamin = elecs(jelectrode)%ithetamin
      elecs(ielectrode)%ithetamax = elecs(jelectrode)%ithetamax
    enddo
  endif

  !-------------------------------------------------------------------
  ! Calculate top and bottom z-coordinates of the electrodes.
  !-------------------------------------------------------------------
  call calc_electrode_coordinates(par%nrings, par%nel, kguards, elecs)

  !=======================================================
  ! Compute height interval along z-axis, and spacing.
  !=======================================================

! VO VO 13 Jan 2015 Fixed problems of different electrode top and bottom position at different levels.
  igaps_on_z = nint(par%sens%space_elec_guards / (par%sens%heicyl/dble(nz+1)))

  if (igaps_on_z /= 0) then

    ! izspace = the number of points on the gap, i.e., those which will not get any Dirichlet BCs.
    if (par%ifixed_elecgeo == 1) then
      izspace = par%dims%nz/36 / 2**(ilevel-1)
    else
      izspace = igaps_on_z
    endif

    dz1 = (0.25_CUSTOM_REAL*par%sens%heicyl - par%sens%space_elec_guards) / dble(kguards+1-izspace-1)
    dz2 = par%sens%space_elec_guards / dble(izspace+1)
    dz3 = 0.5_CUSTOM_REAL*par%sens%heicyl / dble(3*kguards-kguards-1)
    dz4 = dz2
    dz5 = dz1

    ! lower guard.
    do k=0,kguards+1-izspace-1
      dz(k) = dz1
      z(k)  = dble(k)*dz1
    enddo

    ! lower gap.
    do k=kguards+1-izspace,kguards+1
      dz(k) = dz2
      z(k)  = (0.25_CUSTOM_REAL*par%sens%heicyl - par%sens%space_elec_guards) + dble(k-kguards+izspace)*dz2
    enddo

    ! electrode.
    do k=kguards+2,3*kguards
      dz(k) = dz3
      z(k)  = 0.25_CUSTOM_REAL*par%sens%heicyl + dble(k-kguards-1)*dz3
    enddo

    ! upper gap.
    do k=3*kguards+1,3*kguards+izspace+1
      dz(k) = dz4
      z(k)  = 0.75_CUSTOM_REAL*par%sens%heicyl + dble(k-3*kguards)*dz4
    enddo

    ! upper guard.
    do k=3*kguards+izspace+2,nz+1
      dz(k) = dz5
      z(k)  = (0.75_CUSTOM_REAL*par%sens%heicyl + par%sens%space_elec_guards) + dble(k-3*kguards-izspace-1)*dz5
    enddo
! VO VO 13 Jan 2015 End of the fix.

  else
! VO VO 9 Jan 2015 Fixed problems of different electrode top and bottom position at different levels.
    izspace = 0

    ! The guard height.
    heiguard = 0.25_CUSTOM_REAL*par%sens%heicyl

    dz1 = heiguard/dble(kguards+1)
    dz2 = (par%sens%heicyl-2._CUSTOM_REAL*heiguard) / dble(3*kguards-kguards-1)
    dz3 = dz1

    do k=0,kguards+1
      dz(k) = dz1
      z(k)  = dble(k)*dz1
    enddo

    do k=kguards+2,3*kguards
      dz(k) = dz2
      z(k)  = heiguard + dble(k-kguards-1)*dz2
    enddo

    do k=3*kguards+1,nz+1
      dz(k) = dz3
      z(k)  = (par%sens%heicyl - heiguard) + dble(k-3*kguards)*dz3
    enddo
! VO VO 9 Jan 2015 End of the fix.
  endif

  !------------------------------------------------------------------
  ! Calculate the z grid coordinates for the guards.
  guards%lower_izmax = kguards+1-izspace-1
  guards%upper_izmin = 3*kguards+izspace+1

end subroutine set_geometry

!===========================================================================================
! Compute grid coordinates.
! This grid is used for model initialization and visualization.
!===========================================================================================
subroutine compute_grid_coordinates(nr, ntheta, nz, r, theta, z, xgrid, ygrid, zgrid)
  integer, intent(in) :: nr, ntheta, nz
  real(kind=CUSTOM_REAL), intent(in) :: r(0:)
  real(kind=CUSTOM_REAL), intent(in) :: theta(0:)
  real(kind=CUSTOM_REAL), intent(in) :: z(0:)

  real(kind=CUSTOM_REAL), intent(out) :: xgrid(0:,0:,0:)
  real(kind=CUSTOM_REAL), intent(out) :: ygrid(0:,0:,0:)
  real(kind=CUSTOM_REAL), intent(out) :: zgrid(0:,0:,0:)

  ! Local variables.
  integer :: i,j,k
  real(kind=CUSTOM_REAL), allocatable :: x(:,:)
  real(kind=CUSTOM_REAL), allocatable :: y(:,:)

  allocate(x(0:2*nr+1,0:ntheta+1))
  allocate(y(0:2*nr+1,0:ntheta+1))

  ! x-y-z coordinates of r-theta points
  do j=0,ntheta+1
    do i=0,2*nr+1
      x(i,j) = r(i)*cos(theta(j))
      y(i,j) = r(i)*sin(theta(j))
    enddo
  enddo

  ! x-y-z coordinates of r-theta points 1 over 2 points
  ! DG DG and DK DK: create a cutplane in (r,theta) direction, and extend
  ! it in the z direction
  ! DG DG and DK DK: staggered grid, thus only use corner points
  ! VO VO 6 Jan 2015 fixed according to the new, correct radii indices mapping i-->2*i-1
  do k=0,nz+1
    do j=0,ntheta+1
      xgrid(0,j,k) = 0._CUSTOM_REAL
      ygrid(0,j,k) = 0._CUSTOM_REAL
      zgrid(0,j,k) = z(k)

      do i=1,nr+1
        xgrid(i,j,k) = x(2*i-1,j)
        ygrid(i,j,k) = y(2*i-1,j)
        zgrid(i,j,k) = z(k)
      enddo
    enddo
  enddo

  deallocate(x)
  deallocate(y)

end subroutine compute_grid_coordinates

!===========================================================================================
! Calculate the volume of the cells which is used to calculate the sensitivity or
! to scale the columns of the sensitivity matrix and the model.
!===========================================================================================
subroutine calculate_cell_volumes(dims, volume, volume_avg, r, theta, z, myrank, nbproc)
  type(t_dimensions), intent(in) :: dims
  real(kind=CUSTOM_REAL), intent(in) :: r(0:)
  real(kind=CUSTOM_REAL), intent(in) :: theta(0:)
  real(kind=CUSTOM_REAL), intent(in) :: z(0:)
  integer, intent(in) :: myrank, nbproc

  real(kind=CUSTOM_REAL), intent(out) :: volume(:,:,:)
  real(kind=CUSTOM_REAL), intent(out) :: volume_avg

  integer :: i,j,k,k1
  real(kind=CUSTOM_REAL) :: dz,dtheta
  integer :: n_avg

  volume_avg = 0._CUSTOM_REAL
  n_avg = 0

  do k=1,dims%nzlocal + (myrank+1)/nbproc
    k1 = myrank*dims%nzlocal+k
    dz = z(k1)-z(k1-1)

    do j=1,dims%ntheta
      dtheta = abs(theta(j+1)-theta(j))

      do i=1,dims%nr+1
        if (i == 1) then
          ! Initialize also these elements (not used). Note that r(-1) = 0 by definition.
          volume(i,j,k) = 0.5_CUSTOM_REAL * r(2*i-1)**2 * dtheta * dz
        else
          volume(i,j,k) = 0.5_CUSTOM_REAL * (r(2*i-1)**2 - r(2*i-3)**2) * dtheta * dz
        endif

        volume_avg = volume_avg + volume(i,j,k)
        n_avg = n_avg + 1
      enddo
    enddo
  enddo

  ! TODO: Gather volume_avg for parallel version.
  volume_avg = volume_avg / dble(n_avg)

end subroutine calculate_cell_volumes

!===========================================================================================
! Shift grid coordinates to the nodes where the permittivity is defined.
! This is needed to for example visualize the sensitivity data which is calculated on those nodes.
!===========================================================================================
subroutine compute_shifted_grid_coordinates(nr,ntheta,nz,r,theta,z,xgrid,ygrid,zgrid)
  integer, intent(in) :: nr, ntheta, nz
  real(kind=CUSTOM_REAL), intent(in) :: r(0:)
  real(kind=CUSTOM_REAL), intent(in) :: theta(0:)
  real(kind=CUSTOM_REAL), intent(in) :: z(0:)

  real(kind=CUSTOM_REAL), intent(out) :: xgrid(0:,0:,0:)
  real(kind=CUSTOM_REAL), intent(out) :: ygrid(0:,0:,0:)
  real(kind=CUSTOM_REAL), intent(out) :: zgrid(0:,0:,0:)

  real(kind=CUSTOM_REAL), allocatable :: r_shifted(:)
  real(kind=CUSTOM_REAL), allocatable :: theta_shifted(:)
  real(kind=CUSTOM_REAL), allocatable :: z_shifted(:)
  integer :: i,j,k

  allocate(r_shifted(0:2*nr+1))
  allocate(theta_shifted(0:ntheta+1))
  allocate(z_shifted(0:nz+1))

  ! Shift the radius r.
  ! Permittivity permit(i,:,:) is located between phi(i-1,:,:) and phi(i,:,:).
  do i=1,2*nr+1
    r_shifted(i) = r(i-1)
  enddo
  ! Since we have a mapping i-->2*i-1 and therefore r(-1)=0.
  r_shifted(0) = 0._CUSTOM_REAL

  ! Shift the angle theta.
  ! Permittivity permit(:,j,:) is located between phi(:,j,:) and phi(:,j+1,:).
  do j=1,ntheta
    theta_shifted(j) = theta(j) + 0.5_CUSTOM_REAL*(theta(j+1)-theta(j))
  enddo

  ! Apply periodic BCs.
  theta_shifted(0) = theta_shifted(ntheta) + 2._CUSTOM_REAL*PI
  theta_shifted(ntheta+1) = theta_shifted(1) - 2._CUSTOM_REAL*PI

  ! Shift the height z.
  ! Permittivity permit(:,:,k) is located between phi(:,:,k-1) and phi(:,:,k).
  do k=1,nz+1
    z_shifted(k) = z(k) - 0.5_CUSTOM_REAL*(z(k)-z(k-1))
  enddo
  ! We do not have permittivity defined between z(0)=0 and z(-1).
  ! So just define z_shifted(0) to do not have undefined elements.
  z_shifted(0) = 0._CUSTOM_REAL

  call compute_grid_coordinates(nr,ntheta,nz,r_shifted,theta_shifted,z_shifted,xgrid,ygrid,zgrid)

  deallocate(z_shifted)
  deallocate(theta_shifted)
  deallocate(r_shifted)

end subroutine compute_shifted_grid_coordinates

!===========================================================================================
! Increase twice the theta-dimension for the electrode ielectrode and its gaps.
subroutine refine_geometry(ielectrode,ntheta,theta,nel,elecs,&
                           ntheta0,theta0,elecs0,ithetamap)

  integer, intent(in) :: ielectrode,ntheta,nel,ntheta0
  real(kind=CUSTOM_REAL), intent(inout) :: theta0(0:), theta(0:)
  type(t_electrode), intent(inout) :: elecs(:)
  type(t_electrode), intent(inout) :: elecs0(:)
  integer, intent(inout) :: ithetamap(0:)

  integer, save :: originalValuesStored = 0
  integer, save :: npelec,npgaps

  integer :: j,l,el
  logical :: toInsert

  if (originalValuesStored == 0) then
    theta0(:) = theta(:)
    elecs0(:) = elecs(:)

    npelec = get_npelec(ntheta0, nel)
    npgaps = get_npgaps(ntheta0, nel)

    originalValuesStored = 1
  endif

  l=0
  do j=1,ntheta0+1
    toInsert = .false.

    if (ielectrode > 1 .and. ielectrode < 12) then
      if (j > elecs0(ielectrode-1)%ithetamax .and. j <= elecs0(ielectrode+1)%ithetamin) toInsert = .true.
    else if (ielectrode == 1) then
      if ((j > elecs0(ielectrode)%ithetamin .and. j <= elecs0(ielectrode+1)%ithetamin) &
          .or. (j > elecs0(nel)%ithetamax)) toInsert = .true.
    else if (ielectrode == 12) then
      if (j > elecs0(ielectrode-1)%ithetamax) toInsert = .true.
    endif

    if (toInsert) then
    ! Insert an extra grid dimension between theta(j-1) and theta(j).
      l=l+1
      theta(l) = 0.5_CUSTOM_REAL*(theta0(j-1)+theta0(j))
    endif

    l=l+1
    theta(l) = theta0(j)

    ithetamap(j) = l

    ! Sanity check.
    if (j == ntheta0+1) then
      if (l /= ntheta+1) then
        print *, "Error in mesh refinement! j=", j, ", l=", l
        stop
      endif
    endif
  enddo

  do el=1,nel
    if (el < ielectrode) then
        elecs(el)%ithetamin = elecs0(el)%ithetamin
        elecs(el)%ithetamax = elecs0(el)%ithetamax
    else if (el == ielectrode) then
      if (ielectrode == 1) then
        elecs(el)%ithetamax = elecs0(el)%ithetamax + (npelec-1)
      else  if (ielectrode > 1) then
        elecs(el)%ithetamin = elecs0(el)%ithetamin + npgaps+1
        elecs(el)%ithetamax = elecs0(el)%ithetamax + npgaps + npelec
      endif
    else if (el > ielectrode) then
      if (ielectrode == 1) then
        elecs(el)%ithetamin = elecs0(el)%ithetamin + npgaps + npelec
        elecs(el)%ithetamax = elecs0(el)%ithetamax + npgaps + npelec
      else
        elecs(el)%ithetamin = elecs0(el)%ithetamin + 2*(npgaps+1) + (npelec-1)
        elecs(el)%ithetamax = elecs0(el)%ithetamax + 2*(npgaps+1) + (npelec-1)
      endif
    endif
  enddo

  ! Periodic boundaries.
  theta(0) = theta(ntheta) + 2.0_CUSTOM_REAL*PI
  ithetamap(0) = ithetamap(ntheta0)

end subroutine refine_geometry

end module geometry

max_stars_repo_name: RichardScottOZ/Tomofast-x

max_issues_repo_name: RichardScottOZ/Tomofast-x

max_forks_repo_name: RichardScottOZ/Tomofast-x


text:
*DECK DBESK0
      DOUBLE PRECISION FUNCTION DBESK0 (X)
C***BEGIN PROLOGUE  DBESK0
C***PURPOSE  Compute the modified (hyperbolic) Bessel function of the
C            third kind of order zero.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C10B1
C***TYPE      DOUBLE PRECISION (BESK0-S, DBESK0-D)
C***KEYWORDS  FNLIB, HYPERBOLIC BESSEL FUNCTION,
C             MODIFIED BESSEL FUNCTION, ORDER ZERO, SPECIAL FUNCTIONS,
C             THIRD KIND
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DBESK0(X) calculates the double precision modified (hyperbolic)
C Bessel function of the third kind of order zero for double
C precision argument X.  The argument must be greater than zero
C but not so large that the result underflows.
C
C Series for BK0        on the interval  0.          to  4.00000E+00
C                                        with weighted error   3.08E-33
C                                         log weighted error  32.51
C                               significant figures required  32.05
C                                    decimal places required  33.11
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DBESI0, DBSK0E, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DBESK0
      DOUBLE PRECISION X, BK0CS(16), XMAX, XMAXT, XSML, Y,
     1  D1MACH, DCSEVL, DBESI0, DBSK0E
      LOGICAL FIRST
      SAVE BK0CS, NTK0, XSML, XMAX, FIRST
      DATA BK0CS(  1) / -.3532739323 3902768720 1140060063 153 D-1    /
      DATA BK0CS(  2) / +.3442898999 2462848688 6344927529 213 D+0    /
      DATA BK0CS(  3) / +.3597993651 5361501626 5721303687 231 D-1    /
      DATA BK0CS(  4) / +.1264615411 4469259233 8479508673 447 D-2    /
      DATA BK0CS(  5) / +.2286212103 1194517860 8269830297 585 D-4    /
      DATA BK0CS(  6) / +.2534791079 0261494573 0790013428 354 D-6    /
      DATA BK0CS(  7) / +.1904516377 2202088589 7214059381 366 D-8    /
      DATA BK0CS(  8) / +.1034969525 7633624585 1008317853 089 D-10   /
      DATA BK0CS(  9) / +.4259816142 7910825765 2445327170 133 D-13   /
      DATA BK0CS( 10) / +.1374465435 8807508969 4238325440 000 D-15   /
      DATA BK0CS( 11) / +.3570896528 5083735909 9688597333 333 D-18   /
      DATA BK0CS( 12) / +.7631643660 1164373766 7498666666 666 D-21   /
      DATA BK0CS( 13) / +.1365424988 4407818590 8053333333 333 D-23   /
      DATA BK0CS( 14) / +.2075275266 9066680831 9999999999 999 D-26   /
      DATA BK0CS( 15) / +.2712814218 0729856000 0000000000 000 D-29   /
      DATA BK0CS( 16) / +.3082593887 9146666666 6666666666 666 D-32   /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DBESK0
      IF (FIRST) THEN
         NTK0 = INITDS (BK0CS, 16, 0.1*REAL(D1MACH(3)))
         XSML = SQRT(4.0D0*D1MACH(3))
         XMAXT = -LOG(D1MACH(1))
         XMAX = XMAXT - 0.5D0*XMAXT*LOG(XMAXT)/(XMAXT+0.5D0)
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'DBESK0',
     +   'X IS ZERO OR NEGATIVE', 2, 2)
      IF (X.GT.2.0D0) GO TO 20
C
      Y = 0.D0
      IF (X.GT.XSML) Y = X*X
      DBESK0 = -LOG(0.5D0*X)*DBESI0(X) - 0.25D0 + DCSEVL (.5D0*Y-1.D0,
     1  BK0CS, NTK0)
      RETURN
C
 20   DBESK0 = 0.D0
      IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DBESK0',
     +   'X SO BIG K0 UNDERFLOWS', 1, 1)
      IF (X.GT.XMAX) RETURN
C
      DBESK0 = EXP(-X) * DBSK0E(X)
C
      RETURN
      END

max_stars_repo_name: ygeorgi/MESS

max_issues_repo_name: ygeorgi/MESS

max_forks_repo_name: ygeorgi/MESS


text:
C 	Copyright restrictions apply - see stsdas$copyright.stsdas
C
        SUBROUTINE YRCCS9(CCS9,DET,FGWA,POS,PEDGRE,DESCRP,ISTAT)
*
*  Module number:
*
*  Module name: YRCCS9
*
*  Keyphrase:
*  ----------
*       read scattered light parameter table
*
*  Description:
*  ------------
*       This routine reads the table containing the beginning and
*       ending diode numbers of the regions of no sensitivity
*       from which the contribution of scattered light is measured.
*       Rows are selected by detector and fgwa_id.
*
*  FORTRAN name: YRCCS9.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       CCS9                    I       table containing CCS9 coeffecients
*
*  Subroutines Called:
*  -------------------
*  CDBS
*       ymsput
*  SDAS:
*       uttopn, utpgti, utcfnd, utrgt*, uttclo
*
*  History:
*  --------
*  Version      Date        Author          Description
*	1	Mar 94	H. Bushouse	Designed and coded
*	2	Mar 94	H. Bushouse	Mods to handle PEDIGREE keywords
*-------------------------------------------------------------------------------
*
* Input parameters
*
*       CCS9 - table name
*       DET  - detector name ('AMBER' or 'BLUE')
*       FGWA - grating mode
*
* Output parameters
*       POS - first and last diode values of no sensitivity range
*	PEDGRE - CCS9 PEDIGREE keyword
*	DESCRP - CCS9 DESCRIP  keyword
*       ISTAT - error status
*
*-----------------------------------------------------------------------------
C
C     INPUT/OUTPUT DECLARATIONS
C
	CHARACTER*64 CCS9
	CHARACTER*5 DET
	CHARACTER*3 FGWA
	CHARACTER*68 PEDGRE, DESCRP
	REAL POS(2)
	INTEGER ISTAT
C
C     BEGIN IRAF77.INC
C     FILE I/O ACCESS MODES
C
	INTEGER   RDONLY
	PARAMETER (RDONLY = 1)
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
	INTEGER STDOUT
	PARAMETER (STDOUT = 1)
	INTEGER STDERR
	PARAMETER (STDERR = 2)
C
C     THESE MAY BE READ BY UTPGTI BUT MAY NOT BE SET:
C
C                                       NUMBER OF ROWS WRITTEN TO
	INTEGER   TBNROW
	PARAMETER (TBNROW = 21)
C     END IRAF77.INC
C
C local variables
C
	INTEGER IDIN,COLIDS(4),ROW,NROWS
	CHARACTER*80 CONTXT
	CHARACTER*15 COLNAM(4)
	CHARACTER*5 DET1
	CHARACTER*3 FGWA1
	LOGICAL NULL(2)
	DATA COLNAM/'DETECTOR','FGWA_ID','RANGE_BEG','RANGE_END'/
C---------------------------------------------------------------------------
C
C Initialize range to zero
C
	POS(1) = 0.0
	POS(2) = 0.0
C
C Open table
C
        CALL UTTOPN(CCS9,RDONLY,IDIN,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR opening CCS9 table '//CCS9
                GO TO 998
        ENDIF
C
C get PEDIGREE and DESCRIP header keywords
C
        CALL UTHGTT(IDIN,'PEDIGREE',PEDGRE,ISTAT)
        CALL UTHGTT(IDIN,'DESCRIP',DESCRP,ISTAT)
C
C get number of rows
C
        CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading CCS9 table '//CCS9
                GO TO 999
        ENDIF
C
C Get column ids.
C
        CALL UTCFND(IDIN,COLNAM,4,COLIDS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR locating needed columns in CCS9 table '//
     *                  CCS9
                GO TO 999
        ENDIF
C
C Loop on rows
C
        DO 500 ROW=1,NROWS
C
C check detector
C
                CONTXT='ERROR reading CCS9 table '//CCS9
                CALL UTRGTT(IDIN,COLIDS(1),1,ROW,DET1,NULL,ISTAT)
                IF(ISTAT.NE.0)GO TO 999
                IF(DET1.NE.DET) GO TO 500
C
C Check fgwa_id
C
                CALL UTRGTT(IDIN,COLIDS(2),1,ROW,FGWA1,NULL,ISTAT)
                IF(ISTAT.NE.0)GO TO 999
                IF(FGWA1.NE.FGWA) GO TO 500
C
C read limits
C
                CALL UTRGTR(IDIN,COLIDS(3),2,ROW,POS,NULL,ISTAT)
                IF(ISTAT.NE.0) GO TO 999
		GO TO 600
500	CONTINUE
C
C If we got here no applicable row was found
C
	WRITE(CONTXT,98)DET,FGWA
98	FORMAT('ERROR: No row found in CCS9 for ',A5,1X,A3)
	GO TO 999
600	CONTINUE
C
C If we got here we found a row
C
        CALL UTTCLO(IDIN,ISTAT)
        ISTAT=0
        GO TO 1000
999     CALL UTTCLO(IDIN,ISTAT)
998     CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
Describe Users/LanceStark here.

20100823 13:51:59 nbsp Hi Lance, and welcome to the wiki!  Its always great seeing business owners/representatives conducting good customer service over the wiki.  Thanks for reaching out! Users/TomGarberson


text:
      MODULE DE4MODULE
        INTEGER,SAVE,POINTER  ::MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
        INTEGER,SAVE,POINTER  ::NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM
        INTEGER,SAVE,POINTER  ::NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
        REAL   ,SAVE,POINTER  ::ACCLDE4,HCLOSEDE4,DELTL
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::IUPPNT
        INTEGER,          SAVE, POINTER, DIMENSION(:,:,:) ::IEQPNT
        REAL,             SAVE, POINTER, DIMENSION(:,:)   ::AU
        REAL,             SAVE, POINTER, DIMENSION(:,:)   ::AL
        REAL,             SAVE, POINTER, DIMENSION(:)     ::D4B
        REAL,             SAVE, POINTER, DIMENSION(:)     ::HDCGDE4
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::LRCHDE4
      TYPE DE4TYPE
        INTEGER,POINTER  ::MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
        INTEGER,POINTER  ::NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM
        INTEGER,POINTER  ::NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
        REAL   ,POINTER  ::ACCLDE4,HCLOSEDE4,DELTL
        INTEGER,           POINTER, DIMENSION(:,:)   ::IUPPNT
        INTEGER,           POINTER, DIMENSION(:,:,:) ::IEQPNT
        REAL,              POINTER, DIMENSION(:,:)   ::AU
        REAL,              POINTER, DIMENSION(:,:)   ::AL
        REAL,              POINTER, DIMENSION(:)     ::D4B
        REAL,              POINTER, DIMENSION(:)     ::HDCGDE4
        INTEGER,           POINTER, DIMENSION(:,:)   ::LRCHDE4
      END TYPE
      TYPE(DE4TYPE), SAVE ::DE4DAT(10)
      END MODULE DE4MODULE

      SUBROUTINE DE47AR(INDE4,MXITER,IGRID)
C     ******************************************************************
C     ALLOCATE STORAGE IN X ARRAY FOR D4 ARRAYS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,   ONLY:IOUT,NCOL,NROW,NLAY
      USE DE4MODULE
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
      ALLOCATE(MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR)
      ALLOCATE(NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM)
      ALLOCATE(NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW)
      ALLOCATE(ACCLDE4,HCLOSEDE4,DELTL)
C
C1------PRINT A MESSAGE IDENTIFYING DE4 PACKAGE.
      WRITE(IOUT,1) INDE4
    1 FORMAT(1X,/1X,'DE4 -- D4 DIRECT SOLUTION PACKAGE',
     1    ', VERSION 7, 5/2/2005  INPUT READ FROM UNIT ',I4)
C
C2------CALCULATE DEFAULT VALUES FOR MXUP, MXLOW, AND MXBW.
C2------ALSO SET VALUES FOR ID4DIR (DIRECTION OF EQUATION ORDERING) AND
C2------ID4DIM (MAXIMUM NUMBER OF HEAD COEFFICIENTS FOR ONE EQUATION).
C2------ID4DIM=5 for a 2-D problem; ID4DIM=7 for 3-D.
      NODES=NCOL*NROW*NLAY
      NHALFU=(NODES-1)/2 +1
      NHALFL=NODES-NHALFU
      ID4DIM=7
      DELTL = 0.
      NBWL = 0
      NUPL = 0
      NLOWL = 0
      IF(NLAY.LE.NCOL .AND. NLAY.LE.NROW) THEN
         IF(NLAY.EQ.1) ID4DIM=5
         IF(NCOL.GE.NROW) THEN
            ID4DIR=1
            NBWGRD=NROW*NLAY+1
         ELSE
            ID4DIR=2
            NBWGRD=NCOL*NLAY+1
         END IF
      ELSE IF(NROW.LE.NCOL .AND. NROW.LE.NLAY) THEN
         IF(NROW.EQ.1) ID4DIM=5
         IF(NCOL.GE.NLAY) THEN
            ID4DIR=3
            NBWGRD=NROW*NLAY+1
         ELSE
            ID4DIR=4
            NBWGRD=NROW*NCOL+1
         END IF
      ELSE
         IF(NCOL.EQ.1) ID4DIM=5
         IF(NROW.GE.NLAY) THEN
            ID4DIR=5
            NBWGRD=NCOL*NLAY+1
         ELSE
            ID4DIR=6
            NBWGRD=NCOL*NROW+1
         END IF
      END IF
C
C3------READ AND PRINT COMMENTS, ITMX, MXUP, MXLOW, MXBW.  FOR ANY
C3------ZERO OR NEGATIVE VALUES, SUBSTITUE THE DEFAULT VALUE.
      CALL URDCOM(INDE4,IOUT,LINE)
      LLOC=1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ITMX,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXUP,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXLOW,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXBW,R,IOUT,INDE4)
      IF(ITMX.LT.1) ITMX=1
      WRITE(IOUT,3) ITMX
    3 FORMAT(1X,'MAXIMUM ITERATIONS (EXTERNAL OR INTERNAL) =',I3)
      IF(MXUP.LT.1) MXUP=NHALFU
      IF(MXLOW.LT.1) MXLOW=NHALFL
      MXEQ=MXUP+MXLOW
      IF(MXBW.LT.1) MXBW=NBWGRD
      WRITE(IOUT,4) MXUP,MXLOW,MXBW
    4 FORMAT(1X,'MAXIMUM EQUATIONS IN UPPER PART OF [A]:',I7,/
     1       1X,'MAXIMUM EQUATIONS IN LOWER PART OF [A]:',I7,/
     2       1X,'MAXIMUM BAND WIDTH OF [AL] PLUS 1:',I5)
C
C4------ALLOCATE SPACE FOR THE DE4 ARRAYS.
      ALLOCATE (AU(ID4DIM,MXUP))
      ALLOCATE (IUPPNT(ID4DIM,MXUP))
      ALLOCATE (AL(MXBW,MXLOW))
      ALLOCATE (IEQPNT(NCOL,NROW,NLAY))
      ALLOCATE (D4B(MXEQ))
      ALLOCATE (LRCHDE4(3,ITMX))
      ALLOCATE (HDCGDE4(ITMX))
C
C1------READ IFREQ, MUTD4, ACCL, HCLOSE, AND IPRD4
      ONE=1.
      ZERO=0.
      READ(INDE4,*) IFREQ,MUTD4,ACCLDE4,HCLOSEDE4,IPRD4
      IF(ACCLDE4.LE.ZERO) ACCLDE4=ONE
      IF(IPRD4.LE.0) IPRD4=999
      IF(MUTD4.LT.0 .OR. MUTD4.GT.3) MUTD4=0
      IF(IFREQ.LT.1 .OR. IFREQ.GT.3) THEN
         WRITE(IOUT,11) IFREQ
   11    FORMAT(1X,/1X,'INVALID VALUE FOR IFREQ PARAMETER:',I8)
         CALL USTOP(' ')
      END IF
C
C2------CHECK TO SEE IF THERE IS ITERATION (ITMX>1).
      IF(ITMX.GT.1) THEN
C
C3------THERE IS ITERATION -- DETERMINE TYPE OF ITERATION BASED ON IFREG
C3------VALUE.
         IF(IFREQ.EQ.3) THEN
C
C3A-----EXTERNAL ITERATION -- SET ITERATION VARIABLES AND PRINT A
C3A-----MESSAGE.
            MXITER=ITMX
            NITERDE4=1
            WRITE(IOUT,51)
   51       FORMAT(1X,/21X,'SOLUTION BY D4 DIRECT SOLVER WITH EXTERNAL',
     1       ' ITERATION',/21X,52('-'),/)
C
C3B-----INTERNAL ITERATION -- SET ITERATION VARIABLES AND PRINT A
C3B-----MESSAGE.
         ELSE
            MXITER=1
            NITERDE4=ITMX
            WRITE(IOUT,81)
   81       FORMAT(1X,/21X,'SOLUTION BY D4 DIRECT SOLVER WITH INTERNAL',
     1       ' ITERATION'/21X,52('-')/)
         END IF
C
C3C-----PRINT ITERATION INFORMATION.
         WRITE(IOUT,91) ITMX,ACCLDE4,HCLOSEDE4,IPRD4
   91    FORMAT(1X,'MAXIMUM ITERATIONS =',I3/
     1       1X,'RELAXATION-ACCELERATION PARAMETER =',F10.6/
     2       1X,'HEAD CHANGE CRITERION FOR CLOSURE =',G15.5/
     3       1X,'D4 PRINTOUT INTERVAL =',I4)
         IF(MUTD4.EQ.1) WRITE(IOUT,92)
   92    FORMAT(1X,
     1   'CONVERGENCE PRINTOUT WILL SHOW ONLY THE NUMBER OF ITERATIONS')
         IF(MUTD4.EQ.2) WRITE(IOUT,93)
   93    FORMAT(1X,'CONVERGENCE PRINTOUT WILL BE SUPPRESSED')
C
C4------NO ITERATION -- SET ITERATION VARIABLES AND PRINT A MESSAGE.
      ELSE
         MXITER=1
         NITERDE4=1
         ACCLDE4=ONE
         WRITE(IOUT,94)
   94    FORMAT(1X,/21X,
     1   'SOLUTION BY D4 DIRECT SOLVER WITH NO ITERATION',/21X,46('-')/)
         IF(MUTD4.EQ.2) WRITE(IOUT,95)
   95    FORMAT(1X,'PRINTOUT OF MAXIMUM HEAD CHANGE WILL BE SUPPRESSED')
      END IF
C
C5------PRINT MESSAGE ABOUT FREQUENCY AT WHICH [A] IS ELIMINATED.
      IF(IFREQ.EQ.3) WRITE(IOUT,102)
  102 FORMAT(1X,'NON-LINEAR PROBLEM -- [A] MATRIX ELIMINATED',
     1   ' EVERY TIME EQUATIONS ARE REFORMULATED')
      IF(IFREQ.NE.3) WRITE(IOUT,103) IFREQ
  103 FORMAT(1X,'LINEAR PROBLEM -- [A] MATRIX ELIMINATED ONLY',
     1   ' WHEN IT CHANGES -- IFREQ=',I1)
C
C6------RETURN.
      CALL DE47PSV(IGRID)
      RETURN
      END
      SUBROUTINE DE47AP(HNEW,IBOUND,AU,AL,IUPPNT,IEQPNT,D4B,MXUP,MXLOW,
     1  MXEQ,MXBW,CR,CC,CV,HCOF,RHS,ACCLDE4,KITER,ITMX,MXITER,NITERDE4,
     2  HCLOSEDE4,IPRD4,ICNVG,NCOL,NROW,NLAY,IOUT,LRCHDE4,HDCGDE4,
     3  IFREQ,KSTP,KPER,DELT,NSTP,ID4DIR,ID4DIM,MUTD4,DELTL,NBWL,NUPL,
     4  NLOWL,NLOW,NEQ,NUP,NBW,IERR)
C     ******************************************************************
C     SOLVE FINITE-DIFFERENCE EQUATIONS FOR ONE EXTERNAL ITERATION.
C     MULTIPLE SOLUTIONS ARE MADE WHEN INTERNALLY ITERATING.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION HNEW(NCOL,NROW,NLAY),AU(ID4DIM,MXUP),AL(MXBW,MXLOW),
     1          IEQPNT(NCOL,NROW,NLAY),IUPPNT(ID4DIM,MXUP),D4B(MXEQ),
     2          CR(NCOL,NROW,NLAY),CC(NCOL,NROW,NLAY),
     3          CV(NCOL,NROW,NLAY),HCOF(NCOL,NROW,NLAY),
     4          RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     5          LRCHDE4(3,ITMX),HDCGDE4(ITMX)
C
      DOUBLE PRECISION HNEW,EE,COND,RR,DDH,RHS,HCOF
C
      DIMENSION CND(6),IEQ(6),IDIR(6,6)
C      SAVE DELTL,NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
      DATA IDIR/1,2,3,4,5,6,
     1          2,1,3,4,6,5,
     2          1,3,2,5,4,6,
     3          3,1,2,5,6,4,
     4          2,3,1,6,4,5,
     5          3,2,1,6,5,4/
C
C      DATA DELTL,NBWL,NUPL,NLOWL/0.,0,0,0/
C     ------------------------------------------------------------------
C
C1------INITIALIZE VARIABLES AND SET FLAG THAT INDICATES IF [A] REQUIRES
C1------ELIMINATION.
      ZERO=0.
      ICNVG=0
      NIT=0
      ITYPE=0
      IF(IFREQ.EQ.0) THEN
         ITYPE=1
      ELSE IF(IFREQ.EQ.2) THEN
         IF(KSTP.NE.1 .AND. DELT.EQ.DELTL) ITYPE=1
      ELSE IF(IFREQ.EQ.1) THEN
         IF((KPER.NE.1 .OR. KSTP.NE.1) .AND. DELT.EQ.DELTL) ITYPE=1
      END IF
      DELTL=DELT
C
C2------DO ONE INTERNAL ITERATION.
   10 NIT=NIT+1
      BIGA=ZERO
      BIG=ZERO
      IBIG=0
      JBIG=0
      KBIG=0
C
C3------GO TO STATEMENT 100 IF [A] DOES NOT REQUIRE ELIMINATION.
      IF (ITYPE.EQ.1) GO TO 100
C
C4------[A] REQUIRES ELIMINATION.
C4A-----CALL MODULE SDE47N TO ORDER EQUATIONS.  IF MXUP OR MXLOW ARE
C4A-----TOO SMALL, PRINT A MESSAGE AND STOP.
      CALL SDE47N(IEQPNT,IBOUND,NCOL,NROW,NLAY,ID4DIR,NUP,NLOW,NEQ)
      IF(NUP.GT.MXUP .OR. NLOW.GT.MXLOW) THEN
         WRITE(IOUT,41) NUP,NLOW
   41    FORMAT(1X,'INSUFFICIENT MEMORY FOR DE4 SOLVER:',/
     1          1X,'MXUP MUST BE AT LEAST',I8,/
     1          1X,'MXLOW MUST BE AT LEAST',I8)
         IERR = 1
         RETURN
      END IF
C
C4B-----INITIALIZE AU.
      DO 50 I=1,NUP
      DO 50 J=2,ID4DIM
      AU(J,I)=ZERO
   50 CONTINUE
C
C5------LOOP THROUGH ALL CELLS CALCULATING COEFFICIENTS AND LOADING
C5------ARRAYS FOR SOLUTION.  THE RESIDUAL MUST ALWAYS BE CALCULATED
C5------AND LOADED IN D4B; IUPPNT, AU, AND AL(1,n) ARE CALCULATED AND
C5------LOADED ONLY IF [A] REQUIRES ELIMINATION.
  100 DO 310 K=1,NLAY
      DO 310 I=1,NROW
      DO 310 J=1,NCOL
      IR=IEQPNT(J,I,K)
      IF(IR.EQ.0) GO TO 310
C
C5A-----CALCULATE AND LOAD D4B.
      DO 110 N=1,6
      IEQ(N)=0
  110 CONTINUE
      RR=RHS(J,I,K)
      EE=HCOF(J,I,K)
      IF(J.NE.1) THEN
         CLF=CR(J-1,I,K)
         COND=CLF
         RR=RR-COND*HNEW(J-1,I,K)
         EE=EE-COND
         CND(1)=CLF
         IEQ(1)=IEQPNT(J-1,I,K)
      END IF
      IF(I.NE.1) THEN
         CBK=CC(J,I-1,K)
         COND=CBK
         RR=RR-COND*HNEW(J,I-1,K)
         EE=EE-COND
         CND(2)=CBK
         IEQ(2)=IEQPNT(J,I-1,K)
      END IF
      IF(K.NE.1) THEN
         CUP=CV(J,I,K-1)
         COND=CUP
         RR=RR-COND*HNEW(J,I,K-1)
         EE=EE-COND
         CND(3)=CUP
         IEQ(3)=IEQPNT(J,I,K-1)
      END IF
      IF(K.NE.NLAY) THEN
         CDN=CV(J,I,K)
         COND=CDN
         RR=RR-COND*HNEW(J,I,K+1)
         EE=EE-COND
         CND(4)=CDN
         IEQ(4)=IEQPNT(J,I,K+1)
      END IF
      IF(I.NE.NROW) THEN
         CFR=CC(J,I,K)
         COND=CFR
         RR=RR-COND*HNEW(J,I+1,K)
         EE=EE-COND
         CND(5)=CFR
         IEQ(5)=IEQPNT(J,I+1,K)
      END IF
      IF(J.NE.NCOL) THEN
         CRT=CR(J,I,K)
         COND=CRT
         RR=RR-COND*HNEW(J+1,I,K)
         EE=EE-COND
         CND(6)=CRT
         IEQ(6)=IEQPNT(J+1,I,K)
      END IF
      D4B(IR)=RR-EE*HNEW(J,I,K)
C
C5B-----IF [A] REQUIRES ELIMINATION, LOAD AL(1,N) AND AU
      IF (ITYPE.EQ.1) GO TO 310
      IF (IR.GT.NUP) THEN
         IRR=IR-NUP
         AL(1,IRR)=EE
      ELSE
         N=1
         DO 305 II=1,6
            M=IDIR(II,ID4DIR)
            L=IEQ(M)
            IF(L.NE.0) THEN
               N=N+1
               IUPPNT(N,IR)=L
               AU(N,IR)=CND(M)
            END IF
  305    CONTINUE
         AU(1,IR)=EE
         IUPPNT(1,IR)=N
      END IF
  310 CONTINUE
C
C6------IF [A] DOES NOT REQUIRE ELIMINATION, SKIP TO STATEMENT 380.
      IF(ITYPE.EQ.1) GO TO 380
C6A-----[A] REQUIRES ELIMINATION -- DETERMINE BAND WIDTH + 1.
      MNN=999999
      MXN=0
      DO 350 I=1,NUP
      L=IUPPNT(1,I)
      IF(L.LT.2) GO TO 350
      N=IUPPNT(2,I)-I
      IF(N.LT.MNN) MNN=N
      N=IUPPNT(L,I)-I
      IF(N.GT.MXN) MXN=N
  350 CONTINUE
      NBW=MXN-MNN+1
C
C6B-----WRITE BAND WIDTH + 1 AND NUMBER OF EQUATIONS IF ANY HAVE CHANGED.
      IF(NUP.NE.NUPL .OR. NLOW.NE.NLOWL .OR. NBW.NE.NBWL) THEN
         WRITE(IOUT,351) NUP,NLOW,NBW
  351    FORMAT(1X,/1X,I7,' UPPER PART EQS.',
     1     I10,' LOWER PART EQS.    BAND WIDTH + 1 =',I5)
         NUPL=NUP
         NLOWL=NLOW
         NBWL=NBW
      END IF
C
C6C-------STOP IF BAND WIDTH EXCEEDS USER-SPECIFIED SIZE.
      IF(NBW.GT.MXBW) THEN
         WRITE(IOUT,353) NBW
  353    FORMAT(1X,'INSUFFICIENT MEMORY FOR DE4 SOLVER:',/
     1          1X,'MXBW MUST BE AT LEAST',I5)
         IERR = 1
         RETURN
      END IF
C
C6D-------INITIALIZE OFF DIAGONAL PART OF AL.
      DO 360 I=1,NLOW
      DO 360 J=2,NBW
      AL(J,I)=ZERO
  360 CONTINUE
C
C7------CALL MODULE SDE47S TO SOLVE EQUATIONS FOR HEAD CHANGE.
  380 CALL SDE47S(AU,AL,IUPPNT,D4B,NUP,NLOW,NEQ,MXBW,NBW,ITYPE,ID4DIM)
C
C8------CALCULATE NEW HEAD FROM HEAD CHANGE AND FIND MAXIMUM CHANGE.
      DO 400 K=1,NLAY
      DO 400 I=1,NROW
      DO 400 J=1,NCOL
      L=IEQPNT(J,I,K)
      IF(L.EQ.0) GO TO 400
      DH=ACCLDE4*D4B(L)
      TCHK=ABS(DH)
      IF(TCHK.GT.BIGA) THEN
         BIGA=TCHK
         BIG=DH
         IBIG=I
         JBIG=J
         KBIG=K
      END IF
      DDH=DH
      HNEW(J,I,K)=HNEW(J,I,K)+DDH
  400 CONTINUE
C
C9------IF THE NUMBER OF INTERNAL ITERATIONS IS 1, GO TO STATEMENT 500
      IF(NITERDE4.EQ.1) GO TO 500
C
C10-----THE NUMBER OF INTERNAL ITERATIONS IS GREATER THAN 1, SO MUST BE
C10-----ITERATING INTERNALLY.  KEEP TRACK OF MAXIMUM HEAD CHANGE AND
C10-----CHECK FOR CONVERGENCE.
      LRCHDE4(1,NIT)=KBIG
      LRCHDE4(2,NIT)=IBIG
      LRCHDE4(3,NIT)=JBIG
      HDCGDE4(NIT)=BIG
      IF(ABS(BIG).LE.HCLOSEDE4) ICNVG=1
C
C10A----IF NOT CONVERGED AND MAXIMUM ITERATIONS HAS NOT BEEN REACHED,
C10A----GO BACK AND DO ANOTHER INTERNAL ITERATION.  SET ITYPE=1 TO
C10A----AVOID REFORMULATION OF [A].
      ITYPE=1
      IF(ICNVG.EQ.0 .AND. NIT.NE.NITERDE4) GO TO 10
C
C10B----INTERNAL ITERATION IS DONE EITHER BECAUSE CONVERGENCE IS REACHED
C10B----OR BECAUSE MAX. ITERATIONS EXCEEDED.  PRINT CONVERGENCE
C10B----INFORMATION UNLESS IMUTD4=2.
      IF(MUTD4.LT.2) THEN
C10B1---PRINT A BLANK LINE IF THIS IS THE FIRST TIME STEP.
         IF(KSTP.EQ.1) WRITE(IOUT,601)
C10B2---PRINT NUMBER OF ITERATIONS.
         WRITE(IOUT,751) NIT,KSTP,KPER
  751    FORMAT(1X,I5,' INTERNAL ITERATIONS FOR TIME STEP',I6, !gsf
     1        ' IN STRESS PERIOD ',I4)
C10B3---IF MUTD4=0 AND
C10B3---IF FAILED TO CONVERGE OR LAST TIME STEP IN STRESS PERIOD OR
C10B3---IPRD4 INTERVAL IS MET, CALL MODULE SDE47P TO PRINT HEAD CHANGE.
         IF((MUTD4.EQ.0) .AND.
     1      (ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. MOD(KSTP,IPRD4).EQ.0))
     2         CALL SDE47P(HDCGDE4,LRCHDE4,NIT,IOUT,NCOL,NROW)
      ELSE IF(MUTD4.EQ.3 .AND. ICNVG.EQ.0) THEN
         CALL SDE47P(HDCGDE4,LRCHDE4,NIT,IOUT,NCOL,NROW)
      END IF
C10B4---RETURN.
      RETURN
C
C11-----THERE ARE NO INTERNAL ITERATIONS, SO THERE MUST EITHER BE
C11-----EXTERNAL ITERATION OR NO ITERATION.  IF THERE IS EXTERNAL
C11-----ITERATION, GO TO STATEMENT 600.
  500 IF(MXITER.NE.1) GO TO 600
C
C12-----NO ITERATION (NITERDE4=1 AND MXITER=1).  SET FLAG TO INDICATE
C12-----CONVERGENCE HAS OCCURRED, AND PRINT MAXIMUM HEAD CHANGE UNLESS
C12-----MUTD4<2.
      ICNVG=1
      IF(MUTD4.LT.2) WRITE(IOUT,501) KSTP,KPER,BIG,KBIG,IBIG,JBIG
  501 FORMAT(1X,/1X,'MAXIMUM HEAD CHANGE IN TIME STEP',I6, !gsf
     1      ' OF STRESS PERIOD ',I4,' =',G15.5,
     2      ' AT LAYER =',I3,', ROW =',I5,', COL=',I5)
      RETURN
C
C13-----EXTERNAL ITERATION.  KEEP TRACK OF MAXIMUM HEAD CHANGE AND SET
C13-----CONVERGENCE FLAG IF CONVERGENCE OCCURRED.
  600 LRCHDE4(1,KITER)=KBIG
      LRCHDE4(2,KITER)=IBIG
      LRCHDE4(3,KITER)=JBIG
      HDCGDE4(KITER)=BIG
      IF(ABS(BIG).LE.HCLOSEDE4) ICNVG=1
C13A----RETURN IF NO CONVERGENCE AND MAX. ITERATIONS NOT EXCEEDED.
      IF(ICNVG.EQ.0 .AND. KITER.NE.MXITER) RETURN
C13B----EXTERNAL ITERATION IS DONE EITHER BECAUSE CONVERGENCE IS REACHED
C13B----OR BECAUSE MAX. ITERATIONS EXCEEDED.  PRINT CONVERGENCE
C13B----INFORMATION UNLESS IMUTD4=2.
      IF(MUTD4.LT.2) THEN
C13B1---PRINT A BLANK LINE IF THIS IS THE FIRST TIME STEP.
         IF(KSTP.EQ.1) WRITE(IOUT,601)
  601    FORMAT(1X)
C13B2---PRINT NUMBER OF ITERATIONS.
         WRITE(IOUT,602) KITER,KSTP,KPER
  602    FORMAT(1X,I5,' EXTERNAL ITERATIONS FOR TIME STEP',I6, !gsf
     1        ' IN STRESS PERIOD ',I4)
C13B3---IF MUTD4=0 AND
C13B3---IF FAILED TO CONVERGE OR LAST TIME STEP IN STRESS PERIOD OR
C13B3---IPRD4 INTERVAL IS MET, CALL MODULE SDE47P TO PRINT HEAD CHANGE.
         IF((MUTD4.EQ.0) .AND.
     1      (ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. MOD(KSTP,IPRD4).EQ.0))
     2         CALL SDE47P(HDCGDE4,LRCHDE4,KITER,IOUT,NCOL,NROW)
      ELSE IF(MUTD4.EQ.3 .AND. ICNVG.EQ.0) THEN
         CALL SDE47P(HDCGDE4,LRCHDE4,KITER,IOUT,NCOL,NROW)
      END IF
C13B4---Return.
      RETURN
C
      END
      SUBROUTINE SDE47N(IEQPNT,IBOUND,NCOL,NROW,NLAY,ID4DIR,NUP,NLOW,
     1          NEQ)
C     ******************************************************************
C     ORDER EQUATIONS USING D4 ORDERING
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION IEQPNT(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY)
C     ------------------------------------------------------------------
C
C1------CALCULATE MAXIMUM PLANE NUMBER AND INITIALIZE EQUATION POINTERS.
      NPLANE=NCOL+NROW+NLAY
      DO 20 K=1,NLAY
      DO 20 I=1,NROW
      DO 20 J=1,NCOL
      IEQPNT(J,I,K)=0
   20 CONTINUE
      NEQ=0
C
C2------ORDER EQUATIONS BASED ON DIRECTION FLAG, ID4DIR
C2------Ordering is done as described by Price, H.S. and Coats, K.H.,
C2------1974,Direct methods in reservoir simulation: Soc. Petrol. Eng.
C2------Jour., June 1974, p. 295-308.
      GO TO (100,200,300,400,500,600) ID4DIR
      CALL USTOP(' ')
C
C3------DIRECTION 1 -- NCOL>or=NROW>or=NLAY
C3A-----Order equations with odd plane numbers.
  100 DO 130 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NLAY) K1=NLAY
      K2=N-NCOL-NROW
      IF(K2.LT.1) K2=1
      DO 130 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NROW) I1=NROW
      I2=N-K-NCOL
      IF(I2.LT.1) I2=1
      DO 130 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(J,I,K).LE.0) GO TO 130
      NEQ=NEQ+1
      IEQPNT(J,I,K)=NEQ
  130 CONTINUE
      NUP=NEQ
C
C3B-----Order equations with even plane numbers.
      DO 140 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NLAY) K1=NLAY
      K2=N-NCOL-NROW
      IF(K2.LT.1) K2=1
      DO 140 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NROW) I1=NROW
      I2=N-K-NCOL
      IF(I2.LT.1) I2=1
      DO 140 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(J,I,K).LE.0) GO TO 140
      NEQ=NEQ+1
      IEQPNT(J,I,K)=NEQ
  140 CONTINUE
      NLOW=NEQ-NUP
      IF ( NLOW.LT.1 ) NLOW = 1  !JDH 20110814
      RETURN
C
C4------DIRECTION 2 NROW>NCOL>or=NLAY
C4A-----Order equations with odd plane numbers.
  200 DO 230 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NLAY) K1=NLAY
      K2=N-NCOL-NROW
      IF(K2.LT.1) K2=1
      DO 230 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NCOL) I1=NCOL
      I2=N-K-NROW
      IF(I2.LT.1) I2=1
      DO 230 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(I,J,K).LE.0) GO TO 230
      NEQ=NEQ+1
      IEQPNT(I,J,K)=NEQ
  230 CONTINUE
      NUP=NEQ
C
C4B-----Order equations with even plane numbers.
      DO 240 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NLAY) K1=NLAY
      K2=N-NCOL-NROW
      IF(K2.LT.1) K2=1
      DO 240 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NCOL) I1=NCOL
      I2=N-K-NROW
      IF(I2.LT.1) I2=1
      DO 240 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(I,J,K).LE.0) GO TO 240
      NEQ=NEQ+1
      IEQPNT(I,J,K)=NEQ
  240 CONTINUE
      NLOW=NEQ-NUP
      RETURN
C
C-------DIRECTION 3 NCOL>or=NLAY>NROW
C5A-----Order equations with odd plane numbers.
  300 DO 330 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NROW) K1=NROW
      K2=N-NCOL-NLAY
      IF(K2.LT.1) K2=1
      DO 330 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NLAY) I1=NLAY
      I2=N-K-NCOL
      IF(I2.LT.1) I2=1
      DO 330 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(J,K,I).LE.0) GO TO 330
      NEQ=NEQ+1
      IEQPNT(J,K,I)=NEQ
  330 CONTINUE
      NUP=NEQ
C
C5B-----Order equations with even plane numbers.
      DO 340 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NROW) K1=NROW
      K2=N-NCOL-NLAY
      IF(K2.LT.1) K2=1
      DO 340 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NLAY) I1=NLAY
      I2=N-K-NCOL
      IF(I2.LT.1) I2=1
      DO 340 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(J,K,I).LE.0) GO TO 340
      NEQ=NEQ+1
      IEQPNT(J,K,I)=NEQ
  340 CONTINUE
      NLOW=NEQ-NUP
      RETURN
C
C6------DIRECTION 4 NLAY>NCOL>or=NROW
C6A-----Order equations with odd plane numbers.
  400 DO 430 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NROW) K1=NROW
      K2=N-NCOL-NLAY
      IF(K2.LT.1) K2=1
      DO 430 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NCOL) I1=NCOL
      I2=N-K-NLAY
      IF(I2.LT.1) I2=1
      DO 430 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(I,K,J).LE.0) GO TO 430
      NEQ=NEQ+1
      IEQPNT(I,K,J)=NEQ
  430 CONTINUE
      NUP=NEQ
C
C6B-----Order equations with even plane numbers.
      DO 440 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NROW) K1=NROW
      K2=N-NCOL-NLAY
      IF(K2.LT.1) K2=1
      DO 440 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NCOL) I1=NCOL
      I2=N-K-NLAY
      IF(I2.LT.1) I2=1
      DO 440 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(I,K,J).LE.0) GO TO 440
      NEQ=NEQ+1
      IEQPNT(I,K,J)=NEQ
  440 CONTINUE
      NLOW=NEQ-NUP
      RETURN
C
C7------DIRECTION 5 NROW>or=NLAY>NCOL
C7A-----Order equations with odd plane numbers.
  500 DO 530 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NCOL) K1=NCOL
      K2=N-NROW-NLAY
      IF(K2.LT.1) K2=1
      DO 530 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NLAY) I1=NLAY
      I2=N-K-NROW
      IF(I2.LT.1) I2=1
      DO 530 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(K,J,I).LE.0) GO TO 530
      NEQ=NEQ+1
      IEQPNT(K,J,I)=NEQ
  530 CONTINUE
      NUP=NEQ
C
C7B-----Order equations with even plane numbers.
      DO 540 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NCOL) K1=NCOL
      K2=N-NROW-NLAY
      IF(K2.LT.1) K2=1
      DO 540 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NLAY) I1=NLAY
      I2=N-K-NROW
      IF(I2.LT.1) I2=1
      DO 540 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(K,J,I).LE.0) GO TO 540
      NEQ=NEQ+1
      IEQPNT(K,J,I)=NEQ
  540 CONTINUE
      NLOW=NEQ-NUP
      RETURN
C
C-------DIRECTION 6 NLAY>NROW>NCOL
C8A-----Order equations with odd plane numbers.
  600 DO 630 N=3,NPLANE,2
      K1=N-2
      IF(K1.GT.NCOL) K1=NCOL
      K2=N-NROW-NLAY
      IF(K2.LT.1) K2=1
      DO 630 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NROW) I1=NROW
      I2=N-K-NLAY
      IF(I2.LT.1) I2=1
      DO 630 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(K,I,J).LE.0) GO TO 630
      NEQ=NEQ+1
      IEQPNT(K,I,J)=NEQ
  630 CONTINUE
      NUP=NEQ
C
C8B-----Order equations with even plane numbers.
      DO 640 N=4,NPLANE,2
      K1=N-2
      IF(K1.GT.NCOL) K1=NCOL
      K2=N-NROW-NLAY
      IF(K2.LT.1) K2=1
      DO 640 K=K1,K2,-1
      I1=N-K-1
      IF(I1.GT.NROW) I1=NROW
      I2=N-K-NLAY
      IF(I2.LT.1) I2=1
      DO 640 I=I1,I2,-1
      J=N-K-I
      IF(IBOUND(K,I,J).LE.0) GO TO 640
      NEQ=NEQ+1
      IEQPNT(K,I,J)=NEQ
  640 CONTINUE
      NLOW=NEQ-NUP
      RETURN
C
      END
      SUBROUTINE SDE47S(AU,AL,IUPPNT,D4B,NUP,NLOW,NEQ,MXBW,NBW,ITYPE,
     1              ID4DIM)
C     ******************************************************************
C     SOLVE EQUATIONS USING GAUSS ELIMINATION ASSUMING D4 ORDERING
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION AU(ID4DIM,NUP),AL(MXBW,NLOW),IUPPNT(ID4DIM,NUP),D4B(NEQ)
C     ------------------------------------------------------------------
C
C1------DEFINE CONSTANTS.
      ONE=1.
      ZERO=0.
      NLOWM1=NLOW-1
C
C2------DON'T ELIMINATE UNLESS NECESSARY.
      IF (ITYPE.EQ.1) GO TO 80
C
C2A-----MUST ELIMINATE -- DO THIS IN TWO PARTS.
C2A-----ELIMINATE THE LEFT SIDE OF THE LOWER PART OF [A] TO FILL [AL].
      DO 40 I=1,NUP
      JJ=IUPPNT(1,I)
      C1=ONE/AU(1,I)
      DO 30 J=2,JJ
      LR=IUPPNT(J,I)
      L=LR-NUP
      C=AU(J,I)*C1
      DO 20 K=J,JJ
      KL=IUPPNT(K,I)-LR+1
      AL(KL,L)=AL(KL,L)-C*AU(K,I)
   20 CONTINUE
      AU(J,I)=C
   30 CONTINUE
   40 CONTINUE
C
C2B-----ELIMINATE [AL].
      DO 70 I=1,NLOWM1
      L=I
      C1=ONE/AL(1,I)
      DO 60 J=2,NBW
      L=L+1
      IF (AL(J,I).EQ.ZERO) GO TO 60
      C=AL(J,I)*C1
      KL=0
      DO 50 K=J,NBW
      KL=KL+1
      IF (AL(K,I).NE.ZERO) AL(KL,L)=AL(KL,L)-C*AL(K,I)
   50 CONTINUE
      AL(J,I)=C
   60 CONTINUE
   70 CONTINUE
C
C3------D4B MUST ALWAYS BE MODIFIED -- MODIFY D4B IN TWO PARTS.
C3A-----MODIFY D4B DUE TO ELIMINATION TO FILL [AL].
   80 DO 100 I=1,NUP
      JJ=IUPPNT(1,I)
      DO 90 J=2,JJ
      LR=IUPPNT(J,I)
      D4B(LR)=D4B(LR)-AU(J,I)*D4B(I)
   90 CONTINUE
      D4B(I)=D4B(I)/AU(1,I)
  100 CONTINUE
C
C3B-----MODIFY D4B DUE TO ELIMINATION OF [AL].
      DO 120 I=1,NLOWM1
      IR=I+NUP
      LR=IR
      DO 110 J=2,NBW
      LR=LR+1
      IF (AL(J,I).NE.ZERO) D4B(LR)=D4B(LR)-AL(J,I)*D4B(IR)
  110 CONTINUE
      D4B(IR)=D4B(IR)/AL(1,I)
  120 CONTINUE
C
C4------BACK SUBSTITUTE LOWER PART.
!      D4B(NEQ)=D4B(NEQ)/AL(1,NEQ-NUP)
      IJDH = NEQ-NUP                  !JDH 20110814
      IF ( IJDH.LT.1 ) THEN           !JDH 20110814
        IJDH = 1                      !JDH 20110814
      END IF                          !JDH 20110814
      D4B(NEQ)=D4B(NEQ)/AL(1,IJDH)    !JDH 20110814
      DO 140 I=1,NLOWM1
      K=NEQ-I
      KL=K-NUP
      L=K
      DO 130 J=2,NBW
      L=L+1
      IF (AL(J,KL).NE.ZERO) D4B(K)=D4B(K)-AL(J,KL)*D4B(L)
  130 CONTINUE
  140 CONTINUE
C
C5------BACK SUBSTITUTE UPPER PART.
      DO 160 I=1,NUP
      K=NUP+1-I
      JJ=IUPPNT(1,K)
      DO 150 J=2,JJ
      L=IUPPNT(J,K)
      D4B(K)=D4B(K)-AU(J,K)*D4B(L)
  150 CONTINUE
  160 CONTINUE
C
C6------RETURN.
      RETURN
      END
      SUBROUTINE SDE47P(HDCGDE4,LRCHDE4,NUMIT,IOUT,NCOL,NROW)
C     ******************************************************************
C     PRINT MAXIMUM HEAD CHANGE DURING EACH D4 ITERATION FOR A TIME STEP
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION HDCGDE4(NUMIT), LRCHDE4(3,NUMIT)
C     ------------------------------------------------------------------
C
      IF (NCOL.LE.999 .AND. NROW.LE.999) THEN
        WRITE(IOUT,5)
    5   FORMAT(1X,/1X,'MAXIMUM HEAD CHANGE FOR EACH ITERATION:',/
     1    1X,/1X,3('  HEAD CHANGE  LAY,ROW,COL'),/1X,79('-'))
        WRITE (IOUT,10) (HDCGDE4(J),(LRCHDE4(I,J),I=1,3),J=1,NUMIT)
   10   FORMAT((2X,3(2X,1PG10.3,' (',I3,',',I3,',',I3,')')))
      ELSE
        WRITE(IOUT,15)
   15   FORMAT(1X,/1X,'MAXIMUM HEAD CHANGE FOR EACH ITERATION:',/
     1    1X,/1X,2('    HEAD CHANGE  LAY,ROW,COL  '),/1X,60('-'))
        WRITE (IOUT,20) (HDCGDE4(J),(LRCHDE4(I,J),I=1,3),J=1,NUMIT)
   20   FORMAT((2X,2(2X,1PG10.3,' (',I3,',',I5,',',I5,')')))
      ENDIF
      WRITE(IOUT,30)
   30 FORMAT(1X,/1X)
C
      RETURN
      END
      SUBROUTINE DE47DA(IGRID)
C  Deallocate DE4 data
      USE DE4MODULE
C
      DEALLOCATE(DE4DAT(IGRID)%MXUP)
      DEALLOCATE(DE4DAT(IGRID)%MXLOW)
      DEALLOCATE(DE4DAT(IGRID)%MXEQ)
      DEALLOCATE(DE4DAT(IGRID)%MXBW)
      DEALLOCATE(DE4DAT(IGRID)%ITMX)
      DEALLOCATE(DE4DAT(IGRID)%ID4DIR)
      DEALLOCATE(DE4DAT(IGRID)%NITERDE4)
      DEALLOCATE(DE4DAT(IGRID)%IFREQ)
      DEALLOCATE(DE4DAT(IGRID)%IPRD4)
      DEALLOCATE(DE4DAT(IGRID)%MUTD4)
      DEALLOCATE(DE4DAT(IGRID)%ID4DIM)
      DEALLOCATE(DE4DAT(IGRID)%NBWL)
      DEALLOCATE(DE4DAT(IGRID)%NUPL)
      DEALLOCATE(DE4DAT(IGRID)%NLOWL)
      DEALLOCATE(DE4DAT(IGRID)%NLOW)
      DEALLOCATE(DE4DAT(IGRID)%NEQ)
      DEALLOCATE(DE4DAT(IGRID)%NUP)
      DEALLOCATE(DE4DAT(IGRID)%NBW)
      DEALLOCATE(DE4DAT(IGRID)%DELTL)
      DEALLOCATE(DE4DAT(IGRID)%ACCLDE4)
      DEALLOCATE(DE4DAT(IGRID)%HCLOSEDE4)
      DEALLOCATE(DE4DAT(IGRID)%IUPPNT)
      DEALLOCATE(DE4DAT(IGRID)%IEQPNT)
      DEALLOCATE(DE4DAT(IGRID)%AU)
      DEALLOCATE(DE4DAT(IGRID)%AL)
      DEALLOCATE(DE4DAT(IGRID)%D4B)
      DEALLOCATE(DE4DAT(IGRID)%HDCGDE4)
      DEALLOCATE(DE4DAT(IGRID)%LRCHDE4)
C
      RETURN
      END
      SUBROUTINE DE47PNT(IGRID)
C  Set pointers to DE4 data for grid.
      USE DE4MODULE
C
      MXUP=>DE4DAT(IGRID)%MXUP
      MXLOW=>DE4DAT(IGRID)%MXLOW
      MXEQ=>DE4DAT(IGRID)%MXEQ
      MXBW=>DE4DAT(IGRID)%MXBW
      ITMX=>DE4DAT(IGRID)%ITMX
      ID4DIR=>DE4DAT(IGRID)%ID4DIR
      NITERDE4=>DE4DAT(IGRID)%NITERDE4
      IFREQ=>DE4DAT(IGRID)%IFREQ
      IPRD4=>DE4DAT(IGRID)%IPRD4
      MUTD4=>DE4DAT(IGRID)%MUTD4
      ID4DIM=>DE4DAT(IGRID)%ID4DIM
      NBWL=>DE4DAT(IGRID)%NBWL
      NUPL=>DE4DAT(IGRID)%NUPL
      NLOWL=>DE4DAT(IGRID)%NLOWL
      NLOW=>DE4DAT(IGRID)%NLOW
      NEQ=>DE4DAT(IGRID)%NEQ
      NUP=>DE4DAT(IGRID)%NUP
      NBW=>DE4DAT(IGRID)%NBW
      DELTL=>DE4DAT(IGRID)%DELTL
      ACCLDE4=>DE4DAT(IGRID)%ACCLDE4
      HCLOSEDE4=>DE4DAT(IGRID)%HCLOSEDE4
      IUPPNT=>DE4DAT(IGRID)%IUPPNT
      IEQPNT=>DE4DAT(IGRID)%IEQPNT
      AU=>DE4DAT(IGRID)%AU
      AL=>DE4DAT(IGRID)%AL
      D4B=>DE4DAT(IGRID)%D4B
      HDCGDE4=>DE4DAT(IGRID)%HDCGDE4
      LRCHDE4=>DE4DAT(IGRID)%LRCHDE4
C
      RETURN
      END
      SUBROUTINE DE47PSV(IGRID)
C  Save pointers to DE4 data.
      USE DE4MODULE
C
      DE4DAT(IGRID)%MXUP=>MXUP
      DE4DAT(IGRID)%MXLOW=>MXLOW
      DE4DAT(IGRID)%MXEQ=>MXEQ
      DE4DAT(IGRID)%MXBW=>MXBW
      DE4DAT(IGRID)%ITMX=>ITMX
      DE4DAT(IGRID)%ID4DIR=>ID4DIR
      DE4DAT(IGRID)%NITERDE4=>NITERDE4
      DE4DAT(IGRID)%IFREQ=>IFREQ
      DE4DAT(IGRID)%IPRD4=>IPRD4
      DE4DAT(IGRID)%MUTD4=>MUTD4
      DE4DAT(IGRID)%ID4DIM=>ID4DIM
      DE4DAT(IGRID)%NBWL=>NBWL
      DE4DAT(IGRID)%NUPL=>NUPL
      DE4DAT(IGRID)%NLOWL=>NLOWL
      DE4DAT(IGRID)%NLOW=>NLOW
      DE4DAT(IGRID)%NEQ=>NEQ
      DE4DAT(IGRID)%NUP=>NUP
      DE4DAT(IGRID)%NBW=>NBW
      DE4DAT(IGRID)%DELTL=>DELTL
      DE4DAT(IGRID)%ACCLDE4=>ACCLDE4
      DE4DAT(IGRID)%HCLOSEDE4=>HCLOSEDE4
      DE4DAT(IGRID)%IUPPNT=>IUPPNT
      DE4DAT(IGRID)%IEQPNT=>IEQPNT
      DE4DAT(IGRID)%AU=>AU
      DE4DAT(IGRID)%AL=>AL
      DE4DAT(IGRID)%D4B=>D4B
      DE4DAT(IGRID)%HDCGDE4=>HDCGDE4
      DE4DAT(IGRID)%LRCHDE4=>LRCHDE4
C
      RETURN
      END

max_stars_repo_name: cneyens/MODFLOW-MAXDRN

max_issues_repo_name: cneyens/MODFLOW-MAXDRN

max_forks_repo_name: cneyens/MODFLOW-MAXDRN


text:
      Subroutine rddata_line_elem(ncid,imax,jmax,nch
     &,sounding,istatus)
c
c routine designed to read 3-d satellite sounding netCDF data with variable
c line and element dimensions. Should work for 2d data with variable dimensions
c if nch = 1
c
      Implicit None
      include 'netcdf.inc'
      Integer   RCODE
c
      Integer   imax,jmax,nch
      Integer   sounding(imax, jmax, nch)
      Integer   sounding_rec(imax)

      Integer i,j,k
      Integer varid,ncid
      Integer istatus
      Integer dim_id_x
      Integer dim_id_y
      Integer dim_id_k
      Integer i4dum

      Integer START(10)
      Integer COUNT(10)
      Character*31 DUMMY
c
C **************************************************************************
c
      istatus = 1
c
c Now ready to read the data
c
      START(1)=1
      COUNT(2)=1
      COUNT(3)=1
c
c get sounding variable id
c
      rcode=NF_INQ_VARID(ncid,'sounding',varid)
      if(rcode.ne.0)then
         write(6,*)'Error getting sounding varid'
         istatus = -1
         return
      endif

      do k = 1,nch  !# of channels in sounding database (= 19).

         write(6,*)'Reading channel ',k

         START(3)=k

         do j = 1,jmax

            COUNT(1)=imax
            START(2)=j

c read line
      rcode=NF_GET_VARA_INT(NCID,varid,START,COUNT,sounding_rec)
            if(rcode.ne.0)then
               print*,'Error reading sounding database ',rcode
               istatus = -1
               print*,'Returning to main, istatus = ',istatus
               return
            endif
c
c load line into the 3d output array
c
            sounding(:,j,k) = sounding_rec(:)

         enddo
      enddo
c -------------------------------------
C
      write(6,*)'Sucessful reading sounding - rddata_line_elem '

      Return
      End

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense