text:
c **********************************************************************
      SUBROUTINE RESP00(kresis,ksave,kgem,kstep,ndof,kst,kenr,ener,
     1     ened,enso,beto,relas,rdamp,rinit,ddise,
     2     dise,vele)
c *********************************************************************
c  DRAIN-2DX ROTATIONAL SPRING ELEMENT - Bilinear return map
c ---------------------------------------------------------------------
c  PURPOSE
c     The purpose of the subroutine resp00() is to update the element
c     state, form the static and damping resisting forces, perform
c     energy calculations, update response quantities, and put
c     element results in /thelm/ for saving or printing.
c ---------------------------------------------------------------------
c  CALLED FROM:    respxx
c ---------------------------------------------------------------------
c  ARGUMENTS
c   Input:
c     kresis      - indicator for calculating resisting forces
c                   (1=static, 2=static and damping)
c     ksave       - indicator for saving element results
c                   (1=yes, 0=no)
c     kgem        - second-order analysis code
c                   (not used for this element)
c     ndof        - number of element DOFs
c     kenr        - energy calculation indicator
c                   (0=none, 1=static, 2=static and dynamic)
c     beto        - initial stiffness damping factor
c     ddise(ndof) - element nodal incremental displacement vector
c     dise(ndof)  - element total nodal displacement vector
c     vele(ndof)  - element nodal velocity vector
c   Output:
c     ener        - change in element elasto-plastic energy
c     ened        - change in element damping energy
c                   (not computed)
c     enso        - change in element second-order energy
c                   (always zero for this element)
c     relas(ndof) - element static resisting force vector
c     rdamp(ndof) - element damping resisting force vector
c                   (not computed)
c     rinit(ndof) - element initial resisting force vector
c                   (not computed)
c   Modify:
c     kstep       - step number in this segment
c     kst         - stiffness formation code
c                   (1=yes, 0=no)
c ---------------------------------------------------------------------
c  DOUBLE PRECISION
      implicit none
c ---------------------------------------------------------------------
c  ARGUMENT TYPES
      integer   kresis,ksave,kgem,kstep,ndof,kst,kenr
      real*8    ener,ened,enso,beto
      real*8    relas(ndof),rdamp(ndof),rinit(ndof)
      real*8    ddise(ndof),dise(ndof),vele(ndof)

c ---------------------------------------------------------------------
c     LABELLED COMMONS
      real*8 E,sigY,Hiso,Hkin
      real*8 ep,alpha,kappa
      real*8 epsP,sigP,tangP
      real*8 tang
      include 'infel00.h'

      integer sgn
      real*8 eps,sig,xsi,f,dGamma

c     Current strain
      eps = dise(2)-dise(1)

c     Elastic predictor
      sig = E * (eps - ep)

c     Stress relative to back stress
      xsi = sig - kappa

c     Yield function
      f = dabs(xsi) - (sigY + Hiso*alpha)

c     Inside yield surface
c					 if (f <= 0.0) then
      if (f.le.0.0) then
         tang = E
c     Outside yield surface ... do return mapping
      else
c     Consistency parameter
         dGamma = f / (E+Hiso+Hkin)

c     Normal to yield surface
c         if (xsi <= 0.d0) then
         if (xsi.le.0.d0) then
            sgn = -1
         else
            sgn = 1
         endif

c     Updated stress
         sig = sig - dGamma*E*sgn

c     Updated plastic strain
         ep = ep + dGamma*sgn

c     Updated back stress
         kappa = kappa + dGamma*Hkin*sgn

c     Updated internal hardening variable
         alpha = alpha + dGamma

c     Elasto-plastic tangent
         tang = E*(Hkin+Hiso) / (E+Hkin+Hiso)

      endif

c     Populate force vectors
      relas(2) =  sig
      relas(1) = -sig

c     Not computed
      rinit(2) = 0.d0
      rinit(1) = 0.d0

c     Not computed
      rdamp(2) = 0.d0
      rdamp(1) = 0.d0

      ener = 0.5d0*(sig+sigP)*(ddise(2)-ddise(1))

c     Not computed
      ened = 0.d0
      enso = 0.d0

      epsP = eps
      sigP = sig

      END

c **********************************************************************
      SUBROUTINE STIF00 (kstt,ktype,ndof,fk)
c *********************************************************************
c  DRAIN-2DX ROTATIONAL SPRING ELEMENT - Bilinear return map
c ---------------------------------------------------------------------
c  PURPOSE
c     The purpose of the subroutine stif00() is to form the stiffness
c     matrix and return it to the variable fk.
c ---------------------------------------------------------------------
c  CALLED FROM:  stifxx
c ---------------------------------------------------------------------
c  ARGUMENTS
c    Input:
c      kstt           - (1=total, 0=change in) stiffness
c      ktype          - 1 = elastic stiffness only,
c                       0 = elastic + geometric stiffness
c                      -1 = geometric stiffness only
c      ndof           - number of element DOFs
c    Output:
c      fk(ndof, ndof) - stiffness matrix
c ---------------------------------------------------------------------
c  DOUBLE PRECISION
      implicit none
c ---------------------------------------------------------------------
c  ARGUMENT TYPES
      integer   kstt,ktype,ndof
      real*8    fk(ndof,ndof)

c ---------------------------------------------------------------------
c  LABELLED COMMONS
      real*8 E,sigY,Hiso,Hkin
      real*8 ep,alpha,kappa
      real*8 epsP,sigP,tangP
      real*8 tang
      include 'infel00.h'

c
c TOTAL STIFFNESS MATRIX
c
      fk(1,1) =  tang
      fk(1,2) = -tang
      fk(2,1) = -tang
      fk(2,2) =  tang

c
c CHANGE IN STIFFNESS MATRIX
c
      if(kstt.EQ.0) then
         fk(1,1) =  tang-tangP
         fk(1,2) = -fk(1,1)
         fk(2,1) =  fk(1,2)
         fk(2,2) =  fk(1,1)
      endif

      tang = tangP
c
c NO GEOMETRIC STIFFNESS
c
      if(ktype.EQ.-1) then
         fk(1,1) = 0.d0
         fk(1,2) = 0.d0
         fk(2,1) = 0.d0
         fk(2,2) = 0.d0
      end if

      RETURN
      END

ext: f

max_stars_repo_name: kuanshi/ductile-fracture


text:
*DECK DLGAMS
      SUBROUTINE DLGAMS (X, DLGAM, SGNGAM)
C***BEGIN PROLOGUE  DLGAMS
C***PURPOSE  Compute the logarithm of the absolute value of the Gamma
C            function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      DOUBLE PRECISION (ALGAMS-S, DLGAMS-D)
C***KEYWORDS  ABSOLUTE VALUE OF THE LOGARITHM OF THE GAMMA FUNCTION,
C             FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DLGAMS(X,DLGAM,SGNGAM) calculates the double precision natural
C logarithm of the absolute value of the Gamma function for
C double precision argument X and stores the result in double
C precision argument DLGAM.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DLNGAM
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***END PROLOGUE  DLGAMS
      DOUBLE PRECISION X, DLGAM, SGNGAM, DLNGAM
C***FIRST EXECUTABLE STATEMENT  DLGAMS
      DLGAM = DLNGAM(X)
      SGNGAM = 1.0D0
      IF (X.GT.0.D0) RETURN
C
      INT = MOD (-AINT(X), 2.0D0) + 0.1D0
      IF (INT.EQ.0) SGNGAM = -1.0D0
C
      RETURN
      END

max_issues_repo_name: GEHC-Surgery/ITK

max_forks_repo_name: GEHC-Surgery/ITK


text:
      SUBROUTINE SDRHT
C
C     SPECIAL FLUX-DATA-RECOVERY MODULE FOR HBDY ELEMENTS IN HEAT
C     TRANSFER.
C
C     DMAP CALLING SEQUENCE.
C
C     SDRHT SIL,USET,UGV,OEF1,SLT,EST,DIT,QGE,DLT,/OEF1X/V,N,TABS $
C
      LOGICAL         HAVIDS,CARDIN,LHBDY,TRANST,FOUND,MCH521
      INTEGER         TABLST(13),Z,BUF(50),SYSBUF,RD,RDREW,WRT,WRTREW,
     1                CLSREW,CLS,SUBR(2),NAME(2),EOR,IDPOS(3),OUTPT,
     2                SLTYPS,LDWORD(16),UG,OEF1,SLT,EST,DIT,QGE,DLT,
     3                OEF1X,BUF1,BUF2,BUF3,CORE,PASS,HBDYTP,MCBUGV(7),
     4                FILE,ELTYPE,ESTWDS,ECPT(100),SLTAT,SLTREC,MCB(7),
     5                EOL,GSIZE
      REAL            RZ(1),RBUF(50),GRIDS(6)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25,SWM*27
      COMMON /XMSSG / UFM,UWM,UIM,SFM,SWM
      COMMON /CONDAS/ CONSTS(5)
      COMMON /SYSTEM/ KSYSTM(65)
      COMMON /NAMES / RD,RDREW,WRT,WRTREW,CLSREW,CLS
      COMMON /ZNTPKX/ AI(4),IROW,EOL
      COMMON /ZZZZZZ/ Z(1)
      COMMON /MACHIN/ MACH
      COMMON /BLANK / TABS
      EQUIVALENCE     (KSYSTM(1),SYSBUF),(KSYSTM(2),OUTPT),
     1                (CONSTS(2),TWOPI ),(CONSTS(3),RADDEG),
     2                (Z(1),RZ(1)), (BUF(1),RBUF(1))
      DATA    TABLST / 4, 1105,11,1, 1205,12,2, 1305,13,3, 1405,14,4  /
      DATA    LENTRY / 14   /, EOR,NOEOR  / 1, 0 /, SUBR/ 4HSDRH,2HT  /
      DATA    IDPOS  / 2,1,5/, HBDYTP/ 52 /,     SLTYPS / 16 /
      DATA    LDWORD / 6,6,4,4,6,6,2,5,5,6,6,7,2,2,5,5  /
      DATA    UG,OEF1, SLT,EST,DIT,QGE,DLT,OEF1X /
     1        103,104, 105,106,107,108,109, 201  /
      DATA    GRIDS  / 1.0, 2.0, 2.0, 3.0, 4.0, 2.0     /
C
C     SET UP CORE AND BUFFERS
C
      BUF1   = KORSZ(Z) - SYSBUF - 2
      BUF2   = BUF1 - SYSBUF - 2
      BUF3   = BUF2 - SYSBUF - 2
      CORE   = BUF3 - 1
      IDREC  = 1
      IELTAB = 1
      PASS   = 0
      HAVIDS = .FALSE.
      CARDIN = .FALSE.
      MCH521 = MACH.EQ.5 .OR. MACH.EQ.21
C
C     OPEN INPUT FORCES -OEF1- AND OUTPUT FORCES -OEF1X-.
C
      CALL OPEN  (*1190,OEF1,Z(BUF1),RDREW)
      CALL OPEN  (*1200,OEF1X,Z(BUF2),WRTREW)
      CALL FNAME (OEF1X,NAME)
      CALL WRITE (OEF1X,NAME,2,EOR)
      CALL FWDREC (*1220,OEF1)
C
C     COPY RECORD PAIRS OF DATA FROM OEF1 TO OEF1X UNTIL HBDY DATA IS
C     DETECTED.
C
   10 LCORE = CORE - IDREC
      IF (LCORE .LT. 300) CALL MESAGE (-8,0,SUBR)
      FILE = OEF1
      CALL READ (*1180,*30,OEF1,Z(IDREC),LCORE,NOEOR,IFLAG)
      WRITE  (OUTPT,20) SWM
   20 FORMAT (A27,' 3063, INPUT FORCES DATA BLOCK HAS INCORRECT DATA.')
      GO TO 1220
C
C     MODIFY ID-RECORD IF THIS IS FOR HBDY ELEMENTS.
C
   30 IF (Z(IDREC+2) .NE. HBDYTP) GO TO 40
      LHBDY = .TRUE.
      GO TO 50
   40 LHBDY = .FALSE.
      GO TO 90
C
C     SET CONSTANTS FROM OEF1 ID RECORD.
C
   50 IF (Z(IDREC)/10 .EQ. 6) GO TO 70
      LOADID = Z(IDREC+7)
      TRANST = .FALSE.
      GO TO 80
   70 LOADID = Z(IDREC+7)
      TIME   = RZ(IDREC+4)
      TRANST = .TRUE.
   80 LWORDS = Z(IDREC+9)
      Z(IDREC+9) = 5
C/////
C     CALL BUG (4HTRAN ,90,TRANST,1)
C/////
   90 CALL WRITE (OEF1X,Z(IDREC),IFLAG,EOR)
      IF (LHBDY) GO TO 120
C
C     NOT AN HBDY ELEMENT TYPE THUS COPY DATA ACROSS.
C
  100 CALL READ  (*1260,*110,OEF1,Z(IDREC),LCORE,NOEOR,IFLAG)
      CALL WRITE (OEF1X,Z(IDREC),LCORE,NOEOR)
      GO TO 100
  110 CALL WRITE (OEF1X,Z(IDREC),IFLAG,EOR)
      GO TO 10
C
C     HBDY ELEMENT DATA ENCOUNTERED.
C
  120 PASS = PASS + 1
C
C     ON FIRST PASS ELEMENT-DATA-TABLE IS FORMED.
C
C     EACH ENTRY WILL CONTAIN,
C
C      1) ELEMENT-ID (HBDY).
C      2) FLUX-RADIATION TERM FOR THIS ELEMENT.
C      3) FLUX-X FROM OEF1 DATA.
C      4) APPLIED LOAD (USING SLT DATA).
C      5) HBDY ELEMENT TYPE (1 TO 6).
C      6) HBDY AREA FACTOR.
C      7) ALPHA VALUE.
C      8) V1(1) *
C      9) V1(2)  *  VECTOR-V1
C     10) V1(3) *
C     11) V2(1) *
C     12) V2(2)  *  VECTOR-V2
C     13) V2(3) *
C     14) OUTPUT ID*10 + DEVICE CODE (FROM OEF1).
C
C     ON PASSES OTHER THAN THE FIRST ONLY THE FLUX-X VALUE IS EXTRACTED
C     FROM OEF1-HBDY-ENTRIES.
C
      JELTAB = IELTAB - 1
C
C     INPUT = ID*10+CODE, NAME1,NAME2,GRD-X,GRD-Y,GRD-Z,FLUX-X,FLUX-Y,
C             FLUX-Z   (TOTAL OF 9 WORDS)
C
  130 CALL READ (*1260,*150,OEF1,BUF,LWORDS,NOEOR,IFLAG)
      IF (PASS .GT. 1) GO TO 140
      Z(JELTAB+1) = BUF(1)/10
C
C     STORE (FLUX-X) AND (OUTPUT ID*10 + DEVICE CODE).
C
  140 Z(JELTAB+ 3) = BUF(7)
      RZ(JELTAB+2) = 0.0
      Z(JELTAB+14) = BUF(1)
      JELTAB = JELTAB + LENTRY
      IF (JELTAB .GT. CORE) CALL MESAGE (-8,0,SUBR)
      GO TO 130
C
C     END OF DATA.
C
  150 IF (PASS .GT. 1) GO TO 160
      NELTAB = JELTAB
      NUMBER = (NELTAB - IELTAB + 1)/LENTRY
C
C     OPEN UG FILE FOR INPUT OF UG VECTORS.
C
      CALL OPEN (*1220,UG,Z(BUF3),RDREW)
      CALL FWDREC (*1220,UG)
      MCBUGV(1) = UG
      CALL RDTRL (MCBUGV)
      GSIZE = MCBUGV(3)
      GO TO 180
  160 IF (JELTAB .EQ. NELTAB) GO TO 180
      WRITE  (OUTPT,170) SWM,JELTAB,NELTAB
  170 FORMAT (A27,' 3064, INCONSISTANT HBDY DATA RECORDS.  ',2I20)
      GO TO 1220
C
C     ALL DATA FROM OEF1 IS AT HAND NOW.
C     FILES ARE CLOSED WITHOUT REWIND.
C
  180 CALL CLOSE (OEF1,CLS)
      CALL CLOSE (OEF1X,CLS)
C
C     ALLOCATE UG VECTOR SPACE ON FIRST PASS.
C
      IF (PASS .NE. 1) GO TO 190
      IUG  = NELTAB + 1
      NUG  = NELTAB + GSIZE
      IUGZ = IUG - 1
      IF (NUG + 5 .GT. CORE) CALL MESAGE (-8,0,SUBR)
C
C     BRING NEXT DISPLACEMENT VECTOR INTO CORE.
C
  190 DO 200 I = IUG,NUG
      RZ(I) = TABS
  200 CONTINUE
C
      CALL INTPK (*220,UG,0,1,0)
  210 CALL ZNTPKI
      KK = IUGZ + IROW
      RZ(KK) = RZ(KK) + AI(1)
      IF (EOL) 210,210,220
C
C     RAISE VECTOR RESULT TO 4TH POWER
C
  220 DO 230 I = IUG,NUG
      RZ(I) = RZ(I)**4
  230 CONTINUE
C
C     IF TRANSIENT PROBLEM SKIP ACCELERATION AND VELOCITY VECTORS
C
      IF (.NOT. TRANST) GO TO 250
      DO 240 I = 1,2
      CALL FWDREC (*250,UG)
  240 CONTINUE
  250 CONTINUE
C/////
C     CALL BUG (4HUG4  ,200,Z(IUG),NUG-IUG+1)
C/////
C
C     IF NONLINEAR PROBLEM, COMPUTE FLUX RADIATION TERMS.
C
C                                            T           4
C     (FLUX-RADIATION         ) = (Q           )(U +TABS)
C                    EL-SUBSET      G,EL-SUBSET   G
C
      FILE = QGE
      CALL OPEN (*261,QGE,Z(BUF1),RDREW)
      IF (PASS .EQ. 1) GO TO 260
      CALL FWDREC (*1260,QGE)
      GO TO 280
  261 IQGID = NUG + 1
      NQGID = NUG
      IDREC = NQGID
      GO TO 311
  260 CALL READ (*1260,*1270,QGE,BUF,-2,NOEOR,IFLAG)
C
C     ON FIRST PASS PICK UP ELEMENT ID LIST.
C
      IQGID = NUG + 1
      CALL READ (*1260,*270,QGE,Z(IQGID),CORE-IQGID,NOEOR,IFLAG)
      CALL MESAGE (-7,0,SUBR)
  270 NQGID = NUG + IFLAG
      IDREC = NQGID + 1
C/////
C     CALL BUG (4HQGID,410,Z(IQGID),NQGID-IQGID+1)
C/////
C
C     EACH FLUX-RADIATION TERM IN THE ELEMENT TABLE IS CREATED BY
C     FORMING THE DOT-PRODUCT OF THE COLUMN OF -QGE- HAVING THE
C     SAME ELEMENT-ID WITH THE -UG- VECTOR IN CORE.
C
  280 DO 310 I = IQGID,NQGID,1
      CALL INTPK (*310,QGE,0,1,0)
C
C     FIND OUT IF ID OF THIS VECTOR IS IN ELEMENT TABLE.
C
      KID = Z(I)
      CALL BISLOC (*300,KID,Z(IELTAB),LENTRY,NUMBER,JPOINT)
      JWORD = IELTAB + JPOINT
      Z(JWORD) = 0
C
C     FORM DOT PRODUCT
C
  290 CALL ZNTPKI
      K = IUGZ + IROW
      RZ(JWORD) = RZ(JWORD) - AI(1)*RZ(K)
      IF (EOL) 290,290,310
C
C     ID OF THIS COLUMN NOT IN ELEMENT TABLE
C
  300 CALL ZNTPKI
      IF (EOL) 300,300,310
  310 CONTINUE
C/////
C     CALL BUG (4HELTB ,440,Z(IELTAB),NELTAB-IELTAB+1)
C/////
      CALL CLOSE (QGE,CLSREW)
C
C     ON FIRST PASS, EST IS PASSED AND HBDY ELEMENTS CALLED.
C
  311 CONTINUE
      IF (PASS .GT. 1) GO TO 380
      FILE = EST
      CALL GOPEN (EST,Z(BUF1),RDREW)
      NEXT = IELTAB
C
C     READ THE ELEMENT TYPE
C
  320 CALL READ (*350,*1270,EST,ELTYPE,1,NOEOR,IFLAG)
      IF (ELTYPE .EQ. HBDYTP) GO TO 330
      CALL FWDREC (*1260,EST)
      GO TO 320
C
C     HBDY ELEMENT-SUMMARY-TABLE DATA FOUND.
C
  330 ESTWDS = 53
  340 CALL READ (*1260,*1270,EST,ECPT,ESTWDS,NOEOR,IFLAG)
C
C     CHECK TO SEE IF THIS ELEMENT IS IN OUTPUT SET.
C
      IF (Z(NEXT) - ECPT(1)) 350,370,340
  350 WRITE  (OUTPT,360) SWM,Z(NEXT)
  360 FORMAT (A27,' 3065, THERE IS NO EST DATA FOR HBDY ELEMENT ID =',
     1        I10)
      GO TO 1220
C
C     THIS ELEMENT IS IN TABLE.
C
  370 CALL HBDY (ECPT,ECPT,2,RBUF,BUF)
C
C     PLANT HBDY OUTPUTS INTO TABLE.
C
      Z(NEXT+ 4) = ECPT(2)
      Z(NEXT+ 5) = BUF(2)
      Z(NEXT+ 6) = ECPT(17)
      Z(NEXT+ 7) = BUF(11)
      Z(NEXT+ 8) = BUF(12)
      Z(NEXT+ 9) = BUF(13)
      Z(NEXT+10) = BUF(14)
      Z(NEXT+11) = BUF(15)
      Z(NEXT+12) = BUF(16)
      NEXT = NEXT + LENTRY
      IF (NEXT .LT. NELTAB) GO TO 340
      CALL CLOSE (EST,CLSREW)
C
C     LOAD SET PROCESSING IF LOAD-SET-ID IS NON-ZERO.
C
 380  IF (LOADID) 944,944,390
C
C     OPEN SLT FOR LOAD DATA.
C
C
  390 FILE = SLT
      CALL OPEN (*1250,SLT,Z(BUF1),RDREW)
      IF (HAVIDS) CALL FWDREC (*1260,SLT)
      IF (HAVIDS) GO TO 810
      HAVIDS = .TRUE.
      ILDID  = NQGID + 1
      NLDID  = ILDID - 1
C
C     IDS OF LOAD SETS NOT IN CORE THUS BRING IN IDS FROM HEADER RECORD.
C
      IMAST  = ILDID
      NSETS  = 0
      NLDSET = 3
      CALL READ (*1260,*1270,SLT,BUF,-2,NOEOR,IFLAG)
  400 IF (NLDID+5 .GT. CORE) CALL MESAGE (-8,0,SUBR)
      CALL READ (*1260,*410,SLT,Z(NLDID+1),1,NOEOR,IFLAG)
      NSETS = NSETS + 1
      Z(NLDID +2) = 1
      Z(NLDID +3) = Z(NLDID+1)
      RZ(NLDID+4) = 1.0
      Z(NLDID +5) = NSETS
      NLDID = NLDID + 5
      GO TO 400
C
C     IF TRANSIENT PROBLEM THEN DLT OPERATIONS BEGIN
C
  410 IF (.NOT.TRANST) GO TO 800
      FILE = DLT
      CALL OPEN (*1250,DLT,Z(BUF2),RDREW)
      CALL READ (*1260,*1270,DLT,BUF,3,NOEOR,IFLAG)
      M = BUF(3)
      FOUND = .FALSE.
      IF (M .LE. 0) GO TO 430
      DO 420 I = 1,M
      CALL READ (*1260,*1270,DLT,BUF,1,NOEOR,IFLAG)
      IF (BUF(1) .EQ. LOADID) FOUND = .TRUE.
  420 CONTINUE
C
C     NOW READ RLOAD1, RLOAD2, TLOAD1, AND TLOAD2 IDS.
C
  430 IRTIDS = NLDID + 1
      CALL READ (*1260,*440,DLT,Z(IRTIDS),CORE-IRTIDS,NOEOR,IFLAG)
      CALL MESAGE (-7,0,SUBR)
      GO TO 1220
  440 NRTIDS = NLDID + IFLAG
C
C     IF LOADID WAS FOUND AMONG THE DLOAD IDS, SEARCH IS NOW MADE IN
C     RECORD 1 OF THE DLT FOR THAT ID, AND ITS SUB-IDS.
C
      JJ1   = ILDID
      JJ2   = NLDID
      ILDID = NRTIDS + 1
      NLDID = NRTIDS + 2
      Z(ILDID  ) = LOADID
      Z(ILDID+1) = 0
      IF (.NOT. FOUND) GO TO 520
C
C     READ A MASTER DLOAD SET-ID.
C
  450 CALL READ (*1260,*1270,DLT,BUF,2,NOEOR,IFLAG)
      IF (BUF(1) .EQ. LOADID) GO TO 470
C
C     SKIP SUB-ID DATA OF THIS MASTER
C
  460 CALL READ (*1260,*1270,DLT,BUF,2,NOEOR,IFLAG)
      IF (BUF(2)) 450,460,460
C
C     MASTER-ID FOUND.  BUILD LOAD-SET-ID TABLE.
C
  470 FACTOR = RBUF(2)
  472 IF (NLDID+11 .LE. CORE) GO TO 480
      CALL MESAGE (8,0,SUBR)
      GO TO 1220
  480 CALL READ (*1260,*1270,DLT,BUF,2,NOEOR,IFLAG)
      IF (BUF(2)) 540,540,490
  490 Z(ILDID +1) = Z(ILDID+1) + 1
      Z(NLDID +1) = BUF(2)
      RZ(NLDID+2) = 0.0
      Z(NLDID +3) = 0
      RZ(NLDID+4) = RBUF(1)*FACTOR
      Z(NLDID +5) = 0
      NLDID = NLDID + 11
      GO TO 472
C
C     LOADID NOT AMONG DLOADS FOR THIS TRANSIENT PROBLEM
C
  520 Z(ILDID+1) = 1
      IF (NLDID+13 .LE. CORE) GO TO 530
      CALL MESAGE (8,0,SUBR)
      GO TO 1220
  530 Z(NLDID +1) = Z(ILDID)
      RZ(NLDID+2) = 0.0
      Z(NLDID +3) = 0
      RZ(NLDID+4) = 1.0
      Z(NLDID +5) = 0
      NLDID = NLDID + 11
C
C     IF THERE ARE ANY DLOAD CARDS AT ALL THEN BALANCE OF (OR ALL OF)
C     RECORD 1 IS NOW SKIPPED.
C
  540 IF (M .GT. 0) CALL FWDREC (*1260,DLT)
C
C     NOW PICKING UP DATA NEEDED OF DYNAMIC LOAD SET RECORDS.
C
      K1 = ILDID + 2
      K2 = NLDID
      DO 580 I = IRTIDS,NRTIDS
C
C     READ THE LOAD TYPE
C
      CALL READ (*1260,*1270,DLT,BUF,2,NOEOR,IFLAG)
      IF (BUF(1).NE.3 .AND. BUF(1).NE.4) GO TO 570
C
C     CHECK AND SEE IF THIS TLOAD ID IS AMONG THE SUB-IDS
C
      DO 550 J = K1,K2,11
      IF (Z(J) .EQ. Z(I)) GO TO 560
  550 CONTINUE
      GO TO 570
C
C     YES THIS RECORD IS NEEDED.  THUS PUT ITS DATA IN TABLE.
C
  560 Z(J+4) = BUF(1)
C
C     SLT ID INTO TABLE
C
      Z(J) = -BUF(2)
C
C     SET SLT RECORD NUMBER
C
      K = 0
      DO 565 L = JJ1,JJ2,5
      K = K + 1
      IF (Z(L) .EQ. BUF(2)) GO TO 566
  565 CONTINUE
      K = 0
  566 Z(J+2) = K
      CALL READ (*1260,*1270,DLT,Z(J+5),6,EOR,IFLAG)
      IF (BUF(1) .EQ. 3) Z(J+6) = 0
      GO TO 580
  570 CALL FWDREC (*1260,DLT)
  580 CONTINUE
C
C     CHECK IS NOW MADE TO INSURE ALL SUB-IDS RECEIVED DLT DATA.
C
C
C     SET SLT IDS POSITIVE
C
      DO 581 I = K1,K2,11
      Z(I) = IABS(Z(I))
  581 CONTINUE
      DO 610 I = K1,K2,11
      IF (Z(I+4)) 590,590,610
C
C     ERROR
C
  590 WRITE  (OUTPT,600) UWM,Z(I)
  600 FORMAT (A25,' 3066, THERE IS NO TLOAD1 OR TLOAD2 DATA FOR LOAD-',
     1       'ID =',I9)
  610 CONTINUE
C
      CALL CLOSE (DLT,CLSREW)
      NLDSET = 11
C
C     SORT SUB-ID TABLE ON SLT RECORD NUMBERS.
C
      CALL SORT (0,0,11,3,Z(K1),K2-K1+1)
C/////
C     CALL BUG (4HTABL,640,Z(K1),K2-K1+1)
C
C     CONSTRUCTION OF TABLE-ID LIST.
C
      ITABID = NLDID + 1
      NTABID = NLDID + 1
C
C     FIRST GET TABLE ID-S PRESENT IN THE SUB-ID TABLE.
C
      DO 680 I = K1,K2,11
C
C     CHECK FOR OTHER THAN TLOAD1 TYPE CARD
C
      IF (Z(I+4) .NE. 3) GO TO 680
C
C     CHECK FOR ID IN TABLE.
C
      IF (NTABID .LE. ITABID) GO TO 660
      DO 650 J = ITABID,NTABID
      IF (Z(I+5) .EQ. Z(J)) GO TO 680
  650 CONTINUE
  660 NTABID = NTABID + 1
      IF (NTABID .GT. CORE) CALL MESAGE (-8,0,SUBR)
      Z(NTABID) = Z(I+5)
  680 CONTINUE
C
C     NOW PASS SLT AND GET ANY TABLE IDS PRESENT IN QVECT PORTION OF
C     RECORDS WE WILL BE USING.  (SLT IS CURRENTLY POSITIONED AT FIRST
C     RECORD.)
C
      SLTAT = 1
      FILE  = SLT
      DO 780 I = K1,K2,11
      NGO = Z(I+2) - SLTAT
      IF (NGO) 780,710,690
  690 DO 700 J = 1,NGO
      CALL FWDREC (*1260,SLT)
  700 CONTINUE
      SLTAT = SLTAT + NGO
C
C     LOOK FOR QVECT CARDS.
C
  710 CALL READ (*1260,*770,SLT,BUF,2,NOEOR,IFLAG)
      ITYPE  = BUF(1)
      NCARDS = BUF(2)
      NWORDS = LDWORD(ITYPE)
      IF (ITYPE .NE. 16) GO TO 760
C
C     QVECT CARDS FOUND
C
      IF (NCARDS .LE. 0) GO TO 710
      DO 750 J = 1,NCARDS
      CALL READ (*1260,*1270,SLT,BUF,NWORDS,NOEOR,IFLAG)
      DO 740 K = 2,4
      L = NUMTYP(BUF(K))
      IF (MCH521 .AND. BUF(K).GT.16000 .AND. BUF(K).LE.99999999) L= 1
      IF (BUF(K).LE.0 .OR. L.NE.1) GO TO 740
C
C     TABLE ID FOUND.  ADD TO LIST IF NOT YET IN.
C
      IF (NTABID .LE. ITABID) GO TO 730
      DO 720 L = ITABID,NTABID
      IF (BUF(K) .EQ. Z(L)) GO TO 740
  720 CONTINUE
  730 NTABID = NTABID + 1
      IF (NTABID .GT. CORE) CALL MESAGE (-8,0,SUBR)
      Z(NTABID) = BUF(K)
  740 CONTINUE
  750 CONTINUE
      GO TO 710
  760 IF (ITYPE .GT. 16) GO TO 1170
      NWDCRD = -NWORDS*NCARDS
      CALL READ (*1260,*1270,SLT,BUF,NWDCRD,NOEOR,IFLAG)
      GO TO 710
  770 SLTAT = SLTAT + 1
  780 CONTINUE
      NUMTAB   = NTABID - ITABID
      Z(ITABID)= NUMTAB
      NUMTAB   = Z(ITABID)
C
C     TABLE-ID LIST COMPLETE. NOW SORT IT AND PRIME TAB ROUTINE.
C
      CALL REWIND (SLT)
      CALL FWDREC (*1260,SLT)
C/////
C     CALL BUG (4HTBID,555,Z(ITABID),NTABID-ITABID+1)
C/////
      IDIT = NTABID + 1
      NDIT = IDIT
      LZ   = CORE - IDIT
      IF (LZ .GT. 10) GO TO 790
      CALL MESAGE (8,0,SUBR)
      GO TO 1220
  790 CONTINUE
      IF (NUMTAB .EQ. 0) GO TO 792
      CALL SORT (0,0,1,1,Z(ITABID+1),NUMTAB)
      CALL PRETAB(DIT,Z(IDIT),Z(IDIT),Z(BUF2),LZ,LUSED,Z(ITABID),TABLST)
      NDIT = IDIT + LUSED
C/////
C     CALL BUG (4HDITS,557,Z(IDIT),NDIT-IDIT+1)
C/////
  792 CONTINUE
      IDREC = NDIT + 1
      GO TO 810
C
C     DETERMINE IF -LOADID- IS IN LIST OF LOAD SET IDS.
C
  800 IDREC = NLDID + 1
C/////
C     CALL BUG (4HLD1   ,360,Z(ILDID),NLDID-ILDID+1)
C/////
  810 CONTINUE
      NMAST = NLDID
      J = ILDID
  820 IF (J .LE. NLDID) GO TO 910
C
C     LOAD SET ID LIST EXHAUSTED.
C     BRING IN ANY LOAD CARDS IF NOT YET IN.
C
      IF (CARDIN) GO TO 920
C
C     THE LOAD-SET-ID TABLE HAS THE FOLLOWING FORMAT.
C
C     MASTER ID                                ******       Z(ILDID)
C     NUMBER OF SUB-IDS FOR THIS MASTER              *
C     SUB-ID             **  3-WORDS  ***             *
C     SCALE FACTOR = F(T)  * ONLY IF     *             *
C     SLT RECORD NUMBER  **  STATICS      * 11 WORDS   * REPEATS FOR
C     CONSTANT SCALE FACTOR               * REPEATS    * EACH MASTER
C     TYPE OF TLOAD =(3 OR 4)             * FOR EACH   * ID PRESENT
C     TYPE3 TABLE ID (OR) TYPE4 T1        * SUB-ID     *
C           0                   T2        * OF THIS    *
C           0                   OMEGA     * MASTER     *
C           0                   PHI       * ID         *
C           0                   N        *             *
C           0                   ALPHA ***              *
C              .                                       *
C              .                                      *
C              .                                     *
C                                              ******
C             ...                    ...
C             ...                    ...
C             ...                    ...                    Z(NLDID)
C
C
      CARDIN = .TRUE.
C
C     FORWARD SLT TO LOAD CARD RECORD.
C
      IF (NSETS) 850,850,830
  830 DO 840 I = 1,NSETS
      CALL FWDREC (*1250,SLT)
  840 CONTINUE
C
C     READ AND ENTER MASTER ID INTO TABLE
C
  850 IF (NLDID+2 .GT. CORE) CALL MESAGE (-8,0,SUBR)
      CALL READ (*1260,*900,SLT,Z(NLDID+1),2,NOEOR,IFLAG)
      SCALE  = RZ(NLDID+2)
      NLDID  = NLDID + 2
      JCOUNT = NLDID
      Z(JCOUNT) = 0
C
C     READ THE (SID, SCALE-FACTOR)  PAIRS FOR THIS ID.
C
  860 IF (NLDID+3 .GT. CORE) CALL MESAGE (-8,0,SUBR)
      CALL READ (*1260,*1270,SLT,Z(NLDID+1),2,NOEOR,IFLAG)
      IF (Z(NLDID+1) .EQ. -1) GO TO 890
C
C     MULTIPLY SUBID SCALE FACTOR BY MASTER SCALE FACTOR.
C
      RZ(NLDID+2) = RZ(NLDID+2)*SCALE
C
C     DETERMIND SLT RECORD NUMBER OF THIS SUB ID.
C
      KREC = 1
      DO 870 I = IMAST,NMAST,NLDSET
      IF (Z(NLDID+1) .EQ. Z(I)) GO TO 880
      KREC = KREC + 1
  870 CONTINUE
      KREC = 0
  880 Z(NLDID+3) = KREC
      NLDID = NLDID + 3
      Z(JCOUNT) = Z(JCOUNT) + 1
      GO TO 860
C
C     SORT ALL SUB-ID 3 WORD GROUPS ON SLT RECORD NUMBER.
C
  890 CALL SORT (0,0,3,3,Z(JCOUNT+1),NLDID-JCOUNT)
      GO TO 850
C
C     REPOSITION SLT TO BEGINNING OF FIRST SLT RECORD
C
  900 IDREC = NLDID + 1
C/////
C     CALL BUG (4HLDID,460,Z(ILDID),NLDID - ILDID+1)
C/////
      CALL REWIND (SLT)
      CALL FWDREC (*1260,SLT)
      GO TO 820
C
C     CONTINUE SEARCH FOR LOADID
C
  910 IF (LOADID .EQ. Z(J)) GO TO 940
C
C     POSITION -J- TO NEXT LOAD-SET-ID IN TABLE
C
      J = J + NLDSET*Z(J+1) + 2
      GO TO 820
C
C     -LOADID- NOT FOUND ANYWHERE.
C
  920 WRITE  (OUTPT,930) UWM,LOADID
  930 FORMAT (A25,' 3067, LOAD SET ID =',I9,' IS NOT PRESENT.')
      GO TO 1220
C
C     MATCH ON MASTER ID HAS BEEN FOUND.
C
  940 NLOADS = Z(J+1)
      ILOAD  = J + 2
      NLOAD  = ILOAD + NLDSET*NLOADS - 1
C
C     PROCESS ALL THE LOAD RECORDS FOR THIS MASTER-ID
C
      SLTAT = 1
C/////
C     CALL BUG (4HLOAD ,500,Z(ILOAD),NLOAD-ILOAD+1)
C/////
C
C     INITIALIZE APPLIED LOAD TO 0.0 FOR ALL ELEMENTS IN TABLE
C
 944  CONTINUE
      DO 945 I = IELTAB,NELTAB,LENTRY
      RZ(I+3) = 0.0
  945 CONTINUE
      IF (LOADID .LE. 0) GO TO 1140
      DO 1130 I = ILOAD,NLOAD,NLDSET
      FACTOR = RZ(I+1)
      IF (.NOT.TRANST) GO TO 960
C
C     FACTOR HAS TO BE FOUND AS F(TIME)
C
      IF (Z(I+4) .EQ. 4) GO TO 950
      CALL TAB (Z(I+5),TIME,YVALUE)
      FACTOR = RZ(I+3)*YVALUE
      GO TO 960
  950 TT = TIME - RZ(I+5)
      IF (TT .EQ. 0.0) GO TO 951
      IF (TT.LE.0.0 .OR. TIME.GE.RZ(I+6)) GO TO 955
      FACTOR = RZ(I+3)*EXP(RZ(I+10)*TT)*(TT**RZ(I+9))*COS(TWOPI*
     1         RZ(I+7)*TT + RZ(I+8)/RADDEG)
      GO TO 960
  951 IF (RZ(I+9) .NE. 0.0) GO TO 955
      FACTOR = COS(TWOPI*RZ(I+7))
      GO TO 960
  955 FACTOR = 0.0
  960 SLTREC = Z(I+2)
      IF (SLTREC.LE.0 .OR. FACTOR.EQ.0.0) GO TO 1130
C
C     POSITION SLT TO RECORD DESIRED.
C
  980 NGO = SLTREC - SLTAT
      IF (NGO) 990,1020,1000
C
C     NEED TO BACK UP ON SLT.
C
  990 CALL BCKREC (SLT)
      SLTAT = SLTAT - 1
      GO TO 980
C
C     NEED TO GO FORWARD ON SLT
C
 1000 DO 1010 J = 1,NGO
      CALL FWDREC (*1260,SLT)
 1010 CONTINUE
      SLTAT = SLTAT + NGO
C
C     SLT IS NOW POSITIONED TO LOAD RECORD DESIRED.
C
C
C     GENERATE LOADS FOR THOSE ELEMENTS IN THE TABLE USING ONLY QBDY1,
C     QBDY2, AND QVECT CARDS.
C
 1020 CALL READ (*1260,*1130,SLT,BUF,2,NOEOR,IFLAG)
      ITYPE = BUF(1)
      IF (ITYPE .LE. SLTYPS) GO TO 1040
      WRITE  (OUTPT,1030) SWM,ITYPE
 1030 FORMAT (A27,' 3068, UNRECOGNIZED CARD TYPE =',I9,
     1       ' FOUND IN -SLT- DATA BLOCK.')
      GO TO 1220
 1040 NCARDS = BUF(2)
      IF (NCARDS) 1020,1020,1050
 1050 NWORDS = LDWORD(ITYPE)
      IF (ITYPE.GE.14 .AND. ITYPE.LE.16) GO TO 1060
      NWDCRD = -NWORDS*NCARDS
      CALL READ (*1260,*1270,SLT,BUF,NWDCRD,NOEOR,IFLAG)
      GO TO 1020
 1060 ITYPE = ITYPE - 13
      JID = IDPOS(ITYPE)
      DO 1120 K = 1,NCARDS
C
C     READ A QBDY1, QBDY2, OR QVECT ENTRY.
C
      CALL READ (*1260,*1270,SLT,BUF,NWORDS,NOEOR,IFLAG)
C
C     CHECK FOR ID IN THE TABLE (OTHERWISE SKIP).
C
      CALL BISLOC (*1120,BUF(JID),Z(IELTAB),LENTRY,NUMBER,JPOINT)
      KK = IELTAB + JPOINT
C
C     THIS ELEMENT IS IN TABLE, THUS COMPUTE AND SUM IN THE LOAD.
C
      GO TO (1070,1080,1090), ITYPE
 1070 RZ(KK+2) = RZ(KK+2) + RZ(KK+4)*RBUF(1)*FACTOR
      GO TO 1120
C
 1080 KTYPE    = Z(KK+3)
      RZ(KK+2) = RZ(KK+2) + FACTOR*RZ(KK+4)*(RBUF(2)+RBUF(3) + RBUF(4) +
     1           RBUF(5))/GRIDS(KTYPE)
      GO TO 1120
C
C
C     CALL TAB IF E1,E2,E3 OF QVECT DATA ARE TABLE ID-S IMPLYING
C     TIME DEPENDENCE
C
 1090 IF (.NOT.TRANST) GO TO 1099
      DO 1094 KKK = 2,4
      L = NUMTYP(BUF(KKK))
      IF (MCH521 .AND. BUF(KKK).GT.16000 .AND. BUF(KKK).LE.99999999)
     1    L = 1
      IF (BUF(KKK).LE.0 .OR. L.NE.1) GO TO 1094
      CALL TAB (BUF(KKK),TIME,YVALUE)
      RBUF(KKK) = YVALUE
 1094 CONTINUE
 1099 KTYPE = Z(KK+3)
      C = RBUF(2)*RZ(KK+6) + RBUF(3)*RZ(KK+7) + RBUF(4)*RZ(KK+8)
      IF (KTYPE .EQ. 6) GO TO 1110
      IF (C) 1100,1100,1120
 1100 RZ(KK+2) = RZ(KK+2) - C*RZ(KK+4)*RZ(KK+5)*RBUF(1)*FACTOR
      GO TO 1120
 1110 RZ(KK+2) = RZ(KK+2) + FACTOR*RZ(KK+4)*RBUF(1)*RZ(KK+5)*SQRT(C*C +
     1  (RBUF(2)*RZ(KK+9) + RBUF(3)*RZ(KK+10) + RBUF(4)*RZ(KK+11))**2)
C
 1120 CONTINUE
C
 1130 CONTINUE
      CALL CLOSE (SLT,CLSREW)
C/////
C     CALL BUG (4HTELT,670,Z(IELTAB),NELTAB-IELTAB+1)
C/////
C
C     ELEMENT TABLE IS NOW COMPLETE FOR OUTPUT.
C
 1140 FILE = OEF1X
      CALL OPEN (*1250,OEF1X,Z(BUF1),WRT)
      DO 1160 I = IELTAB,NELTAB,LENTRY
      BUF( 1) = Z(I+13)
      RBUF(2) = RZ(I+3)
      RBUF(3) = RZ(I+2)
      RBUF(4) = RZ(I+1)
      RBUF(5) = RBUF(2) + RBUF(3) + RBUF(4)
      CALL WRITE (OEF1X,BUF(1),5,NOEOR)
 1160 CONTINUE
      CALL WRITE (OEF1X,0,0,EOR)
      FILE = OEF1
      CALL OPEN (*1250,OEF1,Z(BUF2),RD)
      GO TO 10
 1170 WRITE (OUTPT,1030)ITYPE
      GO TO 1220
C
C     ALL PROCESSING COMPLETE.
C
 1180 MCB(1) = OEF1
      CALL RDTRL (MCB)
      MCB(1) = OEF1X
      CALL WRTTRL (MCB)
      GO TO 1220
C
C     ERROR CONDITIONS.
C
 1190 RETURN
C
 1200 WRITE  (OUTPT,1210) UWM
 1210 FORMAT (A25,' 3069, OUTPUT DATA BLOCK FOR FORCES IS PURGED.')
 1220 CALL CLOSE (OEF1,CLSREW)
      CALL CLOSE (OEF1X,CLSREW)
      CALL CLOSE (UG,CLSREW)
      CALL CLOSE (EST,CLSREW)
      CALL CLOSE (SLT,CLSREW)
      CALL CLOSE (DLT,CLSREW)
      GO TO 1190
 1250 N = 1
      GO TO 1280
 1260 N = 2
      GO TO 1280
 1270 N = 3
 1280 CALL MESAGE (N,FILE,SUBR)
      GO TO 1220
      END


text:
! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!> Module eccodes\n
!> This is the \b Fortran90 interface for ecCodes
module eccodes

  use grib_api
  implicit none
  include "eccodes_constants.h"
  include "eccodes_visibility.h"

  real(8),    parameter,public :: CODES_MISSING_DOUBLE = -1.D+100
  integer(4), parameter,public :: CODES_MISSING_LONG   = 2147483647

  integer, parameter,public :: CODES_PRODUCT_ANY  = 0
  integer, parameter,public :: CODES_PRODUCT_GRIB = 1
  integer, parameter,public :: CODES_PRODUCT_BUFR = 2

  !> Create a new message in memory from an integer or character array containting the coded message.
  !>
  !> The message can be accessed through its ID and it will be available\n
  !> until @ref codes_release is called. A reference to the original coded\n
  !> message is kept in the new message structure.
  !>
  !> In case of error, if the status parameter (optional) is not given, the program will
  !> exit with an error message.\n Otherwise the error message can be
  !> gathered with @ref codes_get_error_string.
  !>
  !>
  !> \b Examples: \ref grib_copy_message.f90 "grib_copy_message.f90"
  !>
  !> @param id          ID of the message loaded in memory
  !> @param message     array containing the coded message
  !> @param status      CODES_SUCCESS if OK, integer value on error
  interface codes_new_from_message
      module procedure codes_new_from_message_int4
      module procedure codes_new_from_message_char
  end interface codes_new_from_message

  !> Get a value of specified index from an array key.
  !>
  !> Given an ID and key name as input a value corresponding to the given index
  !> is returned. The index is zero based i.e. the first element has
  !> zero index, the second element index one and so on.
  !> If the parameter index is an array all the values correspondig to the indexes
  !> list is returned.
  !> The ID references to a message loaded in memory.
  !>
  !> In case of error, if the status parameter (optional) is not given, the program will
  !> exit with an error message.\n Otherwise the error message can be
  !> gathered with @ref codes_get_error_string.
  !>
  !> \b Examples: \ref grib_nearest.f90 "grib_nearest.f90"
  !>
  !> @see codes_new_from_file, codes_release, codes_get
  !>
  !> @param[in] id          ID of the message loaded in memory
  !> @param[in] key         key name
  !> @param[in] index       index can be a scalar or array of integer(4)
  !> @param[out] value      value can be a scalar or array of integer(4),real(4),real(8)
  !> @param[out] status     CODES_SUCCESS if OK, integer value on error
  interface codes_get_element
    module procedure codes_get_real4_element, &
                     codes_get_real8_element, &
                     codes_get_real4_elements, &
                     codes_get_real8_elements
  end interface codes_get_element

  !> Find the nearest point/points of a given latitude/longitude point.
  !>
  !> The value in the nearest point (or the four nearest points) is returned as well as the
  !> zero based index (which can be used in @ref codes_get_element)
  !> and its distance from the given point using the following
  !> formula radius * acos( sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2) ).
  !>
  !> If the is_lsm flag is .true. the input field gribid is considered as
  !> a land sea mask and the nearest land point is returned.\n
  !> The nearest land point among the four neighbours is:
  !> - the nearest point with land sea mask value >= 0.5
  !> - the nearest without any other condition if all the four have land sea mask value <0.5.
  !>
  !> Arrays (real(8)) of latitude/longitude can be provided to find with one call
  !> the values,indexes and distances for all the lat/lon points listed in the arrays.
  !>
  !> If a single latitude/longitude point is provided and outlat,outlon,value,distance,index
  !> are defined as arrays with four elements the lat/lon coordinates and values, distances
  !> and indexes of the four nearest points are returned.
  !>
  !> In case of error, if the status parameter (optional) is not given, the program will
  !> exit with an error message.\n Otherwise the error message can be
  !> gathered with @ref codes_get_error_string.
  !>
  !> \b Examples: \ref grib_nearest.f90 "grib_nearest.f90"
  !>
  !> @param[in] gribid     id of the grib loaded in memory
  !> @param[in] is_lsm      .true. if the nearest land point is required otherwise .false.
  !> @param[in] inlat       latitude of the point in degrees
  !> @param[in] inlon       longitudes of the point in degrees
  !> @param[out] outlat     latitude of the nearest point in degrees
  !> @param[out] outlon     longitude of the nearest point in degrees
  !> @param[out] distance   distance between the given point and its nearest (in km)
  !> @param[out] index      zero based index
  !> @param[out] value      value of the field in the nearest point
  !> @param[out] status     CODES_SUCCESS if OK, integer value on error
  interface codes_grib_find_nearest
    module procedure codes_grib_find_nearest_single, &
                     codes_grib_find_nearest_four_single, &
                     codes_grib_find_nearest_multiple
  end interface codes_grib_find_nearest

  !> Get latitude/longitude and data values.
  !>
  !> Latitudes, longitudes, data values arrays are returned.
  !> They must be properly allocated by the caller and their required
  !> dimension can be obtained with \ref codes_get_size or by getting (with \ref codes_get)
  !> the value of the integer key "numberOfPoints".
  !>
  !> In case of error, if the status parameter (optional) is not given, the program will
  !> exit with an error message.\n Otherwise the error message can be
  !> gathered with @ref codes_get_error_string.
  !>
  !> \b Examples: \ref grib_get_data.f90 "grib_get_data.f90"
  !>
  !> @param[in] gribid       id of the grib loaded in memory
  !> @param[out] lats        latitudes array with dimension "size"
  !> @param[out] lons        longitudes array with dimension "size"
  !> @param[out] values      data values array with dimension "size"
  interface codes_grib_get_data
    module procedure codes_grib_get_data_real4, &
                     codes_grib_get_data_real8
  end interface codes_grib_get_data

max_stars_repo_name: mkrupcale/eccodes

max_issues_repo_name: mkrupcale/eccodes

max_forks_repo_name: mkrupcale/eccodes


!***********************************/
!	Name:途中から再計算を行うためのデータを吐き出すプログラム
!	Alias:Output4Resume
!	Description:プログラムと同じディレクトリにResume_***th_Resultという名前のファイルが生成される
!	Type:Geom,iStep
!	Input:Geom,iStep
!	Output:
!	Note:1次元にしか対応していない
!	Author:Akitaka Toyota
!	Date:2017.10.26
!	Update:2017.11.09
!	Other:
!***********************************/
subroutine Output4Resume(Geom,iStep)
    use StructVar_Mod
    use LoopVar_Mod
    implicit none

    type(Geometry),intent(in) :: Geom
    integer,intent(in) :: iStep

    integer :: iProgress
    double precision :: Tmp
    double precision, allocatable :: TmpQuantity(:)
    character(len=64), target :: cStep,cDensity,cVelocity,cPressure,cFileName
!    call Conserve2Primitive(Geom,CC)

    allocate(TmpQuantity(Geom%CellNumber(1)))
    write(cStep,*) iStep
    iProgress = int(dble(iStep)/10.0d0)

    cFileName = "Resume_"//trim(adjustl(cStep))//"th_Result"

    open(unit=2,file=trim(adjustl(cFileName)),status='unknown')
            write(cDensity, '("Density/Density",i6.6,".txt")') iStep
            open(unit=1,file=cDensity,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0 * Geom%Width(1,1,1) * dble(iCenterX),TmpQuantity(iCenterX)
            end do
            write(2,*) ""

            write(cVelocity, '("Velocity/Velocity",i6.6,".txt")') iStep
            open(unit=1,file=cVelocity,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0*Geom%Width(1,1,1)*dble(iCenterX), TmpQuantity(iCenterX)
            end do
            write(2,*) ""

            write(cPressure, '("Pressure/Pressure",i6.6,".txt")') iStep
            open(unit=1,file=cPressure,status='unknown')
                do iCenterX=1, Geom%CellNumber(1)
                    read(1,*) Tmp,TmpQuantity(iCenterX)
                end do
            close(1)

            do iCenterX=1, Geom%CellNumber(1)
                write(2,*) Geom%Bound(1,1) + 2.0d0*Geom%Width(1,1,1)*dble(iCenterX), TmpQuantity(iCenterX)
            end do
    close(2)
return
end subroutine Output4Resume


text:
Jump to Timeline #Navigation Navigation

  The Mauve Decade, or Gay Nineties started out with a major Earthquakes earthquake that shattered Winters.  The earthquake measured at a magnitude of 6.6.
  The Davisville Advertiser newspaper is briefly revived in 1892.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
The Davis City of Davis will give each household (not apartment dwellers) a compost bin for $10, along with great literature on how to compost and recycle.

The compost bin used to be given for free, but as of July 2010 there is a $10 fee are in place to help defray some of the cost of these items. The backyard compost bins have almost doubled in price within the past few years.  Due to budget issues, the City Recycling Program is no longer able to cover the entire cost of the compost bins. The $10 fee for compost bins reduces the cost and enables the City to continue distributing these compost bins to residents.

The bin consists of a 5 x 3 sheet of heavy black plastic with holes in it so the organic material can drain, breathe and decompose. The plastic sheet comes rolled up, so it is easy to carry on a bicycles bike. When you set it up in your yard, you can make the bin as small as you need simply by overlapping more of the plastic.

According to http://www.davisrecycling.org DavisRecycling: The City of Davis offers a yearround compost correspondence course to Davis residents in singlefamily homes. To enroll, call the Public works Department (7575686) to request a compost packet mailed to your home. Once you receive the packet, read the materials inside, fill out the enclosed quiz and return the quiz to the Public Works Department. A city staffperson will contact you to let you know that your quiz was received and that your compost bin is available for pickup. A $10 fee will be charged for a compost bin.

Visions of the Bin

Actually, these pictures show the old compost bins that the City used to give out. These compost bins were discontinued by the manufacturer and are no longer available. The ones currently given out are very similar, but with smaller holes. They are also made from recycled plastic.

20100424 15:12:24 nbsp Any idea how I can get one of these? Users/JenniferCook

20100424 15:48:11 nbsp Id start with the Community Services Department.  If they dont do it themselves, Im sure they can point you in the right direction.  Be sure to post on here when you get details! Users/TomGarberson

20100424 16:01:26 nbsp It would have helped for me to do 5 seconds of research before I posted this comment. The compost page has the answer! Users/JenniferCook


text:
!***************************************************************************************************
!>
!  Module for [[LSQR]].
!
!### History
!  * Jacob Williams : 8 Nov 2019 : created module
   module lsqr_module

   use lsqr_kinds
   use lsqpblas_module

   implicit none

   private

   type,abstract,public :: lsqr_solver
      !! main class to access the [[lsqr]] solver.
      !!
      !! You can use this class directory by extending it
      !! and specifying `aprod`, or you can use the
      !! [[lsqr_solver_ez]] class that has an easier
      !! interface.
      private
   contains
      private
      procedure(aprod_func),deferred,public :: aprod !! User function to access the sparse matrix `A`.
      procedure,public :: lsqr  !! main solver routine
      procedure,public :: acheck
      procedure,public :: xcheck
   end type lsqr_solver

   type,public,extends(lsqr_solver) :: lsqr_solver_ez
      !! a simplier version of [[lsqr_solver]] where
      !! the `aprod` function is provided internally.
      !! To use, first call the `initialize` method
      !! to set the matrix and other inputs.
      private

      integer :: m = 0 !! number of rows in `A` matrix
      integer :: n = 0 !! number of columns in `A` matrix
      integer :: num_nonzero_elements = 0 !! number of nonzero elements in `A` matrix
      integer,dimension(:),allocatable  :: irow !! sparsity row indices
      integer,dimension(:),allocatable  :: icol !! sparsity column indices
      real(wp),dimension(:),allocatable :: a    !! sparse `A` matrix

      real(wp) :: atol    = zero  !! relative error in definition of `A`
      real(wp) :: btol    = zero  !! relative error in definition of `b`
      real(wp) :: conlim  = zero  !! An upper limit on `cond(Abar)`, the apparent
                                  !! condition number of the matrix `Abar`.
      integer  :: itnlim  = 100   !! max iterations
      integer  :: nout    = 0     !! output unit for printing

      ! used in aprod_ez:
      real(wp),dimension(:),allocatable :: Ax   !! `A*x` (dimension m)
      real(wp),dimension(:),allocatable :: Aty  !! `A(transpose)*y` (dimension n)

      real(wp),dimension(:),allocatable :: v  !! workspace array (dimension n)
      real(wp),dimension(:),allocatable :: w  !! workspace array (dimension n)

   contains
      private
      procedure,public :: initialize => initialize_ez  !! Constructor. Must be call first.
      procedure,public :: solve => solve_ez
      procedure,public :: aprod => aprod_ez !! internal routine
   end type lsqr_solver_ez

   abstract interface
      subroutine aprod_func ( me, mode, m, n, x, y )
         !! User function to access the sparse matrix `A`.
         import :: wp, lsqr_solver
         implicit none
         class(lsqr_solver),intent(inout) :: me
         integer,intent(in) :: mode !! * If `mode = 1`, compute `y = y + A*x`.
                                    !!   `y` should be altered without changing x.
                                    !! * If `mode = 2`, compute `x = x + A(transpose)*y`.
                                    !!   `x` should be altered without changing `y`.
         integer,intent(in) :: m    !! number of rows in `A` matrix
         integer,intent(in) :: n    !! number of columns in `A` matrix
         real(wp),dimension(:),intent(inout) :: x
         real(wp),dimension(:),intent(inout) :: y
      end subroutine aprod_func
   end interface

   contains
!***************************************************************************************************

!*******************************************************************************
!>
!  Constructor for [[lsqr_solver_ez]].

   subroutine initialize_ez(me,m,n,a,irow,icol,atol,btol,conlim,itnlim,nout)

   implicit none

   class(lsqr_solver_ez),intent(out) :: me
   integer,intent(in)                :: m       !! number of rows in `A` matrix
   integer,intent(in)                :: n       !! number of columns in `A` matrix
   integer,dimension(:),intent(in)   :: irow    !! row indices of nonzero elements of `A`
   integer,dimension(:),intent(in)   :: icol    !! column indices of nonzero elements of `A`
   real(wp),dimension(:),intent(in)  :: a       !! nonzero elements of `A`
   real(wp),intent(in),optional      :: atol    !! relative error in definition of `A`
   real(wp),intent(in),optional      :: btol    !! relative error in definition of `b`
   real(wp),intent(in),optional      :: conlim  !! An upper limit on `cond(Abar)`, the apparent
                                                !! condition number of the matrix `Abar`.
   integer,intent(in),optional       :: itnlim  !! max iterations
   integer,intent(in),optional       :: nout    !! output unit for printing

   ! check for consistent inputs:
   if (any(size(a)/=[size(irow),size(icol)])) error stop 'invalid a,icol,irow sizes in initialize_ez'
   if (any(irow>m)) error stop 'invalid irow or m in initialize_ez'
   if (any(icol>n)) error stop 'invalid icol or n in initialize_ez'

   me%num_nonzero_elements = size(irow)
   me%m     = m
   me%n     = n
   me%irow  = irow
   me%icol  = icol
   me%a     = a

   ! optional inputs:
   if (present(atol))   me%atol   = atol
   if (present(btol))   me%btol   = btol
   if (present(conlim)) me%conlim = conlim
   if (present(itnlim)) me%itnlim = itnlim
   if (present(nout))   me%nout   = nout

   end subroutine initialize_ez
!*******************************************************************************

!*******************************************************************************
!>
!  The internal `aprod` function for the [[lsqr_solver_ez]] class.

    subroutine aprod_ez ( me, mode, m, n, x, y )

    implicit none

    class(lsqr_solver_ez),intent(inout) :: me
    integer,intent(in) :: mode  !! * If `mode = 1`, compute `y = y + A*x`.
                                !!   `y` should be altered without changing x.
                                !! * If `mode = 2`, compute `x = x + A(transpose)*y`.
                                !!   `x` should be altered without changing `y`.
    integer,intent(in) :: m     !! number of rows in `A` matrix
    integer,intent(in) :: n     !! number of columns in `A` matrix
    real(wp),dimension(:),intent(inout) :: x  !! [n]
    real(wp),dimension(:),intent(inout) :: y  !! [m]

    integer :: i    !! counter
    integer :: r    !! row index
    integer :: c    !! column index

    if (m/=me%m .or. n/=me%n) error stop 'lsqr_solver_ez class not properly initialized'

    select case (mode)

    case(1)    ! y = y + A*x

        !   A    x   Ax
        !  ---   -   -
        !  X0X   X   X
        !  0X0 * X = X
        !  00X   X   X
        !  00X       X

        ! A*x:
        if (.not. allocated(me%Ax)) allocate(me%Ax(me%m))
        me%Ax = zero
        do i = 1, me%num_nonzero_elements
            r = me%irow(i)
            c = me%icol(i)
            me%Ax(r) = me%Ax(r) + me%a(i)*x(c)
        end do

        y = y + me%Ax

    case(2)   ! x = x + A(transpose)*y

        !   A     Y   ATy
        !  ---    -   -
        !  X000   Y   X
        !  0X00 * Y = X
        !  X0XX   Y   X
        !         Y

        ! A(transpose)*y
        if (.not. allocated(me%Aty)) allocate(me%Aty(me%n))
        me%Aty = zero
        do i = 1, me%num_nonzero_elements
            r = me%irow(i)
            c = me%icol(i)
            me%Aty(c) = me%Aty(c) + me%a(i)*y(r)
        end do

        x = x + me%Aty

    case default
       error stop 'invalid mode in aprod_ez'
    end select

    end subroutine aprod_ez
!***************************************************************************************************

!***************************************************************************************************
!>
!  Wrapper for [[lsqr]] for the easy version of the class.

   subroutine solve_ez( me, b, damp, x, istop, &
                        se, itn, anorm, acond, rnorm, arnorm, xnorm)

   implicit none

   class(lsqr_solver_ez),intent(inout)   :: me
   real(wp),dimension(me%m),intent(in)   :: b
   real(wp),intent(in)                   :: damp
   real(wp),dimension(me%n),intent(out)  :: x       !! the computed solution `x`.
   integer,intent(out)                   :: istop   !! exit code (see [[lsqr]]).
   real(wp),dimension(me%n),intent(out),optional :: se
   integer,intent(out) ,optional         :: itn
   real(wp),intent(out),optional         :: anorm
   real(wp),intent(out),optional         :: acond
   real(wp),intent(out),optional         :: rnorm
   real(wp),intent(out),optional         :: arnorm
   real(wp),intent(out),optional         :: xnorm

   real(wp),dimension(:),allocatable :: u  !! copy of `b` for call to [[lsqr]]
   real(wp),dimension(:),allocatable :: se_
   logical  :: wantse   !! if `se` is to be returned
   integer  :: itn_
   real(wp) :: anorm_,acond_,rnorm_,arnorm_,xnorm_

   ! get optional inputs:
   wantse = present(se)
   if (wantse) then
      allocate(se_(me%n))
   else
      allocate(se_(1)) ! not needed
   end if
   allocate(u(me%m))
   if (.not. allocated(me%v)) allocate(me%v(me%n))
   if (.not. allocated(me%w)) allocate(me%w(me%n))

   u = b    ! make a copy for input to lsqr (since it will be overwritten)

   ! main routine:
   call me%lsqr(me%m, me%n, damp, wantse, &
                u, me%v, me%w, x, se_, &
                me%atol, me%btol, me%conlim, me%itnlim, me%nout, &
                istop, itn_, anorm_, acond_, rnorm_, arnorm_, xnorm_)

   ! optional outputs:
   if (wantse)          se     = se_
   if (present(itn))    itn    = itn_
   if (present(anorm))  anorm  = anorm_
   if (present(acond))  acond  = acond_
   if (present(rnorm))  rnorm  = rnorm_
   if (present(arnorm)) arnorm = arnorm_
   if (present(xnorm))  xnorm  = xnorm_

   end subroutine solve_ez
!***************************************************************************************************

!***************************************************************************************************
!>
!  LSQR finds a solution \(x\) to the following problems:
!
!  1. Unsymmetric equations -- solve:  $$ \mathbf{A} \cdot \mathbf{x} = \mathbf{b} $$
!  2. Linear least squares  -- solve (in the least-squares sense):  $$ \mathbf{A} \cdot \mathbf{x} = \mathbf{b} $$
!  3. Damped least squares  -- solve (in the least-squares sense): $$ \left[ \begin{array}{c} \mathbf{A}\\ damp \cdot \mathbf{I} \end{array} \right] \cdot \mathbf{x} = \left[ \begin{array}{c} \mathbf{b}\\ \mathbf{0} \end{array} \right] $$
!
!  where A is a matrix with m rows and n columns, b is an
!  m-vector, and `damp` is a scalar.  (All quantities are real.)
!  The matrix A is intended to be large and sparse.  It is accessed
!  by means of subroutine calls to `aprod`.
!
!  The rhs vector b is input via u, and subsequently overwritten.
!
!  Note:  LSQR uses an iterative method to approximate the solution.
!  The number of iterations required to reach a certain accuracy
!  depends strongly on the scaling of the problem.  Poor scaling of
!  the rows or columns of A should therefore be avoided where
!  possible.
!
!  For example, in problem 1 the solution is unaltered by
!  row-scaling.  If a row of A is very small or large compared to
!  the other rows of A, the corresponding row of ( A  b ) should be
!  scaled up or down.
!
!  In problems 1 and 2, the solution x is easily recovered
!  following column-scaling.  Unless better information is known,
!  the nonzero columns of A should be scaled so that they all have
!  the same Euclidean norm (e.g., 1.0).
!
!  In problem 3, there is no freedom to re-scale if damp is
!  nonzero.  However, the value of damp should be assigned only
!  after attention has been paid to the scaling of A.
!
!  The parameter damp is intended to help regularize
!  ill-conditioned systems, by preventing the true solution from
!  being very large.  Another aid to regularization is provided by
!  the parameter acond, which may be used to terminate iterations
!  before the computed solution becomes very large.
!
!  Note that x is not an input parameter.
!  If some initial estimate x0 is known and if damp = 0,
!  one could proceed as follows:
!
!  1. Compute a residual vector     \( r_0 = b - A \cdot x_0 \).
!  2. Use LSQR to solve the system  \( A \cdot \Delta x = r_0 \).
!  3. Add the correction dx to obtain a final solution \( x = x_0 + \Delta x \).
!
!  This requires that x0 be available before and after the call
!  to LSQR.  To judge the benefits, suppose LSQR takes k1 iterations
!  to solve A*x = b and k2 iterations to solve A*dx = r0.
!  If x0 is "good", norm(r0) will be smaller than norm(b).
!  If the same stopping tolerances atol and btol are used for each
!  system, k1 and k2 will be similar, but the final solution x0 + dx
!  should be more accurate.  The only way to reduce the total work
!  is to use a larger stopping tolerance for the second system.
!  If some value btol is suitable for A*x = b, the larger value
!  btol*norm(b)/norm(r0)  should be suitable for A*dx = r0.
!
!  Preconditioning is another way to reduce the number of iterations.
!  If it is possible to solve a related system M*x = b efficiently,
!  where M approximates A in some helpful way
!  (e.g. M - A has low rank or its elements are small relative to
!  those of A), LSQR may converge more rapidly on the system
!        A*M(inverse)*z = b,
!  after which x can be recovered by solving M*x = z.
!
!  NOTE: If A is symmetric, LSQR should not be used!
!  Alternatives are the symmetric conjugate-gradient method (cg)
!  and/or SYMMLQ.
!  SYMMLQ is an implementation of symmetric cg that applies to
!  any symmetric A and will converge more rapidly than LSQR.
!  If A is positive definite, there are other implementations of
!  symmetric cg that require slightly less work per iteration
!  than SYMMLQ (but will take the same number of iterations).
!
!### Notation
!
! The following quantities are used in discussing the subroutine
! parameters:
!
!```
!     Abar   =  (   A    ),         bbar  =  ( b )
!               ( damp*I )                   ( 0 )
!
!     r      = b  -  A*x,           rbar  = bbar  -  Abar*x
!
!     rnorm  = sqrt( norm(r)**2  +  damp**2 * norm(x)**2 )
!            = norm( rbar )
!
!     relpr  = the relative precision of floating-point arithmetic
!              on the machine being used.  On most machines,
!              relpr is about 1.0e-7 and 1.0d-16 in single and double
!              precision respectively.
!```
!
! LSQR  minimizes the function `rnorm` with respect to `x`.
!
!### References
!
!  * C.C. Paige and M.A. Saunders,  LSQR: An algorithm for sparse
!    linear equations and sparse least squares,
!    ACM Transactions on Mathematical Software 8, 1 (March 1982),
!    pp. 43-71.
!
!  * C.C. Paige and M.A. Saunders,  Algorithm 583, LSQR: Sparse
!    linear equations and least-squares problems,
!    ACM Transactions on Mathematical Software 8, 2 (June 1982),
!    pp. 195-209.
!
!  * C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh,
!    Basic linear algebra subprograms for Fortran usage,
!    ACM Transactions on Mathematical Software 5, 3 (Sept 1979),
!    pp. 308-323 and 324-325.
!
!### LSQR development
!
!  * 22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583.
!  * 15 Sep 1985: Final F66 version.  LSQR sent to "misc" in netlib.
!  * 13 Oct 1987: Bug (Robert Davies, DSIR).  Have to delete
!                     if ( (one + dabs(t)) <= one ) GO TO 200
!                  from loop 200.  The test was an attempt to reduce
!                  underflows, but caused w(i) not to be updated.
!  * 17 Mar 1989: First F77 version.
!  * 04 May 1989: Bug (David Gay, AT&T).  When the second beta is zero,
!                  rnorm = 0 and
!                  test2 = arnorm / (anorm * rnorm) overflows.
!                  Fixed by testing for rnorm = 0.
!  * 05 May 1989: Sent to "misc" in netlib.
!  * 14 Mar 1990: Bug (John Tomlin via IBM OSL testing).
!                  Setting rhbar2 = rhobar**2 + dampsq can give zero
!                  if rhobar underflows and damp = 0.
!                  Fixed by testing for damp = 0 specially.
!  * 15 Mar 1990: Converted to lower case.
!  * 21 Mar 1990: d2norm introduced to avoid overflow in numerous
!                  items like  c = sqrt( a**2 + b**2 ).
!  * 04 Sep 1991: wantse added as an argument to LSQR, to make
!                  standard errors optional.  This saves storage and
!                  time when se(*) is not wanted.
!  * 13 Feb 1992: istop now returns a value in [1,5], not [1,7].
!                  1, 2 or 3 means that x solves one of the problems
!                  Ax = b,  min norm(Ax - b)  or  damped least squares.
!                  4 means the limit on cond(A) was reached.
!                  5 means the limit on iterations was reached.
!  * 07 Dec 1994: Keep track of dxmax = max_k norm( phi_k * d_k ).
!                  So far, this is just printed at the end.
!                  A large value (relative to norm(x)) indicates
!                  significant cancellation in forming
!                  x  = D*f  = sum( phi_k * d_k ).
!                  A large column of D need NOT be serious if the
!                  corresponding phi_k is small.
!  * 27 Dec 1994: Include estimate of alfa_opt in iteration log.
!                  alfa_opt is the optimal scale factor for the
!                  residual in the "augmented system", as described by
!                  A. Bjorck (1992),
!                  Pivoting and stability in the augmented system method,
!                  in D. F. Griffiths and G. A. Watson (eds.),
!                  "Numerical Analysis 1991",
!                  Proceedings of the 14th Dundee Conference,
!                  Pitman Research Notes in Mathematics 260,
!                  Longman Scientific and Technical, Harlow, Essex, 1992.
!  * 12 Nov 2019 : Jacob Williams : significant refactoring into modern Fortran.
!
!### Author
! * Michael A. Saunders, Dept of Operations Research, Stanford University
!
!@note The number of iterations required by LSQR will usually decrease
!      if the computation is performed in higher precision.

   subroutine LSQR  ( me, m, n, damp, wantse, &
                     u, v, w, x, se, &
                     atol, btol, conlim, itnlim, nout, &
                     istop, itn, anorm, acond, rnorm, arnorm, xnorm)

   class(lsqr_solver),intent(inout) :: me
   integer,intent(in)    :: m          !! the number of rows in `A`.
   integer,intent(in)    :: n          !! the number of columns in `A`.
   real(wp),intent(in)   :: damp       !! The damping parameter for problem 3 above.
                                       !! (damp should be 0.0 for problems 1 and 2.)
                                       !! If the system `A*x = b` is incompatible, values
                                       !! of `damp` in the range 0 to `sqrt(relpr)*norm(A)`
                                       !! will probably have a negligible effect.
                                       !! Larger values of `damp` will tend to decrease
                                       !! the norm of `x` and reduce the number of
                                       !! iterations required by LSQR.
                                       !!
                                       !! The work per iteration and the storage needed
                                       !! by LSQR are the same for all values of `damp`.
   logical,intent(in)    :: wantse     !! A logical variable to say if the array `se(*)`
                                       !! of standard error estimates should be computed.
                                       !! If `m > n`  or  `damp > 0`,  the system is
                                       !! overdetermined and the standard errors may be
                                       !! useful.  (See the first LSQR reference.)
                                       !! Otherwise (`m <= n`  and  `damp = 0`) they do not
                                       !! mean much.  Some time and storage can be saved
                                       !! by setting `wantse = .false.` and using any
                                       !! convenient array for `se(*)`, which won't be
                                       !! touched.
   real(wp),intent(inout):: u(m)       !! The rhs vector `b`.  Beware that `u` is
                                       !! over-written by LSQR.
   real(wp),intent(inout):: v(n)       !! workspace
   real(wp),intent(inout):: w(n)       !! workspace
   real(wp),intent(out)  :: x(n)       !! Returns the computed solution `x`.
   real(wp),dimension(*),intent(out) :: se   !! If `wantse` is true, the dimension of `se` must be
                                             !! `n` or more. `se(*)` then returns standard error
                                             !! estimates for the components of `x`.
                                             !! For each `i`, `se(i)` is set to the value
                                             !! `rnorm * sqrt( sigma(i,i) / t )`,
                                             !! where `sigma(i,i)` is an estimate of the i-th
                                             !! diagonal of the inverse of `Abar(transpose)*Abar`
                                             !! and:
                                             !! * `t = 1      if  m <= n`
                                             !! * `t = m - n  if  m > n  and  damp = 0`
                                             !! * `t = m      if  damp /= 0`
                                             !!
                                             !! If `wantse` is false, `se(*)` will not be touched.
                                             !! The actual parameter can be any suitable array
                                             !! of any length.
   real(wp),intent(in)   :: atol       !! An estimate of the relative error in the data
                                       !! defining the matrix `A`.  For example,
                                       !! if `A` is accurate to about 6 digits, set
                                       !! `atol = 1.0e-6`.
   real(wp),intent(in)   :: btol       !! An estimate of the relative error in the data
                                       !! defining the rhs vector `b`.  For example,
                                       !! if `b` is accurate to about 6 digits, set
                                       !! `btol = 1.0e-6`.
   real(wp),intent(in)              :: conlim   !! An upper limit on `cond(Abar)`, the apparent
                                                !! condition number of the matrix `Abar`.
                                                !! Iterations will be terminated if a computed
                                                !! estimate of `cond(Abar)` exceeds `conlim`.
                                                !! This is intended to prevent certain small or
                                                !! zero singular values of `A` or `Abar` from
                                                !! coming into effect and causing unwanted growth
                                                !! in the computed solution.
                                                !!
                                                !! `conlim` and `damp` may be used separately or
                                                !! together to regularize ill-conditioned systems.
                                                !!
                                                !! Normally, `conlim` should be in the range
                                                !! 1000 to `1/relpr`.
                                                !!
                                                !! Suggested value:
                                                !!
                                                !! * `conlim = 1/(100*relpr)` for compatible systems,
                                                !! * `conlim = 1/(10*sqrt(relpr))` for least squares.
                                                !!
                                                !! Note:  If the user is not concerned about the parameters
                                                !! `atol`, `btol` and `conlim`, any or all of them may be set
                                                !! to zero.  The effect will be the same as the values
                                                !! `relpr`, `relpr` and `1/relpr` respectively.
   integer,intent(in)               :: itnlim   !! An upper limit on the number of iterations.
                                                !! Suggested value:
                                                !! * `itnlim = n/2` for well-conditioned systems
                                                !!   with clustered singular values,
                                                !! * `itnlim = 4*n` otherwise.
   integer,intent(in)               :: nout     !! File number for printed output.  If nonzero,
                                                !! a summary will be printed on file `nout`.
   integer,intent(out)               :: istop   !! An integer giving the reason for termination:
                                                !!
                                                !! * 0 -- `x` = 0  is the exact solution.
                                                !!   No iterations were performed.
                                                !! * 1 -- The equations `A*x = b` are probably
                                                !!   compatible.  `Norm(A*x - b)` is sufficiently
                                                !!   small, given the values of `atol` and `btol`.
                                                !! * 2 -- `damp` is zero.  The system `A*x = b` is probably
                                                !!   not compatible.  A least-squares solution has
                                                !!   been obtained that is sufficiently accurate,
                                                !!   given the value of `atol`.
                                                !! * 3 -- `damp` is nonzero.  A damped least-squares
                                                !!   solution has been obtained that is sufficiently
                                                !!   accurate, given the value of `atol`.
                                                !! * 4 -- An estimate of `cond(Abar)` has exceeded
                                                !!   `conlim`.  The system `A*x = b` appears to be
                                                !!   ill-conditioned.  Otherwise, there could be an
                                                !!   error in subroutine `aprod`.
                                                !! * 5 -- The iteration limit `itnlim` was reached.
   integer,intent(out)   :: itn        !! The number of iterations performed.
   real(wp),intent(out)  :: anorm      !! An estimate of the Frobenius norm of `Abar`.
                                       !! This is the square-root of the sum of squares
                                       !! of the elements of `Abar`.
                                       !! If `damp` is small and if the columns of `A`
                                       !! have all been scaled to have length 1.0,
                                       !! `anorm` should increase to roughly `sqrt(n)`.
                                       !! A radically different value for `anorm` may
                                       !! indicate an error in subroutine `aprod` (there
                                       !! may be an inconsistency between modes 1 and 2).
   real(wp),intent(out)  :: acond      !! An estimate of `cond(Abar)`, the condition
                                       !! number of `Abar`.  A very high value of `acond`
                                       !! may again indicate an error in `aprod`.
   real(wp),intent(out)  :: rnorm      !! An estimate of the final value of `norm(rbar)`,
                                       !! the function being minimized (see notation
                                       !! above).  This will be small if `A*x = b` has
                                       !! a solution.
   real(wp),intent(out)  :: arnorm     !! An estimate of the final value of
                                       !! `norm( Abar(transpose)*rbar )`, the norm of
                                       !! the residual for the usual normal equations.
                                       !! This should be small in all cases.  (arnorm
                                       !! will often be smaller than the true value
                                       !! computed from the output vector `x`.)
   real(wp),intent(out)  :: xnorm      !! An estimate of the norm of the final
                                       !! solution vector `x`.

   logical :: damped
   integer :: i, maxdx, nconv, nstop
   real(wp) :: alfopt, alpha, beta, bnorm, &
               cs, cs1, cs2, ctol, &
               delta, dknorm, dnorm, dxk, dxmax, &
               gamma, gambar, phi, phibar, psi, &
               res2, rho, rhobar, rhbar1, &
               rhs, rtol, sn, sn1, sn2, &
               t, tau, temp, test1, test2, test3, &
               theta, t1, t2, t3, xnorm1, z, zbar
   logical :: print_iter

   logical,parameter :: extra = .true.  !! for extra printing below.

   character(len=*),parameter :: enter = ' Enter LSQR.  '
   character(len=*),parameter :: exit  = ' Exit  LSQR.  '
   character(len=*),dimension(0:5),parameter :: msg = [ 'The exact solution is x = 0                          ',&
                                                        'A solution to Ax = b was found, given atol, btol     ',&
                                                        'A least-squares solution was found, given atol       ',&
                                                        'A damped least-squares solution was found, given atol',&
                                                        'Cond(Abar) seems to be too large, given conlim       ',&
                                                        'The iteration limit was reached                      ' ]

   ! Initialize.
   if (nout /= 0) then
      write(nout,'(//A)') enter//'     Least-squares solution of  Ax = b'
      write(nout,'(A,I7,A,I7,A)') ' The matrix  A  has', m, ' rows   and', n, ' columns'
      write(nout,'(1P,A,E22.14,3X,A,L10)')   ' damp   =', damp, 'wantse =', wantse
      write(nout,'(1P,A,E10.2,15x,A,E10.2)') ' atol   =', atol, 'conlim =', conlim
      write(nout,'(1P,A,E10.2,15x,A,I10)')   ' btol   =', btol, 'itnlim =', itnlim
   end if

   damped = damp > zero
   itn  = 0
   istop = 0
   nstop = 0
   maxdx = 0
   if (conlim > zero) then
      ctol = one / conlim
   else
      ctol = zero
   end if
   anorm = zero
   acond = zero
   dnorm = zero
   dxmax = zero
   res2 = zero
   psi  = zero
   xnorm = zero
   xnorm1 = zero
   cs2    = - one
   sn2  = zero
   z    = zero

   ! Set up the first vectors u and v for the bidiagonalization.
   ! These satisfy  beta*u = b,  alpha*v = A(transpose)*u.
   do i = 1, n
      v(i)  = zero
      x(i)  = zero
   end do

   if ( wantse ) then
      do i = 1, n
         se(i) = zero
      end do
   end if

   alpha = zero
   beta = dnrm2 ( m, u, 1 )

   if (beta > zero) then
      call dscal ( m, (one / beta), u, 1 )
      call me%aprod ( 2, m, n, v, u )
      alpha = dnrm2 ( n, v, 1 )
   end if

   if (alpha > zero) then
      call dscal ( n, (one / alpha), v, 1 )
      call dcopy ( n, v, 1, w, 1 )
   end if

   arnorm = alpha * beta

   if (arnorm /= zero) then

      rhobar = alpha
      phibar = beta
      bnorm = beta
      rnorm = beta

      if (nout /= 0) then
         if ( damped ) then
            write(nout, '(//A)') &
               '   Itn       x(1)           Function     Compatible   LS     Norm Abar Cond Abar'
         else
            write(nout, '(//A)') &
               '   Itn       x(1)           Function     Compatible   LS        Norm A    Cond A'
         end if
         test1  = one
         test2  = alpha / beta

         if ( extra ) then
            write(nout, '(80X,A)') '    phi    dknorm   dxk  alfa_opt'
         end if
         write(nout, '(1P, I6, 2E17.9, 4E10.2, E9.1, 3E8.1)') itn, x(1), rnorm, test1, test2
         write(nout, '(A)') ''
      end if

      do
         ! Main iteration loop.
         itn = itn + 1

         ! Perform the next step of the bidiagonalization to obtain the
         ! next  beta, u, alpha, v.  These satisfy the relations
         ! beta*u  = A*v  -  alpha*u,
         ! alpha*v  = A(transpose)*u  -  beta*v.
         call dscal ( m, (- alpha), u, 1 )
         call me%aprod ( 1, m, n, v, u )
         beta = dnrm2 ( m, u, 1 )

         ! Accumulate  anorm = || Bk || = sqrt( sum of  alpha**2 + beta**2 + damp**2 ).

         temp = d2norm( alpha, beta )
         temp = d2norm( temp , damp )
         anorm = d2norm( anorm, temp )

         if (beta > zero) then
            call dscal ( m, (one / beta), u, 1 )
            call dscal ( n, (- beta), v, 1 )
            call me%aprod ( 2, m, n, v, u )
            alpha = dnrm2 ( n, v, 1 )
            if (alpha > zero) then
               call dscal ( n, (one / alpha), v, 1 )
            end if
         end if

         ! Use a plane rotation to eliminate the damping parameter.
         ! This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
         rhbar1 = rhobar
         if ( damped ) then
            rhbar1 = d2norm( rhobar, damp )
            cs1    = rhobar / rhbar1
            sn1    = damp   / rhbar1
            psi    = sn1 * phibar
            phibar = cs1 * phibar
         end if

         ! Use a plane rotation to eliminate the subdiagonal element (beta)
         ! of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
         rho  = d2norm( rhbar1, beta )
         cs   = rhbar1 / rho
         sn   = beta   / rho
         theta = sn * alpha
         rhobar = - cs * alpha
         phi  = cs * phibar
         phibar = sn * phibar
         tau  = sn * phi

         ! Update  x, w  and (perhaps) the standard error estimates.
         t1   = phi   / rho
         t2   = - theta / rho
         t3   = one   / rho
         dknorm = zero

         if ( wantse ) then
            do i = 1, n
               t      = w(i)
               x(i)   = t1*t  +  x(i)
               w(i)   = t2*t  +  v(i)
               t      = (t3*t)**2
               se(i)  = t     +  se(i)
               dknorm = t     +  dknorm
            end do
         else
            do i = 1, n
               t      = w(i)
               x(i)   = t1*t  +  x(i)
               w(i)   = t2*t  +  v(i)
               dknorm = (t3*t)**2  +  dknorm
            end do
         end if

         ! Monitor the norm of d_k, the update to x.
         ! dknorm = norm( d_k )
         ! dnorm  = norm( D_k ),        where   D_k = (d_1, d_2, ..., d_k )
         ! dxk    = norm( phi_k d_k ),  where new x = x_k + phi_k d_k.
         dknorm = sqrt( dknorm )
         dnorm  = d2norm( dnorm, dknorm )
         dxk    = abs( phi * dknorm )
         if (dxmax < dxk ) then
            dxmax   = dxk
            maxdx   = itn
         end if

         ! Use a plane rotation on the right to eliminate the
         ! super-diagonal element (theta) of the upper-bidiagonal matrix.
         ! Then use the result to estimate  norm(x).
         delta = sn2 * rho
         gambar = - cs2 * rho
         rhs  = phi    - delta * z
         zbar = rhs    / gambar
         xnorm = d2norm( xnorm1, zbar  )
         gamma = d2norm( gambar, theta )
         cs2  = gambar / gamma
         sn2  = theta  / gamma
         z    = rhs    / gamma
         xnorm1 = d2norm( xnorm1, z     )

         ! Test for convergence.
         ! First, estimate the norm and condition of the matrix  Abar,
         ! and the norms of  rbar  and  Abar(transpose)*rbar.
         acond = anorm * dnorm
         res2 = d2norm( res2 , psi    )
         rnorm = d2norm( res2 , phibar )
         arnorm = alpha * abs( tau )

         ! Now use these norms to estimate certain other quantities,
         ! some of which will be small near a solution.

         alfopt = sqrt( rnorm / (dnorm * xnorm) )
         test1 = rnorm /  bnorm
         test2 = zero
         if (rnorm > zero) test2 = arnorm / (anorm * rnorm)
         test3 = one   /  acond
         t1   = test1 / (one  +  anorm * xnorm / bnorm)
         rtol = btol  +  atol *  anorm * xnorm / bnorm

         ! The following tests guard against extremely small values of
         ! atol, btol  or  ctol.  (The user may have set any or all of
         ! the parameters  atol, btol, conlim  to zero.)
         ! The effect is equivalent to the normal tests using
         ! atol = relpr,  btol = relpr,  conlim = 1/relpr.

         t3   = one + test3
         t2   = one + test2
         t1   = one + t1
         if (itn >= itnlim) istop = 5
         if (t3  <= one   ) istop = 4
         if (t2  <= one   ) istop = 2
         if (t1  <= one   ) istop = 1

         ! Allow for tolerances set by the user.

         if (test3 <= ctol) istop = 4
         if (test2 <= atol) istop = 2
         if (test1 <= rtol) istop = 1

         ! See if it is time to print something.
         if (nout /= 0) then

            print_iter = (n     <= 40       ) .or. &
                         (itn   <= 10       ) .or. &
                         (itn   >= itnlim-10) .or. &
                         (mod(itn,10) == 0  ) .or. &
                         (test3 <=  2.0*ctol) .or. &
                         (test2 <= 10.0*atol) .or. &
                         (test1 <= 10.0*rtol) .or. &
                         (istop /=  0       )

            if (print_iter) then
               ! Print a line for this iteration.
               ! "extra" is for experimental purposes.
               if ( extra ) then
                  write(nout, '(1P, I6, 2E17.9, 4E10.2, E9.1, 3E8.1)') &
                           itn, x(1), rnorm, test1, test2, anorm, acond, phi, dknorm, dxk, alfopt
               else
                  write(nout, '(1P, I6, 2E17.9, 4E10.2, E9.1, 3E8.1)') &
                           itn, x(1), rnorm, test1, test2, anorm, acond
               end if
               !if (mod(itn,10) == 0) write(nout, '(A)') ''
            end if

         end if

         ! Stop if appropriate.
         ! The convergence criteria are required to be met on nconv
         ! consecutive iterations, where nconv is set below.
         ! Suggested value:  nconv = 1, 2 or 3.
         if (istop == 0) then
            nstop  = 0
         else
            nconv  = 1
            nstop  = nstop + 1
            if (nstop < nconv  .and.  itn < itnlim) istop = 0
         end if
         if (istop /= 0) exit

      end do
      ! End of iteration loop.

      ! Finish off the standard error estimates.

      if ( wantse ) then
         t = one
         if (m > n)  t = m - n
         if ( damped )  t = m
         t = rnorm / sqrt( t )
         do i = 1, n
            se(i) = t * sqrt( se(i) )
         end do
      end if

   end if

   ! Decide if istop = 2 or 3.
   ! Print the stopping condition.
   if (damped .and. istop == 2) istop = 3
   if (nout /= 0) then
      write(nout, '(//A,5X,A,I2,15X,A,I8)')       exit, 'istop  =', istop, 'itn    =', itn
      write(nout, '(1P,A,5X,A,E12.5,5X,A,E12.5)') exit, 'anorm  =', anorm, 'acond  =', acond
      write(nout, '(1P,A,5X,A,E12.5,5X,A,E12.5)') exit, 'bnorm  =', bnorm, 'xnorm  =', xnorm
      write(nout, '(1P,A,5X,A,E12.5,5X,A,E12.5)') exit, 'rnorm  =', rnorm, 'arnorm =', arnorm
      write(nout, '(1P,A,5X,A,E8.1,A,I8)')        exit, 'max dx =', dxmax, ' occurred at itn ',maxdx
      write(nout, '(1P,A,5X,A,E8.1,A)'   )        exit, '       =', dxmax/(xnorm + 1.0e-20_wp), '*xnorm'
      write(nout, '(A,5X,A)' )                    exit, msg(istop)
   end if

   end subroutine LSQR
!***************************************************************************************************

!***************************************************************************************************
!>
!  Checks the two modes of aprod for [[lsqr]].
!
!  acheck may be called to test the user-written subroutine
!  aprod required by LSQR and CRAIG.  For some m x n matrix A,
!  aprod with mode = 1 and 2 supplies LSQR and CRAIG with products
!  of the form
!     y := y + Ax  and  x := x + A'y
!  respectively, where A' means A(transpose).
!  acheck tries to verify that A and A' refer to the same matrix.
!
!### Method
!  We cook up some "unlikely" vectors x and y of unit length
!  and test if  y'(y + Ax)  =  x'(x + A'y).
!
!### History
! * 04 Sep 1991  Initial design and code.
!   Michael Saunders, Dept of Operations Research,
!   Stanford University
! * 10 Feb 1992  aprod and eps added as parameters.
!   tol defined via power.

   subroutine acheck( me, m, n, nout, eps, &
                     v, w, x, y, inform )

   implicit none

   class(lsqr_solver),intent(inout) :: me
   integer,intent(in)   :: m       !! No. of rows of A.
   integer,intent(in)   :: n       !! No. of columns of A.
   integer,intent(in)   :: nout    !! A file number for printed output.
   integer,intent(out)  :: inform  !! Error indicator.
                                   !! inform = 0 if aprod seems to be
                                   !! consistent.
                                   !! inform = 1 otherwise.
   real(wp),intent(in)  :: eps     !! The machine precision.
   real(wp)             :: v(n)
   real(wp)             :: w(m)
   real(wp)             :: x(n)
   real(wp)             :: y(m)

   real(wp),parameter :: power = 0.5_wp   !! eps**power is used as the tolerance for judging
                                          !! whether `y'(y + Ax)  =  x'(x + A'y)`
                                          !! to sufficient accuracy.
                                          !! power should be in the range (0.25, 0.9) say.
                                          !! For example, power = 0.75 means that we are happy
                                          !! if three quarters of the available digits agree.
                                          !! power = 0.5 seems a reasonable requirement
                                          !! (asking for half the digits to agree).

   integer :: i, j
   real(wp) :: alfa, beta, t, test1, test2, test3, tol

   tol = eps**power
   if (nout /= 0) write(nout, '(//A)') &
      'Enter acheck. Test of aprod for LSQR and CRAIG'

   ! ==================================================================
   ! Cook up some "unlikely" vectors x and y of unit length.
   ! ==================================================================
   t = one
   do j = 1, n
      t    = t + one
      x(j) = sqrt( t )
   end do

   t = one
   do i = 1, m
      t    = t + one
      y(i) = one / sqrt( t )
   end do

   alfa = dnrm2 ( n, x, 1 )
   beta = dnrm2 ( m, y, 1 )
   call dscal ( n, (one / alfa), x, 1 )
   call dscal ( m, (one / beta), y, 1 )

   ! ==================================================================
   ! Test if y'(y + Ax) = x'(x + A'y).
   ! ==================================================================

   ! First set w = y + Ax, v = x + A'y.

   call dcopy ( m, y, 1, w, 1 )
   call dcopy ( n, x, 1, v, 1 )
   call me%aprod ( 1, m, n, x, w )
   call me%aprod ( 2, m, n, v, y )

   ! Now set alfa = y'w, beta = x'v.

   alfa  = ddot  ( m, y, 1, w, 1 )
   beta  = ddot  ( n, x, 1, v, 1 )
   test1 = abs( alfa - beta )
   test2 = one  +  abs( alfa )  +  abs( beta )
   test3 = test1 / test2

   ! See if alfa and beta are essentially the same.

   if ( test3 <= tol ) then
      inform = 0
      if (nout /= 0) write(nout, '(1P,A,1X,E10.1)') &
         'aprod seems OK. Relative error =', test3
   else
      inform = 1
      if (nout /= 0) write(nout, '(1P,A,1X,E10.1)') &
         'aprod seems incorrect. Relative error =', test3
   end if

   end subroutine acheck
!***************************************************************************************************

!***************************************************************************************************
!>
!  Tests if `x` solves a certain least-squares problem.
!
!  xcheck computes residuals and norms associated with the
!  vector x and the least-squares problem solved by LSQR or CRAIG.
!  It determines whether x seems to be a solution to any of three
!  possible systems:
!
!  1.  Ax = b
!  2.  min norm(Ax - b)
!  3.  min norm(Ax - b)^2 + damp^2 * norm(x)^2
!
!### History
! * 07 Feb 1992  Initial design and code.
!   Michael Saunders, Dept of Operations Research,
!   Stanford University.

   subroutine xcheck( me, m, n, nout, anorm, damp, eps, &
                     b, u, v, w, x, &
                     inform, test1, test2, test3 )

   implicit none

   class(lsqr_solver),intent(inout) :: me
   integer,intent(in)   :: m  !! The number of rows in A.
   integer,intent(in)   :: n  !! The number of columns in A.
   integer,intent(in)   :: nout  !! A file number for printed output.
                                 !! If nout = 0, nothing is printed.
   integer,intent(out)  :: inform   !! inform = 0 if b = 0 and x = 0.
                                    !! inform = 1, 2 or 3 if x seems to
                                    !! solve systems 1 2 or 3 above.
   real(wp),intent(in)  :: anorm    !! An estimate of norm(A) or
                                    !! norm( A, delta*I ) if delta > 0.
                                    !! Normally this will be available
                                    !! from LSQR or CRAIG.
   real(wp),intent(in)  :: damp     !! Possibly defines a damped problem.
   real(wp),intent(in)  :: eps      !! Machine precision.
   real(wp),intent(out) :: test1    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(out) :: test2    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(out) :: test3    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(in)  :: b(m)     !! The right-hand side of Ax = b etc.
   real(wp),intent(out) :: u(m)     !! On exit, u = r (where r = b - Ax).
   real(wp),intent(out) :: v(n)     !! On exit, v = A'r.
   real(wp),intent(out) :: w(n)     !! On exit, w = A'r - damp^2 x.
   real(wp),intent(in)  :: x(n)     !! The given estimate of a solution.

   real(wp),parameter :: power = 0.5_wp

   integer :: j
   real(wp) :: bnorm, dampsq, rho1, rho2, sigma1, sigma2, &
               tol, snorm, xnorm, xsnorm
   real(wp),dimension(n) :: xtmp

   dampsq = damp**2
   tol    = eps**power
   xtmp   = x

   ! Compute u = b - Ax via u = -b + Ax, u = -u.
   ! This is usual residual vector r.

   call dcopy ( m, b, 1, u, 1 )
   call dscal ( m, (-one), u, 1 )
   call me%aprod ( 1, m, n, xtmp, u )
   call dscal ( m, (-one), u, 1 )

   ! Compute v = A'u via v = 0, v = v + A'u.

   do j = 1, n
      v(j)  = zero
   end do
   call me%aprod ( 2, m, n, v, u )

   ! Compute w = A'u - damp**2 * x.
   ! This will be close to zero in all cases
   ! if x is close to a solution.

   call dcopy ( n, v, 1, w, 1 )
   if (damp /= zero) then
      do j = 1, n
         w(j)  = w(j)  -  dampsq * x(j)
      end do
   end if

   ! Compute the norms associated with b, x, u, v, w.

   bnorm  = dnrm2 ( m, b, 1 )
   xnorm  = dnrm2 ( n, x, 1 )
   rho1   = dnrm2 ( m, u, 1 )
   sigma1 = dnrm2 ( n, v, 1 )
   if (nout /= 0) then
      write(nout, '(//A)') 'Enter xcheck. Does x solve Ax = b, etc?'
      write(nout, '(1P,A,E10.3)')      ' damp            =', damp
      write(nout, '(1P,A,E10.3)')      ' norm(x)         =', xnorm
      write(nout, '(1P,A,E15.8,A)')    ' norm(r)         =', rho1,   ' = rho1'
      write(nout, '(1P,A,E10.3,5X,A)') ' norm(A''r)       =', sigma1, ' = sigma1'
   end if

   if (damp == zero) then
      rho2   = rho1
      sigma2 = sigma1
   else
      rho2   = sqrt( rho1**2  +  dampsq * xnorm**2 )
      sigma2 = dnrm2 ( n, w, 1 )
      snorm  = rho1 / damp
      xsnorm = rho2 / damp
      if (nout /= 0) then
         write(nout, '(1P/A,E10.3)')      ' norm(s)         =', snorm
         write(nout, '(1P,A,E10.3)')      ' norm(x,s)       =', xsnorm
         write(nout, '(1P,A,E15.8,A)')    ' norm(rbar)      =', rho2,    ' = rho2'
         write(nout, '(1P,A,E10.3,5X,A)') ' norm(Abar''rbar) =', sigma2, ' = sigma2'
      end if
   end if

   ! See if x seems to solve Ax = b or min norm(Ax - b)
   ! or the damped least-squares system.

   if (bnorm == zero  .and.  xnorm == zero) then
      inform = 0
      test1  = zero
      test2  = zero
      test3  = zero
   else
      inform = 4
      test1  = rho1 / (bnorm  +  anorm * xnorm)
      test2  = zero
      if (rho1 > zero) test2  = sigma1 / (anorm * rho1)
      test3  = test2
      if (rho2 > zero) test3  = sigma2 / (anorm * rho2)
      if (test3 <= tol) inform = 3
      if (test2 <= tol) inform = 2
      if (test1 <= tol) inform = 1
   end if

   if (nout /= 0) then
      write(nout, '(/A,I2)')        ' inform          =', inform
      write(nout, '(1P,A,E10.3)')   ' tol             =', tol
      write(nout, '(1P,A,E10.3,A)') ' test1           =', test1, ' (Ax = b)'
      write(nout, '(1P,A,E10.3,A)') ' test2           =', test2, ' (least-squares)'
      write(nout, '(1P,A,E10.3,A)') ' test3           =', test3, ' (damped least-squares)'
   end if

   end subroutine xcheck
!***************************************************************************************************

!***************************************************************************************************
!>
!  Returns \( \sqrt{ a^2 + b^2 } \) with precautions to avoid overflow.
!
!### History
! * 21 Mar 1990: First version.

   pure function d2norm( a, b )

   real(wp) :: d2norm
   real(wp),intent(in) :: a
   real(wp),intent(in) :: b

   real(wp) :: scale

   scale  = abs( a ) + abs( b )
   if (scale == zero) then
      d2norm = zero
   else
      d2norm = scale * sqrt( (a/scale)**2   +  (b/scale)**2 )
   end if

   end function d2norm
!***************************************************************************************************

!***************************************************************************************************
    end module lsqr_module
!***************************************************************************************************

max_stars_repo_name: jacobwilliams/LSQR

max_issues_repo_name: jacobwilliams/LSQR

max_forks_repo_name: jacobwilliams/LSQR