share/colnew/fortran/contrl.f

   1       SUBROUTINE CONTRL (XI, XIOLD, Z, DMZ, RHS, DELZ, DELDMZ,
   2      1           DQZ, DQDMZ, G, W, V, VALSTR, SLOPE, SCALE, DSCALE,
   3      2           ACCUM, IPVTG, INTEGS, IPVTW, NFXPNT, FIXPNT, IFLAG,
   4      3           FSUB, DFSUB, GSUB, DGSUB, GUESS )
   5 C
   6 C**********************************************************************
   7 C
   8 C   purpose
   9 C     this subroutine is the actual driver.  the nonlinear iteration
  10 C     strategy is controlled here ( see [4] ). upon convergence, errchk
  11 C     is called to test for satisfaction of the requested tolerances.
  12 C
  13 C   variables
  14 C
  15 C     check  - maximum tolerance value, used as part of criteria for
  16 C              checking for nonlinear iteration convergence
  17 C     relax  - the relaxation factor for damped newton iteration
  18 C     relmin - minimum allowable value for relax  (otherwise the
  19 C              jacobian is considered singular).
  20 C     rlxold - previous relax
  21 C     rstart - initial value for relax when problem is sensitive
  22 C     ifrz   - number of fixed jacobian iterations
  23 C     lmtfrz - maximum value for ifrz before performing a reinversion
  24 C     iter   - number of iterations (counted only when jacobian
  25 C              reinversions are performed).
  26 C     xi     - current mesh
  27 C     xiold  - previous mesh
  28 C     ipred  = 0  if relax is determined by a correction
  29 C            = 1  if relax is determined by a prediction
  30 C     ifreez = 0  if the jacobian is to be updated
  31 C            = 1  if the jacobian is currently fixed (frozen)
  32 C     iconv  = 0  if no previous convergence has been obtained
  33 C            = 1  if convergence on a previous mesh has been obtained
  34 C     icare  =-1  no convergence occurred (used for regular problems)
  35 C            = 0  a regular problem
  36 C            = 1  a sensitive problem
  37 C            = 2  used for continuation (see description of ipar(10)
  38 C                 in colnew).
  39 C     rnorm  - norm of rhs (right hand side) for current iteration
  40 C     rnold  - norm of rhs for previous iteration
  41 C     anscl  - scaled norm of newton correction
  42 C     anfix  - scaled norm of newton correction at next step
  43 C     anorm  - scaled norm of a correction obtained with jacobian fixed
  44 C     nz     - number of components of  z  (see subroutine approx)
  45 C     ndmz   - number of components of  dmz  (see subroutine approx)
  46 C     imesh  - a control variable for subroutines newmsh and errchk
  47 C            = 1  the current mesh resulted from mesh selection
  48 C                 or is the initial mesh.
  49 C            = 2  the current mesh resulted from doubling the
  50 C                 previous mesh
  51 C
  52 C**********************************************************************
  53 C
  54       IMPLICIT REAL*8 (A-H,O-Z)
  55       DIMENSION XI(1), XIOLD(1), Z(1), DMZ(1), RHS(1)
  56       DIMENSION G(1), W(1), V(1), VALSTR(1), SLOPE(1), ACCUM(1)
  57       DIMENSION DELZ(1), DELDMZ(1), DQZ(1), DQDMZ(1) , FIXPNT(1)
  58       DIMENSION DUMMY(1), SCALE(1), DSCALE(1)
  59       DIMENSION INTEGS(1), IPVTG(1), IPVTW(1)
  60 C
  61       COMMON /COLOUT/ PRECIS, IOUT, IPRINT
  62       COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
  63       COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
  64       COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
  65       COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
  66       COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
  67       COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
  68      1                ROOT(40), JTOL(40), LTOL(40), NTOL
  69 C
  70       EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS
  71 C
  72 C...  constants for control of nonlinear iteration
  73 C
  74       RELMIN = 1.D-3
  75       RSTART = 1.D-2
  76       LMTFRZ = 4
  77 C
  78 C...  compute the maximum tolerance
  79 C
  80       CHECK = 0.D0
  81       DO 10 I = 1, NTOL
  82    10   CHECK = DMAX1 ( TOLIN(I), CHECK )
  83       IMESH = 1
  84       ICONV = 0
  85       IF ( NONLIN .EQ. 0 ) ICONV = 1
  86       ICOR = 0
  87       NOCONV = 0
  88       MSING = 0
  89 C
  90 C...  the main iteration begins here .
  91 C...  loop 20 is executed until error tolerances are satisfied or
  92 C...  the code fails (due to a singular matrix or storage limitations)
  93 C
  94    20      CONTINUE
  95 C
  96 C...       initialization for a new mesh
  97 C
  98            ITER = 0
  99            IF ( NONLIN .GT. 0 )                     GO TO 50
 100 C
 101 C...       the linear case.
 102 C...       set up and solve equations
 103 C
 104            CALL LSYSLV (MSING, XI, XIOLD, DUMMY, DUMMY, Z, DMZ, G,
 105      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 0,
 106      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 107 C
 108 C...       check for a singular matrix
 109 C
 110            IF ( MSING .EQ. 0 )                      GO TO 400
 111    30      IF ( MSING .LT. 0 )                      GO TO 40
 112            IF ( IPRINT .LT. 1 )  WRITE (IOUT,495)
 113            GO TO 460
 114    40      IF ( IPRINT .LT. 1 )  WRITE (IOUT,490)
 115            IFLAG = 0
 116            RETURN
 117 C
 118 C...       iteration loop for nonlinear case
 119 C...       define the initial relaxation parameter (= relax)
 120 C
 121    50      RELAX = 1.D0
 122 C
 123 C...       check for previous convergence and problem sensitivity
 124 C
 125            IF ( ICARE .EQ. 1 .OR. ICARE .EQ. (-1) )  RELAX = RSTART
 126            IF ( ICONV .EQ. 0 )                      GO TO 160
 127 C
 128 C...       convergence on a previous mesh has been obtained.    thus
 129 C...       we have a very good initial approximation for the newton
 130 C...       process.    proceed with one full newton and then iterate
 131 C...       with a fixed jacobian.
 132 C
 133            IFREEZ = 0
 134 C
 135 C...       evaluate right hand side and its norm  and
 136 C...       find the first newton correction
 137 C
 138            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
 139      1          W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1,
 140      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 141 C
 142            IF ( IPRINT .LT. 0 )  WRITE(IOUT,530)
 143            IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNOLD
 144            GO TO 70
 145 C
 146 C...       solve for the next iterate .
 147 C...       the value of ifreez determines whether this is a full
 148 C...       newton step (=0) or a fixed jacobian iteration (=1).
 149 C
 150    60      IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNORM
 151            RNOLD = RNORM
 152            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
 153      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM,
 154      2          3+IFREEZ, FSUB, DFSUB, GSUB, DGSUB, GUESS )
 155 C
 156 C...       check for a singular matrix
 157 C
 158    70      IF ( MSING .NE. 0 )                      GO TO 30
 159            IF ( IFREEZ .EQ. 1 )                     GO TO 80
 160 C
 161 C...       a full newton step
 162 C
 163            ITER = ITER + 1
 164            IFRZ = 0
 165    80      CONTINUE
 166 C
 167 C...       update   z and dmz , compute new  rhs  and its norm
 168 C
 169            DO 90 I = 1, NZ
 170              Z(I) = Z(I) + DELZ(I)
 171    90      CONTINUE
 172            DO 100 I = 1, NDMZ
 173              DMZ(I) = DMZ(I) + DELDMZ(I)
 174   100      CONTINUE
 175            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
 176      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,
 177      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 178 C
 179 C...       check monotonicity. if the norm of  rhs  gets smaller,
 180 C...       proceed with a fixed jacobian; else proceed cautiously,
 181 C...       as if convergence has not been obtained before (iconv=0).
 182 C
 183            IF ( RNORM .LT. PRECIS )                 GO TO 390
 184            IF ( RNORM .GT. RNOLD )                  GO TO 130
 185            IF ( IFREEZ .EQ. 1 )                     GO TO 110
 186            IFREEZ = 1
 187            GO TO 60
 188 C
 189 C...       verify that the linear convergence with fixed jacobian
 190 C...       is fast enough.
 191 C
 192   110      IFRZ = IFRZ + 1
 193            IF ( IFRZ .GE. LMTFRZ )       IFREEZ = 0
 194            IF ( RNOLD .LT. 4.D0*RNORM )  IFREEZ = 0
 195 C
 196 C...       check convergence (iconv = 1).
 197 C
 198            DO 120 IT = 1, NTOL
 199              INZ = LTOL(IT)
 200              DO 120 IZ = INZ, NZ, MSTAR
 201                IF ( DABS(DELZ(IZ)) .GT.
 202      1           TOLIN(IT) * (DABS(Z(IZ)) + 1.D0))  GO TO 60
 203   120      CONTINUE
 204 C
 205 C...       convergence obtained
 206 C
 207            IF ( IPRINT .LT. 1 )  WRITE (IOUT,560) ITER
 208            GO TO 400
 209 C
 210 C...      convergence of fixed jacobian iteration failed.
 211 C
 212   130      IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNORM
 213            IF ( IPRINT .LT. 0 )  WRITE (IOUT,540)
 214            ICONV = 0
 215            RELAX = RSTART
 216            DO 140 I = 1, NZ
 217              Z(I) = Z(I) - DELZ(I)
 218   140      CONTINUE
 219            DO 150 I = 1, NDMZ
 220              DMZ(I) = DMZ(I) - DELDMZ(I)
 221   150      CONTINUE
 222 C
 223 C...       update old mesh
 224 C
 225            NP1 = N + 1
 226            DO 155 I = 1, NP1
 227   155        XIOLD(I) = XI(I)
 228            NOLD = N
 229 C
 230            ITER = 0
 231 C
 232 C...       no previous convergence has been obtained. proceed
 233 C...       with the damped newton method.
 234 C...       evaluate rhs and find the first newton correction.
 235 C
 236   160      IF(IPRINT .LT. 0)  WRITE (IOUT,500)
 237            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
 238      1          W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1,
 239      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 240 C
 241 C...       check for a singular matrix
 242 C
 243            IF ( MSING .NE. 0 )                       GO TO 30
 244 C
 245 C...       bookkeeping for first mesh
 246 C
 247            IF ( IGUESS .EQ. 1 )  IGUESS = 0
 248 C
 249 C...       find initial scaling
 250 C
 251            CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
 252            GO TO 220
 253 C
 254 C...       main iteration loop
 255 C
 256   170      RNOLD = RNORM
 257            IF ( ITER .GE. LIMIT )                   GO TO 430
 258 C
 259 C...       update scaling
 260 C
 261            CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
 262 C
 263 C...       compute norm of newton correction with new scaling
 264 C
 265            ANSCL = 0.D0
 266            DO 180 I = 1, NZ
 267              ANSCL = ANSCL + (DELZ(I) * SCALE(I))**2
 268   180      CONTINUE
 269            DO 190 I = 1, NDMZ
 270              ANSCL = ANSCL + (DELDMZ(I) * DSCALE(I))**2
 271   190      CONTINUE
 272            ANSCL = DSQRT(ANSCL / DFLOAT(NZ+NDMZ))
 273 C
 274 C...       find a newton direction
 275 C
 276            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
 277      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 3,
 278      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 279 C
 280 C...       check for a singular matrix
 281 C
 282            IF ( MSING .NE. 0 )                      GO TO 30
 283 C
 284 C...       predict relaxation factor for newton step.
 285 C
 286            ANDIF = 0.D0
 287            DO 200 I = 1, NZ
 288              ANDIF = ANDIF + ((DQZ(I) - DELZ(I)) * SCALE(I))**2
 289   200      CONTINUE
 290            DO 210 I = 1, NDMZ
 291              ANDIF = ANDIF + ((DQDMZ(I) - DELDMZ(I)) * DSCALE(I))**2
 292   210      CONTINUE
 293            ANDIF = DSQRT(ANDIF/DFLOAT(NZ+NDMZ) + PRECIS)
 294            RELAX = RELAX * ANSCL / ANDIF
 295            IF ( RELAX .GT. 1.D0 )  RELAX = 1.D0
 296   220      RLXOLD = RELAX
 297            IPRED = 1
 298            ITER = ITER + 1
 299 C
 300 C...       determine a new  z and dmz  and find new  rhs  and its norm
 301 C
 302            DO 230 I = 1, NZ
 303              Z(I) = Z(I)  +  RELAX * DELZ(I)
 304   230      CONTINUE
 305            DO 240 I = 1, NDMZ
 306              DMZ(I) = DMZ(I)  +  RELAX * DELDMZ(I)
 307   240      CONTINUE
 308   250      CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DQZ, DQDMZ, G,
 309      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,
 310      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 311 C
 312 C...       compute a fixed jacobian iterate (used to control relax)
 313 C
 314            CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DQZ, DQDMZ, G,
 315      1          W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 4,
 316      2          FSUB, DFSUB, GSUB, DGSUB, GUESS )
 317 C
 318 C...       find scaled norms of various terms used to correct relax
 319 C
 320            ANORM = 0.D0
 321            ANFIX = 0.D0
 322            DO 260 I = 1, NZ
 323              ANORM = ANORM  +  (DELZ(I) * SCALE(I))**2
 324              ANFIX = ANFIX  +  (DQZ(I) * SCALE(I))**2
 325   260      CONTINUE
 326            DO 270 I = 1, NDMZ
 327              ANORM = ANORM  +  (DELDMZ(I) * DSCALE(I))**2
 328              ANFIX = ANFIX  +  (DQDMZ(I) * DSCALE(I))**2
 329   270      CONTINUE
 330            ANORM = DSQRT(ANORM / DFLOAT(NZ+NDMZ))
 331            ANFIX = DSQRT(ANFIX / DFLOAT(NZ+NDMZ))
 332            IF ( ICOR .EQ. 1 )                         GO TO 280
 333            IF (IPRINT .LT. 0)  WRITE (IOUT,520) ITER, RELAX, ANORM,
 334      1           ANFIX, RNOLD, RNORM
 335            GO TO 290
 336   280      IF (IPRINT .LT. 0) WRITE (IOUT,550) RELAX, ANORM, ANFIX,
 337      1           RNOLD, RNORM
 338   290      ICOR = 0
 339 C
 340 C...       check for monotonic decrease in  delz and deldmz.
 341 C
 342            IF (ANFIX.LT.PRECIS .OR. RNORM.LT.PRECIS)  GO TO 390
 343            IF ( ANFIX .GT. ANORM )                    GO TO 300
 344 C
 345 C...       we have a decrease.
 346 C...       if  dqz  and dqdmz  small, check for convergence
 347 C
 348            IF ( ANFIX .LE. CHECK )                    GO TO 350
 349 C
 350 C...       correct the predicted  relax  unless the corrected
 351 C...       value is within 10 percent of the predicted one.
 352 C
 353            IF ( IPRED .NE. 1 )                        GO TO 170
 354   300      IF ( ITER .GE. LIMIT )                     GO TO 430
 355 C
 356 C...       correct the relaxation factor.
 357 C
 358            IPRED = 0
 359            ARG = (ANFIX/ANORM - 1.D0) / RELAX + 1.D0
 360            IF ( ARG .LT. 0.D0 )                       GO TO 170
 361            IF (ARG .LE. .25D0*RELAX+.125D0*RELAX**2)  GO TO 310
 362            FACTOR = -1.D0 + DSQRT (1.D0+8.D0 * ARG)
 363            IF ( DABS(FACTOR-1.D0) .LT. .1D0*FACTOR )  GO TO 170
 364            IF ( FACTOR .LT. 0.5D0 )  FACTOR = 0.5D0
 365            RELAX = RELAX / FACTOR
 366            GO TO 320
 367   310      IF ( RELAX .GE. .9D0 )                     GO TO 170
 368            RELAX = 1.D0
 369   320      ICOR = 1
 370            IF ( RELAX .LT. RELMIN )                   GO TO 440
 371            FACT = RELAX - RLXOLD
 372            DO 330 I = 1, NZ
 373             Z(I) = Z(I)  +  FACT * DELZ(I)
 374   330      CONTINUE
 375            DO 340 I = 1, NDMZ
 376              DMZ(I) = DMZ(I)  +  FACT * DELDMZ(I)
 377   340      CONTINUE
 378            RLXOLD = RELAX
 379            GO TO 250
 380 C
 381 C...       check convergence (iconv = 0).
 382 C
 383   350      CONTINUE
 384            DO 360 IT = 1, NTOL
 385              INZ = LTOL(IT)
 386              DO 360 IZ = INZ, NZ, MSTAR
 387                IF ( DABS(DQZ(IZ)) .GT.
 388      1         TOLIN(IT) * (DABS(Z(IZ)) + 1.D0) )   GO TO 170
 389   360      CONTINUE
 390 C
 391 C...       convergence obtained
 392 C
 393            IF ( IPRINT .LT. 1 )  WRITE (IOUT,560) ITER
 394 C
 395 C...       since convergence obtained, update  z and dmz  with term
 396 C...       from the fixed jacobian iteration.
 397 C
 398            DO 370 I = 1, NZ
 399              Z(I) = Z(I)  +  DQZ(I)
 400   370      CONTINUE
 401            DO 380 I = 1, NDMZ
 402              DMZ(I) = DMZ(I)  +  DQDMZ(I)
 403   380      CONTINUE
 404   390      IF ( (ANFIX .LT. PRECIS .OR. RNORM .LT. PRECIS)
 405      1          .AND. IPRINT .LT. 1 )  WRITE (IOUT,560) ITER
 406            ICONV = 1
 407            IF ( ICARE .EQ. (-1) )  ICARE = 0
 408 C
 409 C...       if full output has been requested, print values of the
 410 C...       solution components   z  at the meshpoints.
 411 C
 412   400      IF ( IPRINT .GE. 0 )                     GO TO 420
 413            DO 410 J = 1, MSTAR
 414              WRITE(IOUT,610) J
 415   410      WRITE(IOUT,620) (Z(LJ), LJ = J, NZ, MSTAR)
 416 C
 417 C...       check for error tolerance satisfaction
 418 C
 419   420      IFIN = 1
 420            IF (IMESH .EQ. 2) CALL ERRCHK (XI, Z, DMZ, VALSTR, IFIN)
 421            IF ( IMESH .EQ. 1 .OR.
 422      1          IFIN .EQ. 0 .AND. ICARE .NE. 2)     GO TO 460
 423            IFLAG = 1
 424            RETURN
 425 C
 426 C...       diagnostics for failure of nonlinear iteration.
 427 C
 428   430      IF ( IPRINT .LT. 1 )  WRITE (IOUT,570) ITER
 429            GO TO 450
 430   440      IF( IPRINT .LT. 1 )  WRITE(IOUT,580) RELAX, RELMIN
 431   450      IFLAG = -2
 432            NOCONV = NOCONV + 1
 433            IF ( ICARE .EQ. 2 .AND. NOCONV .GT. 1 )  RETURN
 434            IF ( ICARE .EQ. 0 )  ICARE = -1
 435 C
 436 C...       update old mesh
 437 C
 438   460      NP1 = N + 1
 439            DO 470 I = 1, NP1
 440   470        XIOLD(I) = XI(I)
 441            NOLD = N
 442 C
 443 C...       pick a new mesh
 444 C...       check safeguards for mesh refinement
 445 C
 446            IMESH = 1
 447            IF ( ICONV .EQ. 0 .OR. MSHNUM .GE. MSHLMT
 448      1                       .OR. MSHALT .GE. MSHLMT )  IMESH = 2
 449            IF ( MSHALT .GE. MSHLMT .AND.
 450      1          MSHNUM .LT. MSHLMT )  MSHALT = 1
 451            CALL NEWMSH (IMESH, XI, XIOLD, Z, DMZ, VALSTR,
 452      1                  SLOPE, ACCUM, NFXPNT, FIXPNT)
 453 C
 454 C...       exit if expected n is too large (but may try n=nmax once)
 455 C
 456            IF ( N .LE. NMAX )                       GO TO 480
 457            N = N / 2
 458            IFLAG = -1
 459            IF ( ICONV .EQ. 0 .AND. IPRINT .LT. 1 )  WRITE (IOUT,590)
 460            IF ( ICONV .EQ. 1 .AND. IPRINT .LT. 1 )  WRITE (IOUT,600)
 461            RETURN
 462   480      IF ( ICONV .EQ. 0 )  IMESH = 1
 463            IF ( ICARE .EQ. 1 )  ICONV = 0
 464       GO TO 20
 465 C     ---------------------------------------------------------------
 466   490 FORMAT(//35H THE GLOBAL BVP-MATRIX IS SINGULAR )
 467   495 FORMAT(//40H A LOCAL ELIMINATION MATRIX IS SINGULAR )
 468   500 FORMAT(/30H FULL DAMPED NEWTON ITERATION,)
 469   510 FORMAT(13H ITERATION = , I3, 15H  NORM (RHS) = , D10.2)
 470   520 FORMAT(13H ITERATION = ,I3,22H  RELAXATION FACTOR = ,D10.2
 471      1       /33H NORM OF SCALED RHS CHANGES FROM ,D10.2,3H TO,D10.2
 472      2       /33H NORM   OF   RHS  CHANGES  FROM  ,D10.2,3H TO,D10.2,
 473      2       D10.2)
 474   530 FORMAT(/27H FIXED JACOBIAN ITERATIONS,)
 475   540 FORMAT(/35H SWITCH TO DAMPED NEWTON ITERATION,)
 476   550 FORMAT(40H RELAXATION FACTOR CORRECTED TO RELAX = , D10.2
 477      1       /33H NORM OF SCALED RHS CHANGES FROM ,D10.2,3H TO,D10.2
 478      2       /33H NORM   OF   RHS  CHANGES  FROM  ,D10.2,3H TO,D10.2
 479      2       ,D10.2)
 480   560 FORMAT(/18H CONVERGENCE AFTER , I3,11H ITERATIONS /)
 481   570 FORMAT(/22H NO CONVERGENCE AFTER , I3, 11H ITERATIONS/)
 482   580 FORMAT(/37H NO CONVERGENCE.  RELAXATION FACTOR =,D10.3
 483      1       ,24H IS TOO SMALL (LESS THAN, D10.3, 1H)/)
 484   590 FORMAT(18H  (NO CONVERGENCE) )
 485   600 FORMAT(50H  (PROBABLY TOLERANCES TOO STRINGENT, OR NMAX TOO
 486      1       ,6HSMALL) )
 487   610 FORMAT( 19H MESH VALUES FOR Z(, I2, 2H), )
 488   620 FORMAT(1H , 8D15.7)
 489       END