text:
      SUBROUTINE STGPD (X, Y, Z, DZ, NP, TRIANG, NT, IWORK)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c  File: STGPD.[f|for] contains STGPD, sTGMOR, sTGLS.
C>> 2007-02-28 STGPD Krogh Massive changed to remove assigned go to's.
c>> 2005-12-07 STGPD Krogh  Removed unused label.
c>> 1996-05-10 STGPD Krogh  MW removed from dim. declarator in C version
c>> 1996-03-30 STGPD Krogh  MIN0 => MIN
C>> 1996-02-02 STGPD CLL
C>> 1996-01-15 STGPD CLL
C>> 1996-01-11 STGPD CLL
C>> 1995-09-26 STGPD CLL Editing for inclusion into MATH77.
C
c  THIS SUBR ESTIMATES FIRST PARTIAL DERIVS AT THE GIVEN DATA POINTS.
c  C.L.LAWSON, JPL, 1976 DEC 21.  EDITED COMMENTS 1979 MAR 5.
C
c  ---------------------------------------------------------------------
C                     Method
c
c      THE PARTIAL DERIVS AT POINT  JP  WILL BE ESTIMATED BY FITTING
c  A QUADRATIC IN X AND Y TO A SET OF  MFIT+1  POINTS CONSISTING OF
c  POINT  JP  AND   MFIT NEARBY POINTS IN THE TRIANGULAR GRID.
c         THE VALUE OF MFIT WILL IN GENERAL BE DIFFERENT AT DIFFERENT
c  POINTS.  MFIT WILL BE IN THE RANGE FROM MFMIN TO MFMAX.  LET NNEB
c  DENOTE THE NO. OF IMMEDIATE NEIGHBORING POINTS TO POINT JB.  IF
c  MFMIN .LE. NNEB .LE. MFMAX  THEN MFIT=NNEB.
c      IF NNEB .LT. MFMIN THEN MFIT = MFMIN.
c      IF NNEB .GT. MFMAX THEN MFIT = MFMAX.
c  ---------------------------------------------------------------------
c                 Subroutine Arguments
c
c  X(), Y(), Z()     [in]
c  DZ(2,)            [out]
c  NP                [in]
c  TRIANG(), NT  [in]
c  (IWORK(I),I=1,NP) [scratch]  INITIALLY SET TO ZERO BY STGPD.
c       STGPD WILL SET IWORK(JP)= 1 WHEN IT ESTIMATES THE PARTIAL
c       DERIVATIVES AT POINT JP.
C
c  ---------------------------------------------------------------------
C                     COMMON BLOCK
c
c  /sTGCM1/ provides communication between STGPD and sTGMOR.
c  ---------------------------------------------------------------------
c                 Internal variables
C
c  AVAIL  POINTS TO AVAILABLE SPACE IN KE().
c  MFMAX  Dimensioning parameter used in the following declarations:
c           INTEGER   KE(MFMAX  ,4), JUSED(MFMAX)
c           REAL             W(MFMAX+5,6), DSQ(MFMAX)
c  MFMIN  LOWER BOUND FOR MFIT.  SET TO A CONSTANT VALUE BY STGPD.
c  MFIT   NO. OF NEIGHBORING POINTS USED WITH EACH POINT JP
C         TO ESTIMATE PARTIALS AT POINT JP.  GENERALLY DIFFERS
C         FOR EACH POINT JP.
c  FIRSTL,LASTL     POINTERS TO FIRST AND LAST ELEMENTS OF THE
c                   LIST CONTAINED IN KE().
c  KE(,) IS A DOUBLY  LINKED LIST DESCRIBING A CYCLE OF EDGES
c       SURROUNDING POINT JP.  AN EDGE IS DEFINED BY ITS TWO ADJACENT
c       TRIANGLES.
c           KE(L,1)   FWD POINTER  (COUNTERCLOCKWISE).
c           KE(L,2)   BKWD POINTER (CLOCKWISE).
c           KE(L,3)   TRIANGLE INSIDE EDGE L.
c           KE(L,4)   TRIANGLE OUTSIDE EDGE L,  OR ELSE ZERO IF EDGE
c                     ON THE BOUNDARY.
c  JUSED()   AN UNORDERED SET OF INDICES OF ROINTS USED IN THE LOCAL
c            FIT THAT ARE NOT IMMEDIATE NEIGHBORS OF POINT JP.
c  X0,Y0,Z0  COORDINATES OF POINT JP.
c  DSQ(L)    SQUARED DISTANCE FROM POINT JP TO THE MIDPOINT OF SIDE L
c  NEWPT     LOGICAL FLAG.  FOR A NEW POINT JP STGPD SETS
c            NEWPT TRUE.  THIS CAUSES sTGLS TO TRIANGULARIZE THE
c            LEAST SQUARES SYSTEM FOR JP. sTGLS THEN TESTS THE
c            CONDITION OF THE SYSTEM. IF IT IS BAD sTGLS SETS NEWPT
c            FALSE AND RETURNS. ESTPT THEN USES sTGMOR TO GET ONE
c            MORE POINT AND CALLS sTGLS AGAIN, LEAVING NEWPT
c            FALSE.  sTGLS THEN ACCUMULATES THE ADDITIONAL POINT
c            AND AGAIN TESTS THE CONDITION.  EVENTUALLY EITHER
c            ENOUGH POINTS ARE ADDED TO OBTAIN ADEQUATE CONDITION
c            OR ELSE sTGLS APPLIES AN ARBITRARY STABILIZATION AND
c            SO IN EITHER CASE sTGLS FINALLY RETURNS A SOLUTION.
C
c     ------------------------------------------------------------------
c--S replaces "?": ?TGPD, ?TGGET, ?TGMOR, ?TGLS, ?TGCM1
c--&               ?TGMOR, ?TGCM1, ?TGLS,  ?ROT,  ?ROTG
c     ------------------------------------------------------------------
      common  /sTGCM1/ X0, Y0, Z0, DSQ,
     *  AVAIL, MFMIN, MFIT, FIRSTL, LASTL, KE, JUSED, NEWPT
c Note -- if nfmax is changed, mw changes and C code must change
c everywhere mw is used in a dimension declarator.
      integer mfmax
      parameter(mfmax = 16)
      real             X0, Y0, Z0, DSQ(mfmax)
      integer AVAIL,MFMIN,MFIT, FIRSTL, LASTL, KE(mfmax,4), JUSED(mfmax)
      logical NEWPT
      integer mw
      parameter(mw = mfmax+5)
      integer NP
      integer ADD1(3), i, IWORK(np), ip1, ipcent, it, iv1, iv2, ivert
      integer j1, j2, jj1, jj2, jp
      integer kp, l, LIMIT
      integer lolim
      integer NT
      integer SUB1(3)
      integer TNOW, TRI(6), TRIANG(1:*), TRSAVE(6)
c++ Code for .C. is INACTIVE
C      real             dsfcn, DZ(2,NP), W(21,6)
c++ Code for ~.C. is ACTIVE
      real             dsfcn, DZ(2,NP), W(mw,6)
c++ END
      real             X(NP), x1, x2, Y(NP), y1, y2, Z(NP)
      logical HITBDY
C
      data ADD1(1),ADD1(2),ADD1(3),SUB1(1),SUB1(2),SUB1(3)/2,3,1, 3,1,2/
c     ------------------------------------------------------------------
      DSFCN(X1,Y1,X2,Y2) = ((X1+X2)*.5 - X0)**2 + ((Y1+Y2)*.5 - Y0)**2
c     ------------------------------------------------------------------
      MFMIN = 6
      LIMIT = min( MFMAX, NP - 1 )
      LOLIM = min( MFMIN, NP - 1 )
      NEWPT = .true.
C                        SET IWORK() = 0
      do 20 I = 1 , NP
        IWORK(I) = 0
 20   continue
C                                  MAIN LOOP THRU TRIANGLES
C
      IT = 1
      go to 40
 30   IT = IT + 1
 40   if (NT .lt. IT) return
      call STGGET(IT , TRSAVE, TRIANG)
C                                       LOOP THRU VERTICES OF A TRIANGLE
      IVERT = 1
      go to 60
 50   IVERT = IVERT + 1
 60   if (IVERT .gt. 3) go to 30
      JP= TRSAVE(IVERT + 3)
      if (IWORK(JP) .ne. 0) go to 50
      IWORK(JP) = 1
C1001 FORMAT(39H0ESTPD..  ESTIMATING PARTIALS AT POINT ,I5/1X)
C     THE FORTRAN END FOLLOWS..
C              WRITE(*,1001) JP
C           PARTIALS FOR POINT JP
      X0 = X(JP)
      Y0 = Y(JP)
      Z0 = Z(JP)
C                                  SET POINTERS FOR AVAILABLE SPACE.
      do 80 I = 2, MFMAX
        KE(I - 1,1) = I
 80   continue
      AVAIL = 1
      KE(MFMAX,1) = 0
C                                            BUILD FITTING EQUATIONS.
      go to 120
 90   if (MFIT .lt. LOLIM) then
        call sTGMOR ( X,Y,Z,NP, TRIANG,W)
      end if
C                                            SOLVE FITTING EQUATIONS.
C
      call sTGLS(W,MW,MFIT,NEWPT, LIMIT, DZ(1,JP),DZ(2,JP))
C                     GO FOR MORE POINTS IF NEEDED TO IMPROVE CONDITION.
 100  if (.not. NEWPT) then
        call sTGMOR (X,Y,Z,NP, TRIANG,W)
        call sTGLS(W,MW,MFIT,NEWPT, LIMIT, DZ(1,JP),DZ(2,JP))
        go to 100
      end if
C
      go to 50
C          PROCESS RING OF NEIGHBORS
 120  do 130 I = 1,6
        TRI(I) = TRSAVE(I)
 130  continue
C            POP AVAIL TO L
      L = AVAIL
      AVAIL = KE(AVAIL,1)
C                             BUILD FIRST EDGE IN RING OF NEIGHBORS.
      FIRSTL = L
      LASTL =L
      KE(L,1) = 0
      KE(L,2) = 0
      KE(L,3) = IT
      IV1=   ADD1(IVERT)
      KE(L,4)= TRI(IV1)
      J1 = TRI(IV1 + 3)
      IV2=   SUB1(IVERT)
      J2 = TRI(IV2 + 3)
      DSQ(L) = DSFCN(X(J1),Y(J1),X(J2),Y(J2))
C
C                            BUILD EQUAS FOR FIRST TWO NEIGHBORING PTS.
      MFIT = 1
      KP = J1
C                   BUILD ROW
C1000 FORMAT(25H BUILD ROW..  CENTER PT =,I5,12H,  EQUA NO =,I3,
C    *13H  NEARBY PT =,I5)
C        WRITE(*,1000) JP,MFIT,KP
      JUSED(MFIT) = KP
      W(MFIT,4) = X(KP) - X0
      W(MFIT,5) = Y(KP) - Y0
      W(MFIT,6) = Z(KP) - Z0
      W(MFIT,1) = W(MFIT,4)**2
      W(MFIT,2) = W(MFIT,4) * W(MFIT,5)
      W(MFIT,3) = W(MFIT,5)**2
      MFIT = 2
      KP = J2
C                   BUILD ROW
C1000 FORMAT(25H BUILD ROW..  CENTER PT =,I5,12H,  EQUA NO =,I3,
C    *13H  NEARBY PT =,I5)
C        WRITE(*,1000) JP,MFIT,KP
      JUSED(MFIT) = KP
      W(MFIT,4) = X(KP) - X0
      W(MFIT,5) = Y(KP) - Y0
      W(MFIT,6) = Z(KP) - Z0
      W(MFIT,1) = W(MFIT,4)**2
      W(MFIT,2) = W(MFIT,4) * W(MFIT,5)
      W(MFIT,3) = W(MFIT,5)**2
      go to  180
C
C                             MOVE LEFT AND THEN RIGHT AROUND POINT JP
C                             TO BUILD EQUAS FOR IMMEDIATELY
C                             NEIGHBORING POINTS.
C
 160  if (.not.  HITBDY) go to 90
      do 170 I = 1,6
        TRI(I) = TRSAVE(I)
 170  continue
      go to 360
C                    MOVE LEFT
C                                     LEFT MEANS COUNTERCLOCKWISE.
 180  JJ2 =  J2
C                                             *
C      ****************************************
      IPCENT = IVERT
      HITBDY = .false.
      if (MFIT .ge. MFMAX) go to 160
      I     = SUB1(IPCENT)
      TNOW = TRI(I)
      if (TNOW .ne. 0) go to 200
C                                            HAVE HIT BOUNDARY
      HITBDY  = .true.
      go to 160
 200  call STGGET(TNOW,TRI, TRIANG)
      do 210 IPCENT = 1,3
        if (TRI(IPCENT + 3).eq.JP) go to 220
 210  continue
C
C                                       IDENTIFY NEXT POINT
 220  I = SUB1(IPCENT)
      KP= TRI(I + 3)
C
C                   RECORD INFORMATION ABOUT  THE NEW EDGE IN
C                   THE LIST STRUCTURE.
C
C            POP AVAIL TO L
      L = AVAIL
      AVAIL = KE(AVAIL,1)
      KE(L,2)= LASTL
      KE(L,3)= TNOW
      I      = ADD1(IPCENT)
      KE(L,4)= TRI(I)
      KE(LASTL,1) = L
      LASTL  = L
      JJ1 = JJ2
      JJ2 = KP
      DSQ(L) = DSFCN(X(JJ1),Y(JJ1),X(JJ2),Y(JJ2) )
      if (KP .eq. J1) go to 240
C                                       KP IS A NEW NEIGHBORING POINT.
      MFIT = MFIT +1
C                   BUILD ROW
C1000 FORMAT(25H BUILD ROW..  CENTER PT =,I5,12H,  EQUA NO =,I3,
C    *13H  NEARBY PT =,I5)
C        WRITE(*,1000) JP,MFIT,KP
      JUSED(MFIT) = KP
      W(MFIT,4) = X(KP) - X0
      W(MFIT,5) = Y(KP) - Y0
      W(MFIT,6) = Z(KP) - Z0
      W(MFIT,1) = W(MFIT,4)**2
      W(MFIT,2) = W(MFIT,4) * W(MFIT,5)
      W(MFIT,3) = W(MFIT,5)**2
      KE(L,1) = 0
      go to 160
C
C                   KP IS NOT A NEW POINT.  WE HAVE CYCLED ALL THE WAY
C                   AROUND TO THE INITIAL POINT J1.  THIS MEANS THE
C                   ENTIRE RING OF NEIGHBORS HAS BEEN FOUND.
C
 240  KE(L,1) = FIRSTL
      KE(FIRSTL,2) = LASTL
C
C         ABOVE LOOP IS EXITED FOR ONE OF THREE REASONS..
C              (1) MFIT .EQ. MFMAX,
C              (2) HITBDY = .TRUE.,
C           OR (3) COMPLETE RING OF NEIGHBORS HAS BEEN BUILT
C
      go to 160
C
C     **************************
C                              *
C     PROCEDURE ( MOVE RIGHT )
C                              *
C      *************************
C                                      RIGHT MEANS CLOCKWISE.
 360  JJ1 = J1
      IPCENT = IVERT
 370  if (MFIT .ge. MFMAX) go to 90
      TNOW = TRI(IPCENT)
C
C              THE FOLLOWING TEST CAUSES AN EXIT FROM THE LOOP IF THE
C              BOUNDARY ON THE RIGHT HAS BEEN ENCOUNTERED.
C
      if (TNOW .eq. 0) go to 90
      call STGGET(TNOW,TRI, TRIANG)
      do 380 IPCENT = 1,3
        if (TRI(IPCENT + 3) .eq. JP) go to 390
 380  continue
C                                       IDENTIFY NEXT POINT
 390  IP1 = ADD1(IPCENT)
      KP = TRI(IP1 + 3)
C
C                   RECORD INFO ABOUT THE NEW EDGE IN THE LIST STRUCTURE
C
C            POP AVAIL TO L
      L = AVAIL
      AVAIL = KE(AVAIL,1)
      KE(L,1) = FIRSTL
      KE(L,2) = 0
      KE(L,3) = TNOW
      KE(L,4) = TRI(IP1)
      KE(FIRSTL,2) = L
      FIRSTL  = L
      JJ2 = JJ1
      JJ1 = KP
      DSQ(L) = DSFCN(X(JJ1),Y(JJ1),X(JJ2),Y(JJ2) )
C
C                                       KP IS A NEW NEIGHBORING POINT.
      MFIT = MFIT + 1
C                   BUILD ROW
C1000 FORMAT(25H BUILD ROW..  CENTER PT =,I5,12H,  EQUA NO =,I3,
C    *13H  NEARBY PT =,I5)
C        WRITE(*,1000) JP,MFIT,KP
      JUSED(MFIT) = KP
      W(MFIT,4) = X(KP) - X0
      W(MFIT,5) = Y(KP) - Y0
      W(MFIT,6) = Z(KP) - Z0
      W(MFIT,1) = W(MFIT,4)**2
      W(MFIT,2) = W(MFIT,4) * W(MFIT,5)
      W(MFIT,3) = W(MFIT,5)**2
      go to 370
C
C                   ABOVE LOOP IS EXITED WHEN
C                      (1) MFIT .EQ. MFMAX
C                   OR (2) HIT BOUNDARY ON RIGHT
C
      END
c     ==================================================================
      SUBROUTINE sTGMOR(X,Y,Z,NP,TRIANG,W)
C>> 1995-09-26 CLL Editing for inclusion into MATH77.
C
C     C.L.LAWSON, JPL, 1977 JAN 4
C         CHANGED NOV 30, 1977    TT
C         THIS SUBR IS CALLED BY STGPD AS FOLLOWS:
C      (1) IF NEWPT = .TRUE.,
C            THE SUBR EXTENDS THE RING OF NEARBY POINTS AROUND
C          POINT (X0,Y0) UP TO A TOTAL OF MIN(MFMIN,NP-1) NEARBY POINTS
C          TO BE USED IN COMPUTING A LOCAL FIT TO THE SURFACE NEAR POINT
C          (X0,Y0).
C      (2) IF NEWPT = .FALSE.,
C            THE SUBR GETS ONLY ONE MORE NEARBY POINT AND
C          RETURNS.
c     ------------------------------------------------------------------
      common  /sTGCM1/ X0, Y0, Z0, DSQ,
     *  AVAIL, MFMIN, MFIT, FIRSTL, LASTL, KE, JUSED, NEWPT
      integer mfmax
      parameter(mfmax = 16)
      real             X0, Y0, Z0, DSQ(mfmax)
      integer AVAIL,MFMIN,MFIT, FIRSTL, LASTL, KE(mfmax,4), JUSED(mfmax)
      logical NEWPT
C
      logical TSTLP1,TSTLP2,TSTLM1,TSTLM2,AVOK, TEST, USED
      integer ADD1(3), i, im1, ip1, k, kp
      integer l, limit2, lm1, lm2, lmin, lnew
      integer lp1, lp2, lp3
      integer NP
      integer SUB1(3),TNOW, TRI(6),TRIANG(1:*),V1,V2
c++ Code for .C. is INACTIVE
C      real             dmin, dsfcn, W(21,6)
c++ Code for ~.C. is ACTIVE
      integer mw
      parameter(mw = mfmax+5)
      real             dmin, dsfcn, W(mw,6)
c++ END
      real             X(NP), x1, x2, Y(NP), y1, y2, Z(NP)
C
      data ADD1(1),ADD1(2),ADD1(3),SUB1(1),SUB1(2),SUB1(3)/2,3,1, 3,1,2/
C
C     NAMELIST/DEBUG/AVAIL,MFMIN,MFIT,FIRSTL,LASTL, JUSED,
C    *X0,Y0,Z0,DSQ
c     ------------------------------------------------------------------
      DSFCN(X1,Y1,X2,Y2) = ((X1+X2)*.5 - X0)**2 + ((Y1+Y2)*.5 - Y0)**2
c     ------------------------------------------------------------------
C
      if (NEWPT) then
        LIMIT2 = MIN( MFMIN, NP - 1 )
      else
        LIMIT2 = MFIT + 1
C                        WHEN sTGMOR IS CALLED WITH
C                        NEWPT = .FALSE. MFIT MUST BE LESS THAN MFMIN
C                        AND NP-1 SO THAT LIMIT2 CAN SAFELY BE SET
C                        TO MFIT+1 AS ABOVE.
      end if
      AVOK = .true.
      go to 200

C
C     *******************************************
C                                               *
C     PROCEDURE ( DELETE EDGES L,LP1, AND LP2 )
C                                               *
C      ******************************************
C
 10   LM1 = KE(L,2)
      LP3 = KE(LP2,1)
      if (LM1 .ne. 0) KE(LM1,1) = LP3
      if (LP3 .ne. 0) KE(LP3,2) = LM1
      if ( L .eq. FIRSTL) then
        FIRSTL = LP3
      else if ((LP1 .eq. FIRSTL)  .or.  (LP2 .eq. FIRSTL)) then
        FIRSTL = LP3
        LASTL  = LM1
      else if (LP2 .eq. LASTL) then
        LASTL = LM1
      end if
C                   TRANSFER L, LP1, AND LP2 TO AVAILABLE SPACE LIST
C
      KE(LP2,1) = AVAIL
      AVAIL = L
 20   FORMAT(60X,'MORPTS..  DELETING THREE EDGES.',
     *  '  TNOW=', I5)
      WRITE(*,20) TNOW
      go to 200
C                                                                 *
C     DELETE EDGE L.  ADD TWO NEW EDGES. /KP,AVOK/
 50   continue
C             POP AVAIL TO LNEW. SET AVOK
C                                              *
C      *****************************************
      if ( AVAIL .eq. 0) then
        AVOK = .false.
C                                            ERROR CONDITION
C
 60     FORMAT(/' SUBR sTGMOR..  WARNING..  DIMENSION OF KE(,) NOT',
     *    ' LARGE ENOUGH TO FIND MFMIN NEARBY POINTS FOR ESTIMATING',
     *    ' PARTIAL DERIVS.'/
     *    ' PROGRAM WILL CONTINUE USING FEWER NEARBY POINTS.')
        WRITE(*,60)
      else
        LNEW = AVAIL
        AVAIL = KE(AVAIL,1)
      end if

      if (AVOK) then
        KE(L,1) = LNEW
        KE(LNEW,1) = LP1
        KE(LNEW,2) = L
        KE(LNEW,3) = TNOW
        if (LP1 .ne. 0) KE(LP1,2) = LNEW
        if (L .eq. LASTL) LASTL = LNEW
        call STGGET( TNOW,TRI, TRIANG)
        do 80 I = 1,3
          if (TRI(I) .eq. KE(L,3)) go to 90
 80     continue
 90     IP1 = ADD1(I)
        KE(L,4) = TRI(IP1)
        KE(L,3) = TNOW
        IM1 = SUB1(I)
        KE(LNEW,4) = TRI(IM1)
        KP         = TRI(IM1 + 3)
        V1 = TRI(IP1 + 3)
        DSQ(L)    = DSFCN( X(V1),Y(V1), X(KP),Y(KP) )
        V1 = TRI(I + 3)
        DSQ(LNEW) = DSFCN( X(V1),Y(V1), X(KP),Y(KP) )
      end if
C            PROCESS NEW VERTEX
C
C         HERE THE TRIANGLE TNOW HAS ONLY EDGE L IN COMMON WITH THE
C         CURRENT RING  AROUND POINT JP.  THE VERTEX OPPOSITE SIDE L
C         IN TRIANGLE TNOW IS KP.
C            TEST TO SEE IF POINT KP HAS ALREADY BEEN USED TO
C         BUILD AN EQUATION.
C
      USED = .false.
      do 160 K = 1,MFIT
        if (JUSED(K) .eq. KP) then
          USED = .true.
          go to 170
        end if
 160  continue
C
 170  if (USED) go to 200
      MFIT = MFIT + 1
C               BUILD ROW
C1000    FORMAT(' BUILD ROW IN sTGMOR..   ',5X,''   EQUA NO =',I3,
C    *   '  NEARBY PT =',I5)
C        WRITE(*,1000)    MFIT,KP
      JUSED(MFIT) = KP
      W(MFIT,4) = X(KP) - X0
      W(MFIT,5) = Y(KP) - Y0
      W(MFIT,6) = Z(KP) - Z0
      W(MFIT,1) = W(MFIT,4)**2
      W(MFIT,2) = W(MFIT,4) * W(MFIT,5)
      W(MFIT,3) = W(MFIT,5)**2

c
 200  if ((MFIT .ge. LIMIT2) .or. (.not. AVOK)) return
C        WRITE(*,DEBUG)
C1004    FORMAT(/' KE(,)='/(7X,5I5))
C        WRITE(*,1004) (I,(KE(I,J), J=1,4), I=1,16)
C
C                   FIND EDGE WHOSE MIDPOINT IS CLOSEST TO POINT JP.
C                   SKIP BOUNDARY EDGES.
      L = FIRSTL
      LMIN = 0
 220  if (KE(L,4) .ne. 0) then
        if (LMIN .eq. 0) then
          TEST = .true.
        else
          TEST = DSQ(L) .lt. DMIN
        end if
        if (TEST) then
          DMIN = DSQ(L)
          LMIN = L
        end if
      end if
      if (L .eq. LASTL) go to 230
      L = KE(L,1)
      go to 220
 230  L = LMIN
C
C         EDGE L IS THE NEAREST TO POINT JP.
C         INVESTIGATE TRIANGLE ON OPPOSIDE SIDE OF EDGE L.
C
      TNOW = KE(L,4)
      LP1 = KE(L,1)
      TSTLP1 = LP1 .ne. 0
      if (TSTLP1) TSTLP1 = KE(LP1,4) .eq. TNOW
      LM1 = KE(L,2)
      TSTLM1 = LM1 .ne. 0
      if (TSTLM1) TSTLM1 = KE(LM1,4) .eq. TNOW
      if (.not. TSTLP1) go to 270
      LP2 = KE(LP1,1)
      TSTLP2 = LP2 .ne. 0
      if (TSTLP2) TSTLP2 = KE(LP2,4) .eq. TNOW
      if (.not. TSTLP2) go to 250
      go to 10
 250  if (.not. TSTLM1) go to 300
      LP2 = LP1
      LP1 = L
      L  =LM1
      go to 10
 270  if (.not. TSTLM1) go to 50
      LM2 = KE(LM1,2)
      TSTLM2 = LM2 .ne. 0
      if (TSTLM2) TSTLM2 = KE(LM2,4) .eq. TNOW
      if (.not. TSTLM2) go to 280
      LP2 = L
      LP1 = LM1
      L   = LM2
      go to 10
 280  LP1 = L
      L  =LM1
C                        END OF TOP LEVEL CODE.  PROCEDURES FOLLOW..
C
C     ****************************************************
C                                                        *
C     PROCEDURE ( DELETE EDGES L AND LP1.  ADD NEW EDGE)
C                                                        *
C      ***************************************************
C
 300  LP2 = KE(LP1,1)
      KE(L,1) = LP2
      if (LP2 .ne. 0) KE(LP2,2) = L
C
      if (LP1 .eq.  FIRSTL) then
        FIRSTL = LP2
      else if (LP1 .eq. LASTL) then
        LASTL = L
      end if
C
      KE(LP1,1) = AVAIL
      AVAIL = LP1
      call STGGET( TNOW,TRI, TRIANG)
C
      do 310 I = 1,3
        if (TRI(I) .eq. KE(L,3)) go to 320
 310  continue
 320  KE(L,3) = TNOW
      IP1 = ADD1(I)
      KE(L,4) = TRI(IP1)
      if (KE(L,4) .ne. 0) then
        IM1 = SUB1(I)
        V1  = TRI(IP1 + 3)
        V2  = TRI(IM1 + 3)
        DSQ(L) = DSFCN(X(V1),Y(V1),  X(V2),Y(V2) )
      end if
C1003    FORMAT(60X,'MORPTS..  DELETING TWO EDGES, ADDING ONE.',
C    *   ' TNOW=', I5)
C        WRITE(*,1003) TNOW
      go to 200
C
      END
c     ==================================================================
      SUBROUTINE sTGLS (W,MW,MFIT,NEWPT, LIMIT, DZ1,DZ2)
C>> 1995-09-26 CLL Editing for inclusion into MATH77.
C
C     WEIGHT, SCALE, STABILIZE, AND SOLVE A LEAST SQRS PROBLEM.
C
C     C.L.LAWSON, JPL, 1976 FEB 27
C         CHANGED NOV 28, 1977    TT
C     EDITED CODE. NO CHANGE IN ALGORITHM. CLL. 1979 MAR 5.
C
C     ON INPUT, WITH NEWPT TRUE, COLUMNS OF W(,) CONTAIN TRANSLATED
C     VALUES OF X**2, X*Y, Y**2, X, Y, Z.
C
C          THE NUMBER, STAB, IS USED AS A MARQUARDT STABILIZATION
C     VALUE TO DAMP DOWN THE VALUES OF THE THREE SECOND PARTIAL
C     DERIVATIVES WHEN THE AVAILABLE DATA DOES NOT PRODUCE A
C     REASONABLY WELL-CONDITIONED LEAST SQUARES PROBLEM.
C     A LARGER VALUE OF STAB DAMPS THE SECOND PARTIALS MORE.
C
C     THE OUTPUT IS DZ1 AND DZ2.
c     ------------------------------------------------------------------
      real             AVE,WT,WT1,WT2
      save             AVE,WT,WT1,WT2
      integer limit
C
      integer i, IROW, J, JJ, KP, LEN, MFIT, MW
      real             AVESQ, C1, C2, COSINE, DMIN, DZ1, DZ2
c++ Code for .C. is INACTIVE
C      real             SINE, STAB, SUM, W(21,6)
c++ Code for ~.C. is ACTIVE
      real             SINE, STAB, SUM, W(MW,6)
c++ END
      logical   NEWPT
c     ------------------------------------------------------------------
      STAB = 1.E0
      if (.not. NEWPT) go to 200
      SUM = 0.0e0
      KP = MFIT
      do 100 I = 1,KP
        SUM = SUM + W(I,1) + W(I,3)
 100  continue
c                          KP  will be converted to floating point.
      AVESQ = SUM/KP
      AVE   = SQRT(AVESQ)
C                          WE HAVE EXPERIMENTED WITH A VARIETY OF
C                          WEIGHTING FORMULAS AND FOUND NONE SIG-
C                          NIFICANTLY BETTER THAN THE UNIFORM
C                          WEIGHTING..  WT = 1.
      WT = 1.0e0
      WT1 = WT / AVE
      WT2 = WT / AVESQ
      IROW = 1
      go to 130
C                  SCALE ROW IROW
 120  W(IROW,1) = W(IROW,1) * WT2
      W(IROW,2) = W(IROW,2) * WT2
      W(IROW,3) = W(IROW,3) * WT2
      W(IROW,4) = W(IROW,4) * WT1
      W(IROW,5) = W(IROW,5) * WT1
      W(IROW,6) = W(IROW,6) * WT
      IROW = IROW + 1
 130  if (KP .ge. IROW) go to 120
C           GIVENS TRIANGULARIZATION
      do 160 J = 1,6
        if (J + 1 .LE. KP) then
          LEN = 6 - J
          do 150 I = J + 1,KP
            if (W(I,J) .ne. 0.0e0) then
              call SROTG (W(J,J),W(I,J),COSINE,SINE)
              W(I,J) = 0.0e0
              if (LEN .gt. 0)
     *          call SROT(LEN,W(J,J+1),MW,W(I,J+1),MW,COSINE,SINE)
            end if
 150      continue
        end if
 160  continue
      go to 230
C
 200  IROW = MFIT
C                  SCALE NEW ROW
      W(IROW,1) = W(IROW,1) * WT2
      W(IROW,2) = W(IROW,2) * WT2
      W(IROW,3) = W(IROW,3) * WT2
      W(IROW,4) = W(IROW,4) * WT1
      W(IROW,5) = W(IROW,5) * WT1
      W(IROW,6) = W(IROW,6) * WT
C
C             THE NEW ROW WILL BE ACCUMULATED INTO THE TRIANGLE
C             USING GIVENS ROTATIONS.
      do 220 J = 1,6
        LEN = 6 - J
        if (W(IROW,J) .ne. 0.0e0) then
          call SROTG (W(J,J), W(IROW,J), COSINE,SINE)
          W(IROW,J) = 0.0e0
          if (LEN .gt. 0)
     *      call SROT (LEN,W(J,J+1),MW,W(IROW,J+1),MW,COSINE,SINE)
        end if
 220  continue
c
 230  if (KP .ge. 5) then
        DMIN = min( min( min(abs(W(1,1)), abs(W(2,2))),
     *    min(abs(W(3,3)), abs(W(4,4)))),
     *    abs(W(5,5)) )
      else
        DMIN = 0.0e0
      end if
C               TEST CONDITION OF SYSTEM.
      if (DMIN .ge. 0.01e0) go to 300
C
C               SYSTEM IS ILL-CONDITIONED.
C
      if (MFIT .lt. LIMIT) then
C
C               RETURN TO ADD ONE MORE POINT TO THE FITTING PROBLEM
C               IN HOPES THAT THIS WILL IMPROVE THE CONDITION.
C
        NEWPT = .false.
        return
      end if
C
C               CANNOT ADD ANY MORE POINTS.  STABILIZE BY DAMPING
C               THE SECOND PARTIAL DERIVATIVES.
C
      JJ = 1
      go to 280
c
 250  continue
C             THE NEW ROW WILL BE ACCUMULATED INTO THE TRIANGLE
C             USING GIVENS ROTATIONS.
      do 260 J = 1,6
        LEN = 6 - J
        if (W(IROW,J) .ne. 0.0e0) then
          call SROTG (W(J,J), W(IROW,J), COSINE,SINE)
          W(IROW,J) = 0.0e0
          if (LEN .gt. 0)
     *      call SROT (LEN,W(J,J+1),MW,W(IROW,J+1),MW,COSINE,SINE)
        end if
 260  continue

      JJ = JJ + 1
 280  if (JJ .gt. 3) go to 300
      IROW = MFIT + JJ
      do 290 J = 1,6
        W(IROW,J) = 0.0e0
 290  continue
      W(IROW,JJ) = STAB
      go to 250
C
C                 SOLVE FOR DZ1 AND DZ2
C
 300  C2 = W(5,6) / W(5,5)
      C1 = (W(4,6) - C2 * W(4,5)) / W(4,4)
      DZ1 = C1 / AVE
      DZ2 = C2 / AVE
      NEWPT = .true.
C
      return
      end


text:
      SUBROUTINE PCC_SIN(NR,ICORE,AC,BC,CC,R,CDC)

ccccccccccccccccccccccccccccccccccccccccccc
cccccc   Replace the original full core charge CDC/r inside r(icore),
cccccc   with a form ac*sin(bc*r), match the original charge at r(icore)
cccccc   for its value and slop.  The purpose is to generate a smoother
cccccc   core charge to be used in plane wave calculation
cccccc
cccccc   Lin-Wang Wang, Feb. 28, 2001
cccccc    Actually, it should be called: PCC_SIN(xxxx)
cccccc
cccccc
ccccccccc  This subroutine  returns a pseudo core  of
ccccccccc  cdc(r)=ac*r*sin(bc*r)   inside r(icore)
ccccccccc  instead of the original
ccccccccc  cdc(r)=r^2 * exp(ac+bc*r^2+cc*r^4)  inside r(icore)
ccccccccc
ccccccccc  This new pseudo-core is softer than
ccccccccc  the original exp. core, and uses smaller Ecut in
ccccccccc  PEtot calculation. For PEtot calculation, we suggest
ccccccccc  this pseudo-core for "pe" calculation (for no GGA pseudopotential).
ccccccccc  Lin-Wang Wang, Feb. 28, 2001

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      PARAMETER (NN = 5)

      DIMENSION R(NR),CDC(NR),
     .          DR_K(-NN:NN),DDR_K(-NN:NN),
     .          CDC_SCA(ICORE-NN:ICORE+NN)

            tpfive=2.5d0
            pi=4*datan(1.d0)

cccccccccccccccccccccccccccccccccccccccccccc

            cdcp = (cdc(icore+1)-cdc(icore))/(r(icore+1)-r(icore))
            tanb = cdc(icore)/(r(icore)*cdcp-cdc(icore))
            rbold = tpfive
            do 500 i = 1, 50
               rbnew = pi + datan(tanb*rbold)
               if (abs(rbnew-rbold) .lt. .00001D0) then
                  bc = rbnew/r(icore)
                  ac = cdc(icore)/(r(icore)*sin(rbnew))
                  do 490 j = 1, icore
                     cdc(j) = ac*r(j)*sin(bc*r(j))
  490             continue
           write(*,*) "r(icore),ac,bc", r(icore), ac, bc
           write(*,*) "**************"
           write(*,*) "**** USING sin(bc*r) for core charge ***"
           write(*,*) "**************"
c
                  go to 530
c
               else
                  rbold = rbnew
               end if
  500       continue
            write(*,*) "SHOULD NEVER BE HERE, fail to find core charge"
            Write(*,*) "inside PCC_EXP.f, stop"
            stop

  530       continue
            cc=1.d0

           return
           end

ext: f

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!

module mod

type base_t
contains
procedure, pass   :: baseproc_pass => baseproc
procedure, nopass :: baseproc_nopass => baseproc
generic           :: some_proc => baseproc_pass, baseproc_nopass
end type

type, extends(base_t) :: ext_t
end type

contains

logical function baseproc(this)
class(base_t) :: this
select type(this)
type is(base_t)
baseproc = .true.
type is (ext_t)
baseproc = .false.
class default
stop 'baseproc: unexepected type for this'
end select
end function

end module

program p
USE CHECK_MOD
use mod
logical results(2)
logical expect(2)
data results /.false.,.true./
data expect /.true.,.false./
type(base_t) :: t
type(ext_t) :: t2

  results(1) = t%some_proc()
  results(2) = t%some_proc(t2)

  call check(results,expect,2)

end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
      MODULE PCGMODULE
        INTEGER,SAVE,POINTER  ::ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER
        integer,save,pointer  :: mxiterpcg
        REAL   ,SAVE,POINTER  ::HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG
        REAL   ,SAVE,POINTER  ::DAMPPCGT
        INTEGER,SAVE,POINTER  ::IHCOFADD
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:,:,:) ::VPCG
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:,:,:) ::SS
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:,:,:) ::P
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:,:,:) ::HPCG
        REAL,             SAVE, POINTER, DIMENSION(:,:,:) ::CD
        REAL,             SAVE, POINTER, DIMENSION(:,:,:) ::HCSV
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::LHCH
        REAL,             SAVE, POINTER, DIMENSION(:)     ::HCHG
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::LRCHPCG
        REAL,             SAVE, POINTER, DIMENSION(:)     ::RCHG
        INTEGER,          SAVE, POINTER, DIMENSION(:)     ::IT1
      TYPE PCGTYPE
        INTEGER,POINTER  ::ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER
        integer,pointer  :: mxiterpcg
        REAL   ,POINTER  ::HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG
        REAL   ,POINTER  :: DAMPPCGT
        INTEGER,POINTER  ::IHCOFADD
        DOUBLE PRECISION,  POINTER, DIMENSION(:,:,:) ::VPCG
        DOUBLE PRECISION,  POINTER, DIMENSION(:,:,:) ::SS
        DOUBLE PRECISION,  POINTER, DIMENSION(:,:,:) ::P
        DOUBLE PRECISION,  POINTER, DIMENSION(:,:,:) ::HPCG
        REAL,              POINTER, DIMENSION(:,:,:) ::CD
        REAL,              POINTER, DIMENSION(:,:,:) ::HCSV
        INTEGER,           POINTER, DIMENSION(:,:)   ::LHCH
        REAL,              POINTER, DIMENSION(:)     ::HCHG
        INTEGER,           POINTER, DIMENSION(:,:)   ::LRCHPCG
        REAL,              POINTER, DIMENSION(:)     ::RCHG
        INTEGER,           POINTER, DIMENSION(:)     ::IT1
      END TYPE
      TYPE(PCGTYPE), SAVE ::PCGDAT(10)
      END MODULE PCGMODULE

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

      SUBROUTINE PCG7AR(IN,MXITER,IGRID)
C     ******************************************************************
C     ALLOCATE STORAGE FOR PCG ARRAYS AND READ PCG DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,   ONLY:IOUT,NCOL,NROW,NLAY
      USE PCGMODULE,ONLY:ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER,
     1                   HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG,VPCG,SS,P,
     2                   HPCG,CD,HCSV,LHCH,HCHG,LRCHPCG,RCHG,IT1,
     3                   DAMPPCGT,mxiterpcg,
     4                   IHCOFADD  !JDH 20110814
      use utl7module, only: URDCOM, URWORD
C
      CHARACTER*200 LINE
      INTEGER IN,MXITER
      double precision :: r
C     ------------------------------------------------------------------
      ALLOCATE(ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER)
      allocate(mxiterpcg)
      ALLOCATE(HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG,DAMPPCGT)
      ALLOCATE(IHCOFADD)
C
C-------PRINT A MESSAGE IDENTIFYING PCG PACKAGE
      WRITE (IOUT,500)
  500 FORMAT (1X,/1X,'PCG -- CONJUGATE-GRADIENT SOLUTION PACKAGE',
     &        ', VERSION 7, 5/2/2005')
C
C-------READ AND PRINT COMMENTS, MXITER,ITER1 AND NPCOND
      CALL URDCOM(IN,IOUT,LINE)
      LLOC=1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXITER,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ITER1,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NPCOND,R,IOUT,IN)
C  JDH 20110814 - ADDED IHCOFADD
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IHCOFADD,R,-IOUT,IN)
      mxiterpcg = MXITER
      WRITE (IOUT,509) MXITER, ITER1, NPCOND
  509 FORMAT (' MAXIMUM OF ',I6,' CALLS OF SOLUTION ROUTINE',/,
     &        ' MAXIMUM OF ',I6,
     &        ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE',/,
     &        ' MATRIX PRECONDITIONING TYPE :',I5)
      IF(IHCOFADD.NE.0) WRITE(IOUT,510)
  510 FORMAT(' IHCOFADD option -- When an active cell is surrounded ',
     & 'by dry cells, convert',/,
     & ' the cell to dry only if storage and head-dependent boundary ',
     & 'flow are 0.')
C
C-------ALLOCATE SPACE FOR THE PCG ARRAYS
      ALLOCATE (VPCG(NCOL,NROW,NLAY))
      ALLOCATE (SS(NCOL,NROW,NLAY))
      ALLOCATE (P(NCOL,NROW,NLAY))
      ALLOCATE (HPCG(NCOL,NROW,NLAY))
      ALLOCATE (CD(NCOL,NROW,NLAY))
      IF(NPCOND.EQ.2) THEN
         ALLOCATE (HCSV(NCOL,NROW,NLAY))
      ELSE
         ALLOCATE (HCSV(1,1,1))
      END IF
      ITMEM=MXITER*ITER1
      ALLOCATE (HCHG(ITMEM))
      ALLOCATE (LHCH(3,ITMEM))
      ALLOCATE (RCHG(ITMEM))
      ALLOCATE (LRCHPCG(3,ITMEM))
      ALLOCATE (IT1(ITMEM))
C
C-------READ HCLOSEPCG,RCLOSEPCG,RELAXPCG,NBPOL,IPRPCG,MUTPCG
Crgn changed read to check for different value for transient dampening
      READ (IN,*) HCLOSEPCG,RCLOSEPCG,RELAXPCG,
     1    NBPOL,IPRPCG,MUTPCG,DAMPPCG
      IF ( DAMPPCG.LT.0.0 ) THEN
        BACKSPACE IN
        READ (IN,*) HCLOSEPCG,RCLOSEPCG,RELAXPCG,
     1    NBPOL,IPRPCG,MUTPCG,DAMPPCG,DAMPPCGT
          DAMPPCG = -DAMPPCG
          IF (DAMPPCGT.EQ.0.0) DAMPPCGT = 1.0
      ELSE
        IF (DAMPPCG.EQ.0.0) DAMPPCG = 1.0
        DAMPPCGT = DAMPPCG
      END IF
C
C-------PRINT MXITER,ITER1,NPCOND,HCLOSEPCG,RCLOSEPCG,RELAXPCG,NBPOL,IPRPCG,
C-------MUTPCG,DAMPPCG
        WRITE (IOUT,511)
  511   FORMAT (1X,///,36X,'SOLUTION BY THE CONJUGATE-GRADIENT METHOD',
     &        /,35X,43('-'))
      WRITE (IOUT,512) MXITER
  512 FORMAT (1X,19X,'MAXIMUM NUMBER OF CALLS TO PCG ROUTINE =',I9)
      WRITE (IOUT,515) ITER1
  515 FORMAT (1X,23X,'MAXIMUM ITERATIONS PER CALL TO PCG =',I9)
      WRITE (IOUT,520) NPCOND
  520 FORMAT (1X,30X,'MATRIX PRECONDITIONING TYPE =',I9)
      IF (NPCOND.EQ.2) WRITE (IOUT,525)
  525 FORMAT (1X,53X,'THE MATRIX WILL BE SCALED')
      WRITE (IOUT,530) RELAXPCG, NBPOL
  530 FORMAT (1X,7X,'RELAXATION FACTOR (ONLY USED WITH',
     &        ' PRECOND. TYPE 1) =',E15.5,/,1X,
     &        'PARAMETER OF POLYNOMIAL PRECOND.',
     &        ' = 2 (2) OR IS CALCULATED :',I9)
      WRITE (IOUT,535) HCLOSEPCG
  535 FORMAT (1X,24X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)
      WRITE (IOUT,540) RCLOSEPCG
  540 FORMAT (1X,20X,'RESIDUAL CHANGE CRITERION FOR CLOSURE =',E15.5)
      IF (IPRPCG.LE.0) IPRPCG = 999
      WRITE (IOUT,545) IPRPCG, MUTPCG
  545 FORMAT (1X,11X,'PCG HEAD AND RESIDUAL CHANGE PRINTOUT INTERVAL =',
     &        I9,/,1X,4X,
     &        'PRINTING FROM SOLVER IS LIMITED(1) OR SUPPRESSED (>1) =',
     &        I9)
      WRITE (IOUT,550) DAMPPCG,DAMPPCGT
  550 FORMAT (1X,27X,'STEADY-STATE DAMPING PARAMETER =',E15.5
     &       /1X,30X,'TRANSIENT DAMPING PARAMETER =',E15.5)
      NITER = 0
C
      CALL PCG7PSV(IGRID)
      RETURN
      END SUBROUTINE PCG7AR

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

      SUBROUTINE PCG7DA(IGRID)
C  Deallocate PCG DATA
      USE PCGMODULE
C
      CALL PCG7PNT(IGRID)
        DEALLOCATE(ITER1,NPCOND,NBPOL,IPRPCG,MUTPCG,NITER)
        deallocate(mxiterpcg)
        DEALLOCATE(HCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG,DAMPPCGT)
        DEALLOCATE(IHCOFADD)
        DEALLOCATE(VPCG)
        DEALLOCATE(SS)
        DEALLOCATE(P)
        DEALLOCATE(HPCG)
        DEALLOCATE(CD)
        DEALLOCATE(HCSV)
        DEALLOCATE(LHCH)
        DEALLOCATE(HCHG)
        DEALLOCATE(LRCHPCG)
        DEALLOCATE(RCHG)
        DEALLOCATE(IT1)
C
      RETURN
      END SUBROUTINE PCG7DA

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

      SUBROUTINE PCG7PNT(IGRID)
C  Set pointers to PCG data for a grid
      USE PCGMODULE
C
      ITER1=>PCGDAT(IGRID)%ITER1
      NPCOND=>PCGDAT(IGRID)%NPCOND
      NBPOL=>PCGDAT(IGRID)%NBPOL
      IPRPCG=>PCGDAT(IGRID)%IPRPCG
      MUTPCG=>PCGDAT(IGRID)%MUTPCG
      NITER=>PCGDAT(IGRID)%NITER
      HCLOSEPCG=>PCGDAT(IGRID)%HCLOSEPCG
      RCLOSEPCG=>PCGDAT(IGRID)%RCLOSEPCG
      RELAXPCG=>PCGDAT(IGRID)%RELAXPCG
      DAMPPCG=>PCGDAT(IGRID)%DAMPPCG
      DAMPPCGT=>PCGDAT(IGRID)%DAMPPCGT
      IHCOFADD=>PCGDAT(IGRID)%IHCOFADD
      VPCG=>PCGDAT(IGRID)%VPCG
      SS=>PCGDAT(IGRID)%SS
      P=>PCGDAT(IGRID)%P
      HPCG=>PCGDAT(IGRID)%HPCG
      CD=>PCGDAT(IGRID)%CD
      HCSV=>PCGDAT(IGRID)%HCSV
      LHCH=>PCGDAT(IGRID)%LHCH
      HCHG=>PCGDAT(IGRID)%HCHG
      LRCHPCG=>PCGDAT(IGRID)%LRCHPCG
      RCHG=>PCGDAT(IGRID)%RCHG
      IT1=>PCGDAT(IGRID)%IT1
      mxiterpcg => pcgdat(igrid)%mxiterpcg
C
      RETURN
      END SUBROUTINE PCG7PNT

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

      SUBROUTINE PCG7PSV(IGRID)
C  Save pointers to PCG data
      USE PCGMODULE
C
      PCGDAT(IGRID)%ITER1=>ITER1
      PCGDAT(IGRID)%NPCOND=>NPCOND
      PCGDAT(IGRID)%NBPOL=>NBPOL
      PCGDAT(IGRID)%IPRPCG=>IPRPCG
      PCGDAT(IGRID)%MUTPCG=>MUTPCG
      PCGDAT(IGRID)%NITER=>NITER
      PCGDAT(IGRID)%HCLOSEPCG=>HCLOSEPCG
      PCGDAT(IGRID)%RCLOSEPCG=>RCLOSEPCG
      PCGDAT(IGRID)%RELAXPCG=>RELAXPCG
      PCGDAT(IGRID)%DAMPPCG=>DAMPPCG
      PCGDAT(IGRID)%DAMPPCGT=>DAMPPCGT
      PCGDAT(IGRID)%IHCOFADD=>IHCOFADD
      PCGDAT(IGRID)%VPCG=>VPCG
      PCGDAT(IGRID)%SS=>SS
      PCGDAT(IGRID)%P=>P
      PCGDAT(IGRID)%HPCG=>HPCG
      PCGDAT(IGRID)%CD=>CD
      PCGDAT(IGRID)%HCSV=>HCSV
      PCGDAT(IGRID)%LHCH=>LHCH
      PCGDAT(IGRID)%HCHG=>HCHG
      PCGDAT(IGRID)%LRCHPCG=>LRCHPCG
      PCGDAT(IGRID)%RCHG=>RCHG
      PCGDAT(IGRID)%IT1=>IT1
      pcgdat(igrid)%mxiterpcg => mxiterpcg
C
      RETURN
      END SUBROUTINE PCG7PSV

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


text:
subroutine wrm_main(lundia    ,error     ,selmap    ,grdang    ,dtsec     , &
                  & itmapc    ,runtxt    ,trifil    ,wrifou    ,initi     , &
                  & gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: wrm_main.f90 6540 2016-09-12 14:27:12Z ye $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/wrm_main.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Main routine for writing the FLOW HIS file.
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use dfparall, only: inode, master, parll, dfint, nproc
    use datagroups
    use netcdf
    use dffunctionals, only: dffind_duplicate
    !
    use globaldata
    !
    implicit none
    !
    include 'fsm.i'
    include 'tri-dyn.igd'
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                              , pointer :: nmax
    integer                              , pointer :: nmaxus
    integer                              , pointer :: mmax
    integer                              , pointer :: kmax
    integer                              , pointer :: lmax
    integer                              , pointer :: lstsci
    integer                              , pointer :: lsed
    integer                              , pointer :: lsedtot
    integer                              , pointer :: ltur
    integer                              , pointer :: noroco
    integer                              , pointer :: norow
    integer                              , pointer :: nsrc
    integer                              , pointer :: nostat
    integer                              , pointer :: ntruv
    logical                              , pointer :: lfsdu
    logical                              , pointer :: zmodel
    integer(pntrsize)                    , pointer :: ibuff
    integer(pntrsize)                    , pointer :: irocol
    integer(pntrsize)                    , pointer :: kcs
    integer(pntrsize)                    , pointer :: kcu
    integer(pntrsize)                    , pointer :: kcv
    integer(pntrsize)                    , pointer :: kfs
    integer(pntrsize)                    , pointer :: kfu
    integer(pntrsize)                    , pointer :: kfv
    integer(pntrsize)                    , pointer :: kfumin
    integer(pntrsize)                    , pointer :: kfvmin
    integer(pntrsize)                    , pointer :: kfsmin
    integer(pntrsize)                    , pointer :: kfumax
    integer(pntrsize)                    , pointer :: kfvmax
    integer(pntrsize)                    , pointer :: kfsmax
    integer(pntrsize)                    , pointer :: mnksrc
    integer(pntrsize)                    , pointer :: alfas
    integer(pntrsize)                    , pointer :: cfurou
    integer(pntrsize)                    , pointer :: cfvrou
    integer(pntrsize)                    , pointer :: cvalu0
    integer(pntrsize)                    , pointer :: cvalv0
    integer(pntrsize)                    , pointer :: dicww
    integer(pntrsize)                    , pointer :: dp
    integer(pntrsize)                    , pointer :: dps
    integer(pntrsize)                    , pointer :: dpu
    integer(pntrsize)                    , pointer :: dpv
    integer(pntrsize)                    , pointer :: dzs1
    integer(pntrsize)                    , pointer :: enstro
    integer(pntrsize)                    , pointer :: eroll1
    integer(pntrsize)                    , pointer :: evap
    integer(pntrsize)                    , pointer :: ewave1
    integer(pntrsize)                    , pointer :: fxw
    integer(pntrsize)                    , pointer :: fyw
    integer(pntrsize)                    , pointer :: gsqs
    integer(pntrsize)                    , pointer :: guu
    integer(pntrsize)                    , pointer :: gvv
    integer(pntrsize)                    , pointer :: hrms
    integer(pntrsize)                    , pointer :: p1
    integer(pntrsize)                    , pointer :: patm
    integer(pntrsize)                    , pointer :: precip
    integer(pntrsize)                    , pointer :: qxkr
    integer(pntrsize)                    , pointer :: qxkw
    integer(pntrsize)                    , pointer :: qykr
    integer(pntrsize)                    , pointer :: qykw
    integer(pntrsize)                    , pointer :: r1
    integer(pntrsize)                    , pointer :: rbuff
    integer(pntrsize)                    , pointer :: rho
    integer(pntrsize)                    , pointer :: rich
    integer(pntrsize)                    , pointer :: rtur1
    integer(pntrsize)                    , pointer :: s1
    integer(pntrsize)                    , pointer :: sig
    integer(pntrsize)                    , pointer :: taubmx
    integer(pntrsize)                    , pointer :: taubpu
    integer(pntrsize)                    , pointer :: taubpv
    integer(pntrsize)                    , pointer :: taubsu
    integer(pntrsize)                    , pointer :: taubsv
    integer(pntrsize)                    , pointer :: thick
    integer(pntrsize)                    , pointer :: u1
    integer(pntrsize)                    , pointer :: umnldf
    integer(pntrsize)                    , pointer :: v1
    integer(pntrsize)                    , pointer :: vicuv
    integer(pntrsize)                    , pointer :: vicww
    integer(pntrsize)                    , pointer :: vmnldf
    integer(pntrsize)                    , pointer :: vortic
    integer(pntrsize)                    , pointer :: w1
    integer(pntrsize)                    , pointer :: windcd
    integer(pntrsize)                    , pointer :: windu
    integer(pntrsize)                    , pointer :: windv
    integer(pntrsize)                    , pointer :: wphy
    integer(pntrsize)                    , pointer :: wrkb3
    integer(pntrsize)                    , pointer :: wrkb4
    integer(pntrsize)                    , pointer :: wsu
    integer(pntrsize)                    , pointer :: wsv
    integer(pntrsize)                    , pointer :: xcor
    integer(pntrsize)                    , pointer :: xz
    integer(pntrsize)                    , pointer :: ycor
    integer(pntrsize)                    , pointer :: yz
    integer(pntrsize)                    , pointer :: z0ucur
    integer(pntrsize)                    , pointer :: z0vcur
    integer(pntrsize)                    , pointer :: z0urou
    integer(pntrsize)                    , pointer :: z0vrou
    integer(pntrsize)                    , pointer :: namcon
    integer(pntrsize)                    , pointer :: namsrc
    integer                              , pointer :: itdate
    real(fp)                             , pointer :: dt
    real(fp)                             , pointer :: tunit
    real(fp)                             , pointer :: tzone
    logical                              , pointer :: lfbedfrm
    logical                              , pointer :: roller
    logical                              , pointer :: sferic
    logical                              , pointer :: xbeach
    character(20) , dimension(:)         , pointer :: namsed
    character(8)                         , pointer :: dpsopt
    character(4)                         , pointer :: rouflo
    integer       , dimension(:)         , pointer :: order_sta
    integer       , dimension(:)         , pointer :: order_tra
    !
    type (datagroup)                     , pointer :: group
    logical                              , pointer :: first
    real(fp)                             , pointer :: rhow
    integer                              , pointer :: ktemp
    !
    integer                              , pointer :: mfg
    integer                              , pointer :: mlg
    integer                              , pointer :: nfg
    integer                              , pointer :: nlg
    integer                              , pointer :: nmaxgl
    integer                              , pointer :: mmaxgl
    integer       , dimension(:,:)       , pointer :: iarrc
    integer       , dimension(:)         , pointer :: mf
    integer       , dimension(:)         , pointer :: ml
    integer       , dimension(:)         , pointer :: nf
    integer       , dimension(:)         , pointer :: nl
    !
    logical                              , pointer :: mergemap
!
! Global variables
!
    integer                                      :: itmapc !!  Current time counter for the map data file
    integer                                      :: lundia !  Description and declaration in inout.igs
    logical                                      :: error  !!  Flag=TRUE if an error is encountered
    real(fp)                                     :: dtsec  !!  Integration time step [in seconds]
    real(fp)                                     :: grdang !  Description and declaration in tricom.igs
    character(23)                                :: selmap !  Description and declaration in tricom.igs
    character(30) , dimension(10)                :: runtxt !!  Textual description of model input
    character(*)                   , intent(in)  :: trifil !  File name for FLOW NEFIS output files (tri"h/m"-"casl""labl".dat/def)
    logical                                      :: wrifou
    integer                                      :: initi
!
! Local variables
!
    character(10)                                     :: velt          ! Velocity type 'Eulerian' or 'GLM'
    integer                                           :: nostatgl      ! global number of stations (i.e. original number excluding duplicate stations located in the halo regions)
    integer                                           :: nostatto      ! total number of stations (including "duplicate" stations located in halo regions)
    integer                                           :: ntruvgl       ! global number of tracks (i.e. original number excluding duplicate stations located in the halo regions)
    integer                                           :: ntruvto       ! total number of tracks (including "duplicate" stations located in halo regions)
    integer                                           :: irequest
    integer                                           :: ierror
    integer                                           :: ifile
    integer                                           :: filetype
    integer                                           :: fds
    integer                                , external :: open_datdef
    integer                                , external :: clsnef
    character(256)                                    :: filename
    integer(pntrsize)                                 :: velu           ! U velocity array (FSM r-index)
    integer(pntrsize)                                 :: velv           ! V velocity array (FSM r-index)
    character(256)                                    :: version_full
    character(8)                                      :: cdate
    character(10)                                     :: ctime
    character(5)                                      :: czone
    character(1024)                                   :: error_string
    character(16)                                     :: simdat        ! Simulation date representing the flow condition at this date
    character(20)                                     :: rundat        ! Execution date of the simulation
    character(6)                                      :: ftype         ! String containing to which output file version group should be written
    !
    integer                                           :: ii
    integer                                           :: ip
    integer                                           :: istat
    integer                                           :: n
    integer                                           :: m
    character(4)                                      :: part_nr
    integer                                           :: bck_inode
    integer                                           :: bck_nproc
    logical                                           :: bck_parll
    type(dfparalltype)                  , pointer     :: bck_gdparall
    integer      , dimension(:,:)       , allocatable :: ipartition
!
!! executable statements -------------------------------------------------------
!
    wrkb3               => gdp%gdaddress%wrkb3
    wrkb4               => gdp%gdaddress%wrkb4
    nmax                => gdp%d%nmax
    nmaxus              => gdp%d%nmaxus
    mmax                => gdp%d%mmax
    kmax                => gdp%d%kmax
    lmax                => gdp%d%lmax
    lstsci              => gdp%d%lstsci
    lsed                => gdp%d%lsed
    lsedtot             => gdp%d%lsedtot
    ltur                => gdp%d%ltur
    noroco              => gdp%d%noroco
    norow               => gdp%d%norow
    nsrc                => gdp%d%nsrc
    nostat              => gdp%d%nostat
    ntruv               => gdp%d%ntruv
    zmodel              => gdp%gdprocs%zmodel
    ibuff               => gdp%gdr_i_ch%ibuff
    irocol              => gdp%gdr_i_ch%irocol
    kcs                 => gdp%gdr_i_ch%kcs
    kcu                 => gdp%gdr_i_ch%kcu
    kcv                 => gdp%gdr_i_ch%kcv
    kfs                 => gdp%gdr_i_ch%kfs
    kfu                 => gdp%gdr_i_ch%kfu
    kfv                 => gdp%gdr_i_ch%kfv
    kfumin              => gdp%gdr_i_ch%kfumin
    kfvmin              => gdp%gdr_i_ch%kfvmin
    kfsmin              => gdp%gdr_i_ch%kfsmin
    kfumax              => gdp%gdr_i_ch%kfumax
    kfvmax              => gdp%gdr_i_ch%kfvmax
    kfsmax              => gdp%gdr_i_ch%kfsmax
    mnksrc              => gdp%gdr_i_ch%mnksrc
    alfas               => gdp%gdr_i_ch%alfas
    cfurou              => gdp%gdr_i_ch%cfurou
    cfvrou              => gdp%gdr_i_ch%cfvrou
    cvalu0              => gdp%gdr_i_ch%cvalu0
    cvalv0              => gdp%gdr_i_ch%cvalv0
    dicww               => gdp%gdr_i_ch%dicww
    dp                  => gdp%gdr_i_ch%dp
    dps                 => gdp%gdr_i_ch%dps
    dpu                 => gdp%gdr_i_ch%dpu
    dpv                 => gdp%gdr_i_ch%dpv
    dzs1                => gdp%gdr_i_ch%dzs1
    enstro              => gdp%gdr_i_ch%enstro
    eroll1              => gdp%gdr_i_ch%eroll1
    evap                => gdp%gdr_i_ch%evap
    ewave1              => gdp%gdr_i_ch%ewave1
    fxw                 => gdp%gdr_i_ch%fxw
    fyw                 => gdp%gdr_i_ch%fyw
    gsqs                => gdp%gdr_i_ch%gsqs
    guu                 => gdp%gdr_i_ch%guu
    gvv                 => gdp%gdr_i_ch%gvv
    hrms                => gdp%gdr_i_ch%hrms
    p1                  => gdp%gdr_i_ch%p1
    patm                => gdp%gdr_i_ch%patm
    precip              => gdp%gdr_i_ch%precip
    qxkr                => gdp%gdr_i_ch%qxkr
    qxkw                => gdp%gdr_i_ch%qxkw
    qykr                => gdp%gdr_i_ch%qykr
    qykw                => gdp%gdr_i_ch%qykw
    r1                  => gdp%gdr_i_ch%r1
    rbuff               => gdp%gdr_i_ch%rbuff
    rho                 => gdp%gdr_i_ch%rho
    rich                => gdp%gdr_i_ch%rich
    rtur1               => gdp%gdr_i_ch%rtur1
    s1                  => gdp%gdr_i_ch%s1
    sig                 => gdp%gdr_i_ch%sig
    taubmx              => gdp%gdr_i_ch%taubmx
    taubpu              => gdp%gdr_i_ch%taubpu
    taubpv              => gdp%gdr_i_ch%taubpv
    taubsu              => gdp%gdr_i_ch%taubsu
    taubsv              => gdp%gdr_i_ch%taubsv
    thick               => gdp%gdr_i_ch%thick
    u1                  => gdp%gdr_i_ch%u1
    umnldf              => gdp%gdr_i_ch%umnldf
    v1                  => gdp%gdr_i_ch%v1
    vicuv               => gdp%gdr_i_ch%vicuv
    vicww               => gdp%gdr_i_ch%vicww
    vmnldf              => gdp%gdr_i_ch%vmnldf
    vortic              => gdp%gdr_i_ch%vortic
    w1                  => gdp%gdr_i_ch%w1
    windcd              => gdp%gdr_i_ch%windcd
    windu               => gdp%gdr_i_ch%windu
    windv               => gdp%gdr_i_ch%windv
    wphy                => gdp%gdr_i_ch%wphy
    wsu                 => gdp%gdr_i_ch%wsu
    wsv                 => gdp%gdr_i_ch%wsv
    xcor                => gdp%gdr_i_ch%xcor
    xz                  => gdp%gdr_i_ch%xz
    ycor                => gdp%gdr_i_ch%ycor
    yz                  => gdp%gdr_i_ch%yz
    z0ucur              => gdp%gdr_i_ch%z0ucur
    z0vcur              => gdp%gdr_i_ch%z0vcur
    z0urou              => gdp%gdr_i_ch%z0urou
    z0vrou              => gdp%gdr_i_ch%z0vrou
    namcon              => gdp%gdr_i_ch%namcon
    namsrc              => gdp%gdr_i_ch%namsrc
    itdate              => gdp%gdexttim%itdate
    dt                  => gdp%gdexttim%dt
    tunit               => gdp%gdexttim%tunit
    tzone               => gdp%gdexttim%tzone
    lfbedfrm            => gdp%gdbedformpar%lfbedfrm
    roller              => gdp%gdprocs%roller
    sferic              => gdp%gdtricom%sferic
    xbeach              => gdp%gdprocs%xbeach
    namsed              => gdp%gdsedpar%namsed
    dpsopt              => gdp%gdnumeco%dpsopt
    rouflo              => gdp%gdtricom%rouflo
    rhow                => gdp%gdphysco%rhow
    ktemp               => gdp%gdtricom%ktemp
    mergemap            => gdp%gdpostpr%mergemap
    lfsdu               => gdp%gdprocs%lfsdu
    !
    if (wrifou) then
       ifile = FILOUT_FOU
       filename = trifil(1:3) // 'f' // trifil(5:)
    else
       ifile = FILOUT_MAP
       filename = trifil(1:3) // 'm' // trifil(5:)
    endif
    !
    call getdatagroup(gdp, ifile, 'map-const', group)
    first               => group%first
    !
    if (first) then
       mfg                 => gdp%gdparall%mfg
       nfg                 => gdp%gdparall%nfg
       mf                  => gdp%gdparall%mf
       ml                  => gdp%gdparall%ml
       nf                  => gdp%gdparall%nf
       nl                  => gdp%gdparall%nl
       !
       allocate(ipartition(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub), stat=istat)
       if (parll) then
          do ip = 0, nproc-1
             do n = max(gdp%d%nlb,nf(ip)-nfg+1), min(nl(ip)-nfg+1,gdp%d%nub)
                do m = max(gdp%d%mlb,mf(ip)-mfg+1), min(ml(ip)-mfg+1,gdp%d%mub)
                   ipartition(n, m) = ip
                enddo
             enddo
          enddo
       else
          ipartition = 0
       endif
    endif
    !
    if (mergemap) then
        part_nr = '';
    else
        bck_inode = inode
        bck_nproc = nproc
        bck_parll = parll
        bck_gdparall => gdp%gdparall
        !
        write(part_nr,'(A,I3.3)') '-',inode
        !
        inode = 1
        nproc = 1
        parll = .false.
        !
        if (first) then
           gdp%iopartit%mfg = 1
           gdp%iopartit%mlg = mmax
           gdp%iopartit%nfg = 1
           gdp%iopartit%nlg = nmaxus
           gdp%iopartit%mmaxgl = mmax
           gdp%iopartit%nmaxgl = nmaxus
           allocate(gdp%iopartit%nf(1))
           allocate(gdp%iopartit%nl(1))
           allocate(gdp%iopartit%mf(1))
           allocate(gdp%iopartit%ml(1))
           allocate(gdp%iopartit%iarrc(4,1))
           allocate(gdp%iopartit%order_tra(1:ntruv))
           do ii = 1,ntruv
               gdp%iopartit%order_tra(ii) = ii
           enddo
           allocate(gdp%iopartit%order_sta(1:nostat))
           do ii = 1,nostat
               gdp%iopartit%order_sta(ii) = ii
           enddo
        endif
        gdp%gdparall => gdp%iopartit
    endif
    !
    mfg                 => gdp%gdparall%mfg
    mlg                 => gdp%gdparall%mlg
    nfg                 => gdp%gdparall%nfg
    nlg                 => gdp%gdparall%nlg
    mf                  => gdp%gdparall%mf
    ml                  => gdp%gdparall%ml
    nf                  => gdp%gdparall%nf
    nl                  => gdp%gdparall%nl
    iarrc               => gdp%gdparall%iarrc
    mmaxgl              => gdp%gdparall%mmaxgl
    nmaxgl              => gdp%gdparall%nmaxgl
    order_tra           => gdp%gdparall%order_tra
    order_sta           => gdp%gdparall%order_sta
    !
    call dattim(rundat    )
    simdat(1:16)  = 'yyyymmdd  hhmmss'
    simdat(1:4)   = rundat(1:4)
    simdat(5:6)   = rundat(6:7)
    simdat(7:8)   = rundat(9:10)
    simdat(11:12) = rundat(12:13)
    simdat(13:14) = rundat(15:16)
    simdat(15:16) = rundat(18:19)
    !
    call getfullversionstring_flow2d3d(version_full)
    call date_and_time(cdate, ctime, czone)
    !
    filetype = getfiletype(gdp, ifile)
    filename = trim(filename) // trim(part_nr)
    if (filetype == FTYPE_NETCDF) filename = trim(filename)//'.nc'
    !
    if (gdp%gdflwpar%flwoutput%veuler) then
       velu = wrkb3
       velv = wrkb4
       velt = 'Eulerian'
    else
       velu = u1
       velv = v1
       velt = 'GLM'
    endif
    !
    if (parll) then
       !
       ! Recalculates the effective number of stations, filtering out duplicates affected to more
       ! than one partition (i.e. located in halos)
       !
       call dfsync(gdp)
       call dffind_duplicate(lundia, nostat, nostatto, nostatgl, order_sta, gdp)
       !
       ! Recalculates the effective global number of cross sections
       !
       call dfsync(gdp)
       call dffind_duplicate(lundia, ntruv, ntruvto, ntruvgl, order_tra, gdp)
       !
       ! When order_tra points to tra_orgline, both partition-related and re-ordering-related stuff are taken care of
       !
       order_tra => gdp%gdstations%tra_orgline
    elseif (part_nr == '') then
       nostatto = nostat
       nostatgl = nostat
       ntruvto  = ntruv
       ntruvgl  = ntruv
       !
       order_sta => gdp%gdstations%sta_orgline
       order_tra => gdp%gdstations%tra_orgline
    endif
    !
    ierror = 0
    do irequest = REQUESTTYPE_DEFINE, REQUESTTYPE_WRITE
       !
       ! request REQUESTTYPE_DEFINE: define all groups, dimensions, and elements
       !         REQUESTTYPE_WRITE : write the data
       !
       if (irequest == REQUESTTYPE_DEFINE) then
          if (.not.first) cycle
          if (inode /= master) cycle
          call delnef(filename,gdp       )
       endif
       !
       ! create or open the file
       !
       if (inode /= master) then
          ! only the master needs to open the file
       elseif (first .and. irequest == REQUESTTYPE_WRITE) then
          ! the file has already been opened in step 1
       elseif (filetype == FTYPE_NEFIS) then
          if (first .and. irequest == REQUESTTYPE_DEFINE) then
             write(lundia,*) 'Creating new '//trim(filename)//'.dat'
             write(lundia,*) 'Creating new '//trim(filename)//'.def'
          endif
          ierror = open_datdef(filename ,fds      , .false.)
          if (ierror /= 0) then
             write(error_string,'(2a)') 'While trying to open dat/def-file ',trim(filename)
             call prtnefiserr(trim(error_string), gdp)
          endif
       elseif (filetype == FTYPE_NETCDF) then
          if (first .and. irequest == REQUESTTYPE_DEFINE) then
             write(lundia,*) 'Creating new '//trim(filename)
             ierror = nf90_create(filename, or(nf90_clobber,nf90_64bit_offset), fds); call nc_check_err(lundia, ierror, "creating file", filename)
             !
             ! global attributes
             !
             ierror = nf90_put_att(fds, nf90_global,  'Conventions', 'CF-1.6 SGRID-0.3'); call nc_check_err(lundia, ierror, "put_att global Conventions", filename)
             ierror = nf90_put_att(fds, nf90_global,  'institution', trim('Deltares')); call nc_check_err(lundia, ierror, "put_att global institution", filename)
             ierror = nf90_put_att(fds, nf90_global,  'references', trim('www.deltares.nl')); call nc_check_err(lundia, ierror, "put_att global references", filename)
             ierror = nf90_put_att(fds, nf90_global,  'source', trim(version_full)); call nc_check_err(lundia, ierror, "put_att global source", filename)
             ierror = nf90_put_att(fds, nf90_global,  'history', &
                    'This file is created on '//cdate(1:4)//'-'//cdate(5:6)//'-'//cdate(7:8)//'T'//ctime(1:2)//':'//ctime(3:4)//':'//ctime(5:6)//czone(1:5)// &
                    ', '//trim('Delft3D')); call nc_check_err(lundia, ierror, "put_att global history", filename)
          else
             ierror = nf90_open(filename, NF90_WRITE, fds); call nc_check_err(lundia, ierror, "opening file", filename)
          endif
       endif
       if (ierror/=0) goto 9999
       !
       ! time independent data
       !
       if (first) then
          call wrimap(lundia    ,error     ,filename  ,selmap    ,simdat    , &
                    & itdate    ,tzone     ,tunit     ,dt        ,mmax      , &
                    & kmax      ,lmax      ,lstsci    ,ltur      ,nmaxus    , &
                    & noroco    ,norow     ,nostat    ,nsrc      ,ntruv     , &
                    & grdang    ,dpsopt    ,sferic    ,lsed      ,lsedtot   , &
                    & zmodel    ,ch(namsrc),ch(namcon),namsed    , &
                    & i(kcu)    ,i(kcv)    ,i(kcs)    ,i(irocol) , &
                    & r(xcor)   ,r(ycor)   ,r(xz)     ,r(yz)     ,r(alfas)  , &
                    & r(dp)     ,r(thick)  ,r(sig)    ,r(sig)    , &
                    & d(dps)    ,r(dpu)    ,r(dpv)    ,r(gsqs)   ,wrifou    , &
                    & irequest  ,fds       ,iarrc     ,mf        ,ml        , &
                    & nf        ,nl        ,nostatto  ,nostatgl  ,order_sta , &
                    & ntruvto   ,ntruvgl   ,order_tra ,ipartition,gdp       )
          if (error) goto 9999
       endif
       !
       ! data per time step
       !
       if (wrifou) then
          ! no data per time step
       elseif (irequest==REQUESTTYPE_DEFINE .or. .not. first) then
          call wrtmap(lundia    ,error     ,filename  ,selmap    ,itmapc    , &
                    & rhow      ,mmax      ,itdate    ,dtsec     , &
                    & kmax      ,nmaxus    ,lstsci    ,ltur      , &
                    & nsrc      ,zmodel    ,i(kcs)    ,i(kfs)    ,i(kfu)    , &
                    & i(kfv)    ,i(kfumin) ,i(kfvmin) ,i(kfumax) ,i(kfvmax) , &
                    & i(kfsmin) ,i(kfsmax) ,i(mnksrc) ,r(s1)     , &
                    & d(dps)    ,r(dzs1)   ,r(thick)  ,r(windcd)  , &
                    & r(velu)   ,r(velv)   ,r(w1)     ,r(wphy)   ,r(r1)     , &
                    & r(rtur1)  ,r(taubpu) ,r(taubpv) ,r(taubsu) ,r(taubsv) , &
                    & r(vicww)  ,r(dicww)  ,r(rich)   ,r(rho)    ,r(p1)     , &
                    & r(vortic) ,r(enstro) ,r(umnldf) ,r(vmnldf) ,r(vicuv)  , &
                    & r(taubmx) ,r(windu)  ,r(windv)  ,velt      ,r(cvalu0) , &
                    & r(cvalv0) ,r(cfurou) ,r(cfvrou) ,rouflo    ,r(patm)   , &
                    & r(z0ucur) ,r(z0vcur) ,r(z0urou) ,r(z0vrou) ,ktemp     , &
                    & r(precip) ,r(evap)   ,irequest  ,fds       ,iarrc     , &
                    & mf        ,ml        ,nf        ,nl        ,gdp       )
          if (error) goto 9999
          !
          if (lsedtot>0 .or. lfbedfrm .or. lfsdu) then
             call wrsedmgrp(lundia    ,error     ,filename  ,itmapc    ,mmax      , &
                          & kmax      ,nmaxus    ,lsed      ,lsedtot   , &
                          & i(kfsmin) ,i(kfsmax) ,irequest  ,fds       ,iarrc     , &
                          & mf        ,ml        ,nf        ,nl        ,gdp       )
             if (error) goto 9999
          endif
          !
          if (lsedtot > 0) then
             if (initi>=4 .or. gdp%gdmorpar%moroutput%cumavg .or. irequest==REQUESTTYPE_DEFINE) then
                call wrsedmavg(lundia    ,error     ,filename  ,itmapc    ,mmax      , &
                             & nmaxus    ,lsed      ,lsedtot   ,irequest  ,fds       , &
                             & iarrc     ,mf        ,ml        ,nf        ,nl        , &
                             & gdp       )
                if (error) goto 9999
             endif
          endif
          !
          if (roller) then
             call wrrolm(lundia    ,error     ,filename  ,itmapc    ,nmax      , &
                       & mmax      ,nmaxus    ,r(ewave1) ,r(eroll1) ,r(qxkr)   , &
                       & r(qykr)   ,r(qxkw)   ,r(qykw)   ,r(fxw)    ,r(fyw)    , &
                       & r(wsu)    ,r(wsv)    ,r(guu)    ,r(gvv)    ,r(rbuff)  , &
                       & r(hrms)   ,irequest  ,fds       ,iarrc     ,mf        , &
                       & ml        ,nf        ,nl        ,gdp       )
             if (error) goto 9999
          endif
          !
          if (xbeach) then
             call wrxbm(lundia    ,error     ,filename  ,itmapc    ,nmax      , &
                      & mmax      ,nmaxus    ,r(fxw)    ,r(fyw)    , &
                      & r(wsu)    ,r(wsv)    ,r(guu)    ,r(gvv)    ,r(rbuff)  , &
                      & r(hrms)   ,irequest  ,fds       ,iarrc     ,mf        , &
                      & ml        ,nf        ,nl        ,gdp       )
             if (error) goto 9999
          endif
       endif
       !
       if (irequest == REQUESTTYPE_DEFINE) then
          !
          ! upon first request: define all groups, dimensions, and elements on file
          !
          call defnewgrp(fds, ifile, gdp, filename)
          !
          if (filetype == FTYPE_NETCDF) then
             !
             ierror = nf90_enddef(fds); call nc_check_err(lundia, ierror, "enddef", filename)
             if (ierror/=0) goto 9999
             !
          endif
       endif
    enddo
    !
9999 continue
    if (ierror/= 0) error = .true.
    if (inode == master) then
       if (filetype == FTYPE_NETCDF) then
          ierror = nf90_sync(fds); call nc_check_err(lundia, ierror, "sync file", filename)
          ierror = nf90_close(fds); call nc_check_err(lundia, ierror, "closing file", filename)
       elseif (filetype == FTYPE_NEFIS) then
          ierror = clsnef(fds)
          !
          ! wridoc needs to be called AFTER defnewgrp because this routine will
          ! add the map-version not just to the file, but also to the group
          ! administration in memory and calling defnewgrp afterwards would thus
          ! duplicate definition of the map-version group.
          !
          if (first .and. .not.wrifou) then
             ftype = 'map'
             call wridoc(error     ,trifil    ,ftype     ,simdat    ,runtxt    , &
                       & .false.   ,part_nr   ,gdp       )
             !
             if (.not. parll .and. .not. error) then
                call wrplot(filename  ,lundia    ,error     ,mmax      ,nmax      , &
                          & nmaxus    ,i(kcs)    ,i(ibuff)  ,r(xz)     ,r(yz)     , &
                          & sferic    ,gdp       )
             endif
          endif
       endif
       if (ierror/= 0) error = .true.
    endif
    first = .false.
    !
    if (allocated(ipartition)) deallocate(ipartition, stat=istat)
    if (.not.mergemap) then
        inode = bck_inode
        parll = bck_parll
        nproc = bck_nproc
        gdp%gdparall => bck_gdparall
    endif
end subroutine wrm_main

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!=======================================================================
!
!   S E R I A L   R A N D O M   N U M B E R   G E N E R A T O R
!
!=======================================================================
! Here's how to use these functions:
!       Initialization: call Rans( N, StartVal, Seeds )
! This returns an array of N seeds;
! the returned seeds have values that partition the basic
! ranf cycle of 2**46 pseudorandom reals in [0,1] into independent sub-
! sequences. The second argument, StartVal, can be a cycle-starting seed
! for the full-period generator; if it is zero, a special seed will be
! used that produces statistically desirable behavior.
!       Use: call Ranf( Seed, RandNum )
! This returns a pseudorandom real and a seed to be passed back in the
! next invocation. The returned seed carries the state normally hidden
! in imperative generators.
!=======================================================================
module random
  integer,parameter :: IRandNumSize = 4, IBinarySize = 48
  integer,parameter :: Mod4096DigitSize = 12
  integer,parameter :: IZero = 0
  integer,parameter :: NPoissonLimit = 10

  integer,dimension( IRandNumSize ) :: &
       & Multiplier = (/373, 3707, 1442, 647/)
  integer,dimension( IRandNumSize ) :: &
       & DefaultSeed = (/ 3281, 4041, 595, 2376 /)
  real(kind=8),dimension( IRandNumSize ) :: &
       & Divisor = (/281474976710656.0,68719476736.0,16777216.0,4096.0/)

  integer :: IGauss = 0
  real(kind=8) :: GaussBak = 0.0d0

contains
  !=======================================================================
  subroutine ranf( Seed, RandNum )
    implicit none
    integer,dimension( IRandNumSize ) :: Seed
    real(kind=8) :: RandNum

    integer,dimension( IRandNumSize ) :: OutSeed

    RandNum = float( Seed( 4 ) ) / Divisor( 4 ) + &
         &    float( Seed( 3 ) ) / Divisor( 3 ) + &
         &    float( Seed( 2 ) ) / Divisor( 2 ) + &
         &    float( Seed( 1 ) ) / Divisor( 1 )

    call ranfmodmult( Multiplier, Seed, OutSeed )

    Seed = OutSeed

    return
  end subroutine ranf
  !=======================================================================

  !=======================================================================
  subroutine poissdev( Seed, AverNum , PoissNum)
    implicit none
    integer,dimension( IRandNumSize ) :: Seed
    real(kind=8) :: AverNum
    integer :: PoissNum

    real(kind=8) :: Norm, Repar, Proba
    real(kind=8) :: RandNum, GaussNum

    if(AverNum <= DBLE(NPoissonLimit)) then
       Norm=exp(-AverNum)
       Repar=1.0d0
       PoissNum=0
       Proba=1.0d0
       call ranf(Seed,RandNum)
       do while(Repar*Norm <= RandNum .and. PoissNum <= 10*NPoissonLimit )
          PoissNum=PoissNum+1
          Proba=Proba*AverNum/PoissNum
          Repar=Repar+Proba
       end do
    else
       call gaussdev(Seed,GaussNum)
       GaussNum=GaussNum*sqrt(AverNum)-0.5+AverNum
       if(GaussNum<=0.0d0)GaussNum=0.0d0
       PoissNum=nint(GaussNum)
    endif

    return
  end subroutine poissdev
  !=======================================================================

  !=======================================================================
  subroutine gaussdev( Seed, GaussNum )
    implicit none
    integer,dimension( IRandNumSize ) :: Seed
    real(kind=8) :: GaussNum

    real(kind=8) :: fac,rsq,v1,v2

    if (IGauss.eq.IZero) then
       rsq=0.0d0
       do while (rsq.ge.1.0d0.or.rsq.le.0.0d0)
          call ranf(Seed,v1)
          call ranf(Seed,v2)
          v1=2.0d0*v1-1.0d0
          v2=2.0d0*v2-1.0d0
          rsq=v1**2+v2**2
       end do
       fac=sqrt(-2.0d0*log(rsq)/rsq)
       GaussBak=v1*fac
       GaussNum=v2*fac
       IGauss=1
    else
       GaussNum=GaussBak
       IGauss=0
    endif
    return
  END subroutine gaussdev
  !=======================================================================

  !=======================================================================
  subroutine rans( N, StartVal, SeedArray )
    implicit none
    integer :: N
    integer :: StartVal
    integer,dimension( 1, IRandNumSize ) :: SeedArray

    integer,dimension( IRandNumSize ) :: atothek
    integer,dimension( IRandNumSize ) :: K
    integer,dimension( IRandNumSize ) :: InSeed
    integer,dimension( IRandNumSize ) :: OutSeed
    integer,dimension( IBinarySize )  :: KBinary

    integer :: I

    if( StartVal .eq. IZero ) then
       SeedArray( 1, 1 ) = DefaultSeed( 1 )
       SeedArray( 1, 2 ) = DefaultSeed( 2 )
       SeedArray( 1, 3 ) = DefaultSeed( 3 )
       SeedArray( 1, 4 ) = DefaultSeed( 4 )
    else
       SeedArray( 1, 1 ) = abs( StartVal )
       SeedArray( 1, 2 ) = IZero
       SeedArray( 1, 3 ) = IZero
       SeedArray( 1, 4 ) = IZero
    endif

    if( N .eq. 1 ) then
       atothek( 1 ) = Multiplier( 1 )
       atothek( 2 ) = Multiplier( 2 )
       atothek( 3 ) = Multiplier( 3 )
       atothek( 4 ) = Multiplier( 4 )
    else
       call ranfk( N, K )
       call ranfkbinary( K, KBinary )
       call ranfatok( Multiplier, KBinary, atothek )
       do I = 2, N
          InSeed( 1 ) = SeedArray( I-1, 1 )
          InSeed( 2 ) = SeedArray( I-1, 2 )
          InSeed( 3 ) = SeedArray( I-1, 3 )
          InSeed( 4 ) = SeedArray( I-1, 4 )
          call ranfmodmult( InSeed, atothek, OutSeed )
          SeedArray( I, 1 ) = OutSeed( 1 )
          SeedArray( I, 2 ) = OutSeed( 2 )
          SeedArray( I, 3 ) = OutSeed( 3 )
          SeedArray( I, 4 ) = OutSeed( 4 )
       end do
    endif

    return
  end subroutine rans
  !=======================================================================

  !=======================================================================
  subroutine ranfatok( a, Kbinary, atothek )
    implicit none
    integer,dimension( IRandNumSize ) :: a
    integer,dimension( IBinarySize )  :: KBinary
    integer,dimension( IRandNumSize ) :: atothek

    integer,dimension( IRandNumSize ) :: asubi

    integer :: I
    asubi( 1 ) = a( 1 )
    asubi( 2 ) = a( 2 )
    asubi( 3 ) = a( 3 )
    asubi( 4 ) = a( 4 )

    atothek( 1 ) = 1
    atothek( 2 ) = IZero
    atothek( 3 ) = IZero
    atothek( 4 ) = IZero

    do I = 1, 45
       if( KBinary( I ) .ne. IZero ) then
          call ranfmodmult( atothek, asubi, atothek )
       endif
       call ranfmodmult( asubi, asubi, asubi )
    end do

    return
  end subroutine ranfatok
  !=======================================================================

  !=======================================================================
  subroutine ranfk( N, K )
    implicit none
    integer :: N
    integer, dimension( IRandNumSize ) :: K

    integer :: nn, r4, r3, r2, q4, q3, q2, q1

    nn = N + iranfeven( N )

    q4 = 1024 / nn
    r4 = 1024 - (nn * q4)
    q3 = (r4 * 4096) / nn
    r3 = (r4 * 4096) - (nn * q3)
    q2 = (r3 * 4096) / nn
    r2 = (r3 * 4096) - (nn * q2)
    q1 = (r2 * 4096) / nn

    K( 1 ) = q1
    K( 2 ) = q2
    K( 3 ) = q3
    K( 4 ) = q4

    return
  end subroutine ranfk
  !=======================================================================

  !=======================================================================
  subroutine ranfkbinary( K, KBinary )
    implicit none
    integer, dimension( IRandNumSize ) :: K
    integer, dimension( IBinarySize ) ::  KBinary

    integer, dimension( Mod4096DigitSize ) :: Bits
    integer X, I, J

    do I = 1, 4
       X = K( I ) / 2
       Bits( 1 ) = iranfodd( K( I ) )
       do J = 2, Mod4096DigitSize
          Bits( J ) = iranfodd( X )
          X = X / 2
       end do
       do J = 1, Mod4096DigitSize
          KBinary( (I-1)*Mod4096DigitSize + J ) = Bits( J )
       end do
    end do

    return
  end subroutine ranfkbinary
  !=======================================================================

  !=======================================================================
  subroutine ranfmodmult( A, B, C )
    implicit none
    integer, dimension( IRandNumSize ) :: A, B, C

    integer :: j1, j2, j3, j4, k1, k2, k3, k4

    j1 = A( 1 ) * B( 1 )
    j2 = A( 1 ) * B( 2 ) + A( 2 ) * B( 1 )
    j3 = A( 1 ) * B( 3 ) + A( 2 ) * B( 2 ) + A( 3 ) * B( 1 )
    j4 = A( 1 ) * B( 4 ) + A( 2 ) * B( 3 ) + A( 3 ) * B( 2 ) + A( 4 ) * B( 1 )

    k1 = j1
    k2 = j2 + k1 / 4096
    k3 = j3 + k2 / 4096
    k4 = j4 + k3 / 4096

    C( 1 ) = mod( k1, 4096 )
    C( 2 ) = mod( k2, 4096 )
    C( 3 ) = mod( k3, 4096 )
    C( 4 ) = mod( k4, 4096 )

    return
  end subroutine ranfmodmult
  !=======================================================================

  !=======================================================================
  function iranfodd( N )
    integer N
    if( mod( N, 2 ) .eq. 0 ) then
       iranfodd = 0
    else
       iranfodd = 1
    endif

    return
  end function iranfodd
  !=======================================================================
  !=======================================================================
  function iranfeven( N )
    integer N
    if( mod( N, 2 ) .eq. 0 ) then
       iranfeven = 1
    else
       iranfeven = 0
    endif
    return
  end function iranfeven
  !=======================================================================

end module random

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
!
! Copyright (c) 1999, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!       Testing Gen_block with module - must run with  4 processors
!       data initialization or array parameter in module

        module foo
	integer gba(4)
	data gba /1,1,1,1/
	integer, parameter, dimension(4)::gbb=(/1,1,1,1/)
        integer a(4),b(4)
!hpf$ distribute a(gen_block(gba))
!hpf$ distribute b(gen_block(gbb))
	contains
	subroutine s(n)
	integer ca(4),cb(4)
!hpf$ distribute ca(gen_block(gba))
!hpf$ distribute cb(gen_block(gbb))
	ca = a*n
	cb = ca
	b = cb
	end subroutine s
	subroutine t(n)
	call tt(n)
	contains
	subroutine tt(n)
	integer da(4),db(4)
!hpf$ distribute da(gen_block(gba))
!hpf$ distribute db(gen_block(gbb))
	da = a*n
	db = da
	b = db
	end subroutine tt
	end subroutine t

        end module

        module foo2
        use foo
        integer ea(4),eb(4)
!hpf$ distribute ea(gen_block(gba))
!hpf$ distribute eb(gen_block(gbb))
        end module

        program p
        use foo2
	integer result(32),expect(32)
        data expect/1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,3,             &
     &		6,9,12,1,2,3,4,2,4,6,8/
	forall(i=1:4) a(i) = i
	forall(i=1:4) b(i) = i
	forall(i=1:4) ea(i) = i
	forall(i=1:4) eb(i) = i
        !print *, a,b
        !print *, ea,eb
	result(1:4) = a
	result(5:8) = b
	result(9:12) = ea
	result(13:16) = eb
	call s(3)
	!print *, a,b
	result(17:20) = a
	result(21:24) = b
	call t(2)
	!print *, a,b
	result(25:28) = a
	result(29:32) = b
	call check(result,expect,32)
        end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
      SUBROUTINE OREZ (A,NA)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(1)
      DO 1 K=1,NA
1     A(K)=0.0D+0
      RETURN
      END

max_stars_repo_name: aldorona/contur

max_issues_repo_name: aldorona/contur

max_forks_repo_name: aldorona/contur


text:
*
      SUBROUTINE linpol( x1, y1, x2, y2, x, y )
*--------------------------------------------------------------------
* Linearly interpolate between (x1,y1) and (x2,y2), the new value
* at x is y.
*
      REAL a, b, x1, y1, x2, y2, x, y
      a = ( y2-y1 ) / ( x2 -x1 )
      b = y1 - a*x1
      y = a*x + b
      RETURN
      END