text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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 registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      module dhmmra_mod
      contains
      subroutine dhmmra ( lunrep, l_decl, arrpoi, arrtyp, arrbyt,
     &                    arrlen, arrknd, arrdm1, arrdm2, arrdm3,
     &                    arrnam, itota, part )

!     Deltares Software Centre

!>\file Define the real arrays as part of the one overall array
!>    Sets the array pointers in the SYSA common
!>    block. Gives array space of the kind A(pointer)
!>    Declares memory through C-interface if asked
!>    (routine is also called by preprocessor)

!     Created             : June 1998 by Jan van Beek
!     Modified            : May  2010 by Leo Postma
!                           Adds a number of arrays with normal names through
!                           the Fortran allocate feature and the waqmem module

!     Files               : LUNREP - monitoring output file

!     Routines            : SRSTOP, stops execution (on error)

      use waqmem           ! module with the more recently added arrays
      use partition_arrays ! module for computing the pointers into the arrays

      implicit none

      include "omp_lib.h"

!     Parameters          :

!     kind     function         name        description

      integer      , intent(in   ) :: lunrep    ! logical unitnumber output file
      logical      , intent(in   ) :: l_decl    ! Declare memory y/n
!     integer      , intent(in   ) :: noarr     ! dimension of parameter arrays, contained in 'sysn.inc'
      integer      , intent(inout) :: arrpoi(:) ! Pointer in workarray/FMM reference pointer
      integer      , intent(inout) :: arrtyp(:) ! Array type ( INT=,REAL=,CHAR= ), see FMM/NEFIS
      integer      , intent(inout) :: arrbyt(:) ! Number of bytes per element, see FMM/NEFIS
      integer      , intent(inout) :: arrlen(:) ! Length off array
      integer      , intent(inout) :: arrknd(:) ! Kind of array 1=(NOVAR), 2=(NOVAR,NOSEG) or 3=(NOSEG,NOVAR)
      integer      , intent(inout) :: arrdm1(:) ! dimension 1
      integer      , intent(inout) :: arrdm2(:) ! dimension 2
      integer      , intent(inout) :: arrdm3(:) ! dimension 3 ( number of grids mostly )
      character(20), intent(inout) :: arrnam(:) ! Array name
      integer      , intent(inout) :: itota     ! Required array space
      type(memory_partition), intent(inout) :: part ! Private variables for MAKPTR

!     COMMON  /  SYSN   /   System characteristics

      INCLUDE 'sysn.inc'

!     COMMON  /  SYSI   /   Timer characteristics

      INCLUDE 'sysi.inc'

!     COMMON  /  SYSA   /   Pointers in real array workspace

      INCLUDE 'sysa.inc'

!     Local declarations

      integer         i_rar                             ! loop counter
      integer         nr_rar                            ! number of real arrays
      integer         nohor                             ! number of computational volumes in 1 layer
      integer         nsubs                             ! nr of substances for array space declaration
      logical         fluxco                            ! if .true. then flux correction
      logical         steady                            ! if .true. then steady state computation
      logical         iterat                            ! if .true. then iterative solution
      logical         delmat                            ! if .true. then direct Gauss solver
      logical         f_solv                            ! if .true. then GMRES Krilov solver
      logical         triadi                            ! if .true. then ADI like Delft3d-Flow
      logical         balans                            ! if .true. then balances to be computed
      character*20    namarr                            ! help variable for array name
      integer         iartyp                            ! help variable for array type
      integer         iarlen                            ! help variable for array length
      integer         ip                                ! help variable for array pointer
      integer         ip_rar(iasize)                    ! help array to fill the common block / SYSA /
      equivalence   ( ivol   , ip_rar(1) )              ! first entry equivalences with first entry common block
      integer         noth                              ! number of available thread for parallel processing
      integer         ierr                              ! error indicator
      integer         jstart                            ! lower limit Flow arrays method 19 and 20
      integer         nmmaxj                            ! upper limit Flow arrays method 19 and 20

      logical            :: lfound                      ! argument was found
      character(len=256) :: adummy                      ! dummy string
      integer            :: nothreadsarg                ! optional number of threads from delwaq2 commandline arguments
      real               :: rdummy                      ! dummy real
      integer            :: ierr2                       ! error code

      integer   iivol  / 1/, iiarea / 2/, iiflow / 3/, iileng / 4/, iidisp / 5/,
     &          iiconc / 6/, iimass / 7/, iiderv / 8/, iiboun / 9/, iibset /10/,
     &          iibsav /11/, iiwste /12/, iicons /13/, iiparm /14/, iifunc /15/,
     &          iisfun /16/, iidnew /17/, iidiff /18/, iivnew /19/, iivelo /20/,
     &          iiharm /21/, iifarr /22/, iimas2 /23/, iitimr /24/, iivol2 /25/,
     &
     &                       iismas /32/, iiploc /33/, iidefa /34/, iiflux /35/,
     &          iistoc /36/, iiflxd /37/, iiflxi /38/, iiriob /39/, iidspx /40/,
     &          iivelx /41/, iilocx /42/, iidsto /43/, iivsto /44/, iidmpq /45/,
     &          iidmps /46/, iitrra /47/, iinrsp /48/, iivoll /49/,
     &          iir1   /51/, iiqxk  /52/, iiqyk  /53/, iiqzk  /54/, iidifx /55/,
     &          iidify /56/, iidifz /57/, iivola /58/, iivolb /59/, iiguv  /60/,
     &          iigvu  /61/,              iiaak  /63/, iibbk  /64/, iicck  /65/,
     &          iibd3x /66/, iibddx /67/, iibdx  /68/, iibu3x /69/, iibuux /70/,
     &          iibux  /71/, iiwrk1 /72/, iiwrk2 /73/, iiaakl /74/, iibbkl /75/,
     &          iicckl /76/, iiddkl /77/, iiwdmp /78/

!     How many threads ?

!     The '-threads [N]' argument for delwaq2 will turn on parallelism, and override
!     any setting of the number of threads in the input file.
!     No value or zero for [N] will use the maximum number of available threads
      nothreadsarg = 0
      call getcom ( '-threads', 1, lfound, nothreadsarg, rdummy, adummy, ierr2)
      if (lfound) then
         nothrd = nothreadsarg
      else
         call getcom ( '-nothreads', 1, lfound, nothreadsarg, rdummy, adummy, ierr2)
         if (lfound) then
            nothrd = nothreadsarg
         end if
      end if

      if ( nothrd .gt. 0 ) call OMP_SET_NUM_THREADS( nothrd )
      noth = OMP_GET_MAX_THREADS()
      write ( lunrep , 2020 ) noth
      if ( l_decl ) write (    6   , 2030 ) noth

!     Some logicals

      fluxco = intsrt .eq.  5 .or. intsrt .eq. 12 .or. intsrt .eq. 14 .or.
     &         intsrt .eq. 24
      steady = intsrt .eq.  6 .or. intsrt .eq.  7 .or. intsrt .eq.  8 .or.
     &         intsrt .eq.  9 .or. intsrt .eq. 17 .or. intsrt .eq. 18
      iterat = intsrt .eq.  8 .or. intsrt .eq.  9
      delmat = intsrt .eq.  6 .or. intsrt .eq.  7 .or. intsrt .eq. 10
      f_solv = intsrt .eq. 15 .or. intsrt .eq. 16 .or. intsrt .eq. 17 .or.
     &         intsrt .eq. 18 .or. intsrt .eq. 21 .or. intsrt .eq. 22
      triadi = intsrt .eq. 19 .or. intsrt .eq. 20
      balans = btest(intopt,3)

!     Set defaults, no name no length

      nr_rar = iasize                   ! total number of arrays
      do i_rar = 1 , nr_rar
         arrnam(i_rar) = ' '
         arrtyp(i_rar) = rtyp
         arrbyt(i_rar) = 4
         arrknd(i_rar) = 0
         arrdm1(i_rar) = 0
         arrdm2(i_rar) = 0
         arrdm3(i_rar) = 0
         arrlen(i_rar) = 0
      enddo

      arrnam(iivol ) = 'VOLUME'
      arrknd(iivol ) = 2
      arrdm1(iivol ) = 1
      arrdm2(iivol ) = noseg+nseg2
      arrdm3(iivol ) = nogrid

      arrnam(iiarea) = 'AREA  '
      arrknd(iiarea) = 2
      arrdm1(iiarea) = 1
      arrdm2(iiarea) = noq+noq4
      arrdm3(iiarea) = 1

      arrnam(iiflow) = 'FLOW  '
      arrknd(iiflow) = 2
      arrdm1(iiflow) = 1
      arrdm2(iiflow) = noq+noq4
      arrdm3(iiflow) = 1

      arrnam(iileng) = 'LENG  '
      if ( ilflag .eq. 0 ) then
         arrknd(iileng) = 1
         arrdm1(iileng) = 3
         arrdm2(iileng) = 1
      else
         arrknd(iileng) = 2
         arrdm1(iileng) = 2
         arrdm2(iileng) = noq+noq4
      endif
      arrdm3(iileng) = 1

      arrnam(iidisp) = 'DISP  '
      arrknd(iidisp) = 1
      arrdm1(iidisp) = 3
      arrdm2(iidisp) = 1
      arrdm3(iidisp) = 1

      arrnam(iiconc) = 'CONC  '
      arrknd(iiconc) = 2
      if ( steady .and. .not. iterat ) then
         arrdm1(iiconc) = notot
         arrdm2(iiconc) = noseg+nseg2
         arrdm3(iiconc) = nogrid
         nsubs = notot
      else
         arrdm1(iiconc) = notot
         arrdm2(iiconc) = noseg+nseg2
         arrdm3(iiconc) = nogrid
         nsubs = nosys
      endif

      arrnam(iimass) = 'MASS  '
      arrknd(iimass) = 2
      if ( f_solv ) then
         arrdm1(iimass) = notot
         arrdm2(iimass) = noseg+nseg2+nobnd
      else
         arrdm1(iimass) = notot
         arrdm2(iimass) = noseg+nseg2
      endif
      arrdm3(iimass) = nogrid

      arrnam(iiderv) = 'DERIV '
      arrknd(iiderv) = 2
      if ( f_solv .and. steady ) then
         arrdm1(iiderv) = notot
         arrdm2(iiderv) = noseg+nseg2+nobnd
      else
         arrdm1(iiderv) = notot
         arrdm2(iiderv) = noseg+nseg2
      endif
      arrdm3(iiderv) = nogrid

      arrnam(iiboun) = 'BOUND '
      arrknd(iiboun) = 2
      arrdm1(iiboun) = nsubs
      arrdm2(iiboun) = nobnd
      arrdm3(iiboun) = 1

      arrnam(iibset) = 'BSET  '
      arrknd(iibset) = 2
      arrdm1(iibset) = nsubs
      arrdm2(iibset) = nobnd
      arrdm3(iibset) = 1

      arrnam(iibsav) = 'BSAVE '
      arrknd(iibsav) = 2
      arrdm1(iibsav) = nsubs
      arrdm2(iibsav) = nobnd
      arrdm3(iibsav) = 1

      arrnam(iiwste) = 'WASTE '
      arrknd(iiwste) = 2
      arrdm1(iiwste) = notot+2
      arrdm2(iiwste) = nowst
      arrdm3(iiwste) = 1

      arrnam(iicons) = 'CONS  '
      arrknd(iicons) = 1
      arrdm1(iicons) = nocons
      arrdm2(iicons) = 1
      arrdm3(iicons) = 1

      arrnam(iiparm) = 'PARAM '
      arrknd(iiparm) = 2
      arrdm1(iiparm) = nopa
      arrdm2(iiparm) = noseg+nseg2
      arrdm3(iiparm) = nogrid

      arrnam(iifunc) = 'FUNC  '
      arrknd(iifunc) = 1
      arrdm1(iifunc) = nofun
      arrdm2(iifunc) = 1
      arrdm3(iifunc) = 1

      arrnam(iisfun) = 'SFUNC '
      arrknd(iisfun) = 3
      arrdm1(iisfun) = noseg+nseg2
      arrdm2(iisfun) = nosfun
      arrdm3(iisfun) = nogrid

      arrnam(iidnew) = 'DISPNW'
      arrknd(iidnew) = 2
      arrdm1(iidnew) = ndspn
      arrdm2(iidnew) = noq+noq4
      arrdm3(iidnew) = 1

      arrnam(iidiff) = 'DISPER'
      arrknd(iidiff) = 2
      arrdm1(iidiff) = nodisp
      arrdm2(iidiff) = noq+noq4
      arrdm3(iidiff) = 1

      arrnam(iivnew) = 'VELONW'
      arrknd(iivnew) = 2
      arrdm1(iivnew) = nveln
      arrdm2(iivnew) = noq+noq4
      arrdm3(iivnew) = 1

      arrnam(iivelo) = 'VELO  '
      arrknd(iivelo) = 2
      arrdm1(iivelo) = novelo
      arrdm2(iivelo) = noq+noq4
      arrdm3(iivelo) = 1

      arrnam(iiharm) = 'HARMAT'
      arrknd(iiharm) = 1
      arrdm1(iiharm) = nharms
      arrdm2(iiharm) = 1
      arrdm3(iiharm) = 1

      arrnam(iifarr) = 'FARR  '
      arrknd(iifarr) = 1
      arrdm1(iifarr) = nlines
      arrdm2(iifarr) = 1
      arrdm3(iifarr) = 1

      arrnam(iimas2) = 'MASS2 '
      arrknd(iimas2) = 2
      arrdm1(iimas2) = notot
      arrdm2(iimas2) = 5
      arrdm3(iimas2) = 1

      arrnam(iitimr) = 'TIMER '
      if ( fluxco ) then
         arrknd(iitimr) = 2
         arrdm1(iitimr) = notot
         arrdm2(iitimr) = noseg+nseg2
         arrdm3(iitimr) = 1
      elseif ( delmat ) then
         arrknd(iitimr) = 3
         arrdm1(iitimr) = noseg+nseg2
         arrdm2(iitimr) = jtrack*2+1
         arrdm3(iitimr) = 1
      elseif ( iterat ) then
         arrknd(iitimr) = 2
         arrdm1(iitimr) = notot
         arrdm2(iitimr) = noseg+nseg2
         arrdm3(iitimr) = 1
      elseif ( f_solv ) then
         arrknd(iitimr) = 1
         arrdm1(iitimr) = nomat
         arrdm2(iitimr) = 1
         arrdm3(iitimr) = 1
      endif

      arrnam(iivol2) = 'VOL2  '
      if ( f_solv ) then
         arrknd(iivol2) = 3
         arrdm1(iivol2) = noseg+nseg2 + nobnd
         arrdm2(iivol2) = 1
         arrdm3(iivol2) = 1
      else
         arrknd(iivol2) = 3
         arrdm1(iivol2) = noseg+nseg2
         arrdm2(iivol2) = 1
         arrdm3(iivol2) = 1
      endif

      arrnam(iismas) = 'ASMASS'
      if ( balans ) then
         arrknd(iismas) = 4
         arrdm1(iismas) = notot
         arrdm2(iismas) = ndmpar
         arrdm3(iismas) = 6
      endif

      arrnam(iiploc) = 'LOCAL '
      arrknd(iiploc) = 2
      arrdm1(iiploc) = noloc
      arrdm2(iiploc) = noseg+nseg2
      arrdm3(iiploc) = nogrid

      arrnam(iidefa) = 'DEFAUL'
      arrknd(iidefa) = 1
      arrdm1(iidefa) = nodef
      arrdm2(iidefa) = 1
      arrdm3(iidefa) = 1

      arrnam(iiflux) = 'FLUX  '
      arrknd(iiflux) = 2
      arrdm1(iiflux) = nflux
      arrdm2(iiflux) = noseg+nseg2
      arrdm3(iiflux) = nogrid

      arrnam(iistoc) = 'STOCHI'
      arrknd(iistoc) = 4
      arrdm1(iistoc) = notot
      arrdm2(iistoc) = nflux
      arrdm3(iistoc) = 1

      arrnam(iiflxd) = 'FLXDMP'
      if ( balans ) then
         arrknd(iiflxd) = 3
         arrdm1(iiflxd) = ndmps
         arrdm2(iiflxd) = nflux
         arrdm3(iiflxd) = 1
      endif

      arrnam(iiflxi) = 'FLXINT'
      if ( balans ) then
         arrknd(iiflxi) = 3
         arrdm1(iiflxi) = ndmpar
         arrdm2(iiflxi) = nflux
         arrdm3(iiflxi) = 1
      endif

      arrnam(iiriob) = 'RIOBUF'
      arrknd(iiriob) = 1
      arrdm1(iiriob) = nbufmx
      arrdm2(iiriob) = 1
      arrdm3(iiriob) = 1

      arrnam(iidspx) = 'DISPX '
      arrknd(iidspx) = 2
      arrdm1(iidspx) = ndspx
      arrdm2(iidspx) = noq+noq4
      arrdm3(iidspx) = 1

      arrnam(iivelx) = 'VELX  '
      arrknd(iivelx) = 2
      arrdm1(iivelx) = nvelx
      arrdm2(iivelx) = noq+noq4
      arrdm3(iivelx) = 1

      arrnam(iilocx) = 'VLOCX '
      arrknd(iilocx) = 2
      arrdm1(iilocx) = nlocx
      arrdm2(iilocx) = noq+noq4
      arrdm3(iilocx) = 1

      arrnam(iidsto) = 'DSTO  '
      arrknd(iidsto) = 4
      arrdm1(iidsto) = nosys
      arrdm2(iidsto) = ndspx
      arrdm3(iidsto) = 1

      arrnam(iivsto) = 'VSTO  '
      arrknd(iivsto) = 4
      arrdm1(iivsto) = nosys
      arrdm2(iivsto) = nvelx
      arrdm3(iivsto) = 1

      arrnam(iidmpq) = 'DMPQ  '
      arrknd(iidmpq) = 4
      arrdm1(iidmpq) = nosys
      arrdm2(iidmpq) = ndmpq
      arrdm3(iidmpq) = 2

      arrnam(iidmps) = 'DMPS  '
      arrknd(iidmps) = 4
      arrdm1(iidmps) = notot
      arrdm2(iidmps) = ndmps
      arrdm3(iidmps) = 3

      arrnam(iitrra) = 'TRRAAI'
      arrknd(iitrra) = 3
      arrdm1(iitrra) = nosys
      arrdm2(iitrra) = noraai
      arrdm3(iitrra) = 1

      arrnam(iinrsp) = 'INWRSP'
      arrknd(iinrsp) = 1
      arrdm1(iinrsp) = newrsp
      arrdm2(iinrsp) = 1
      arrdm3(iinrsp) = 1

      arrnam(iivoll) = 'VOLUML'
      arrknd(iivoll) = 3
      arrdm1(iivoll) = noseg+nseg2
      arrdm2(iivoll) = 1
      arrdm3(iivoll) = 1

!     The next array's only for TRIADI solvers

      if ( triadi ) then
         nohor = nmax * (mmax+4)

         arrnam(iir1  ) = 'R1    '
         arrknd(iir1  ) = 4
         arrdm1(iir1  ) = nohor*kmax
         arrdm2(iir1  ) = notot
         arrdm3(iir1  ) = 1

         arrnam(iiqxk ) = 'QXK   '
         arrknd(iiqxk ) = 3
         arrdm1(iiqxk ) = nohor*kmax
         arrdm2(iiqxk ) = 1
         arrdm3(iiqxk ) = 1

         arrnam(iiqyk ) = 'QYK   '
         arrknd(iiqyk ) = 3
         arrdm1(iiqyk ) = nohor*kmax
         arrdm2(iiqyk ) = 1
         arrdm3(iiqyk ) = 1

         arrnam(iiqzk ) = 'QZK   '
         arrknd(iiqzk ) = 3
         arrdm1(iiqzk ) = nohor*(kmax+1)
         arrdm2(iiqzk ) = 1
         arrdm3(iiqzk ) = 1

         arrnam(iidifx) = 'DIFX  '
         arrknd(iidifx) = 3
         arrdm1(iidifx) = nohor*kmax
         arrdm2(iidifx) = 1
         arrdm3(iidifx) = 1

         arrnam(iidify) = 'DIFY  '
         arrknd(iidify) = 3
         arrdm1(iidify) = nohor*kmax
         arrdm2(iidify) = 1
         arrdm3(iidify) = 1

         arrnam(iidifz) = 'DIFZ  '
         arrknd(iidifz) = 3
         arrdm1(iidifz) = nohor*(kmax+1)
         arrdm2(iidifz) = 1
         arrdm3(iidifz) = 1

         arrnam(iivola) = 'VOL0  '
         arrknd(iivola) = 3
         arrdm1(iivola) = nohor*kmax
         arrdm2(iivola) = 1
         arrdm3(iivola) = 1

         arrnam(iivolb) = 'VOL1  '
         arrknd(iivolb) = 3
         arrdm1(iivolb) = nohor*kmax
         arrdm2(iivolb) = 1
         arrdm3(iivolb) = 1

         arrnam(iiguv ) = 'GUV   '
         arrknd(iiguv ) = 3
         arrdm1(iiguv ) = nohor
         arrdm2(iiguv ) = 1
         arrdm3(iiguv ) = 1

         arrnam(iigvu ) = 'GVU   '
         arrknd(iigvu ) = 3
         arrdm1(iigvu ) = nohor
         arrdm2(iigvu ) = 1
         arrdm3(iigvu ) = 1

         arrnam(iiaak ) = 'AAK   '
         arrknd(iiaak ) = 3
         arrdm1(iiaak ) = nohor*kmax
         arrdm2(iiaak ) = 1
         arrdm3(iiaak ) = 1

         arrnam(iibbk ) = 'BBK   '
         arrknd(iibbk ) = 3
         arrdm1(iibbk ) = nohor*kmax
         arrdm2(iibbk ) = 1
         arrdm3(iibbk ) = 1

         arrnam(iicck ) = 'CCK   '
         arrknd(iicck ) = 3
         arrdm1(iicck ) = nohor*kmax
         arrdm2(iicck ) = 1
         arrdm3(iicck ) = 1

         arrnam(iibd3x) = 'BD3X  '
         arrknd(iibd3x) = 3
         arrdm1(iibd3x) = nohor*kmax
         arrdm2(iibd3x) = 1
         arrdm3(iibd3x) = 1

         arrnam(iibddx) = 'BDDX  '
         arrknd(iibddx) = 3
         arrdm1(iibddx) = nohor*kmax
         arrdm2(iibddx) = 1
         arrdm3(iibddx) = 1

         arrnam(iibdx ) = 'BDX   '
         arrknd(iibdx ) = 3
         arrdm1(iibdx ) = nohor*kmax
         arrdm2(iibdx ) = 1
         arrdm3(iibdx ) = 1

         arrnam(iibu3x) = 'BU3X  '
         arrknd(iibu3x) = 3
         arrdm1(iibu3x) = nohor*kmax
         arrdm2(iibu3x) = 1
         arrdm3(iibu3x) = 1

         arrnam(iibuux) = 'BUUX  '
         arrknd(iibuux) = 3
         arrdm1(iibuux) = nohor*kmax
         arrdm2(iibuux) = 1
         arrdm3(iibuux) = 1

         arrnam(iibux ) = 'BUX   '
         arrknd(iibux ) = 3
         arrdm1(iibux ) = nohor*kmax
         arrdm2(iibux ) = 1
         arrdm3(iibux ) = 1

         arrnam(iiwrk1) = 'WRK1  '
         arrknd(iiwrk1) = 3
         arrdm1(iiwrk1) = nohor*kmax
         arrdm2(iiwrk1) = 1
         arrdm3(iiwrk1) = 1

         arrnam(iiwrk2) = 'WRK2  '
         arrknd(iiwrk2) = 3
         arrdm1(iiwrk2) = nohor*kmax
         arrdm2(iiwrk2) = 1
         arrdm3(iiwrk2) = 1

         arrnam(iiaakl) = 'AAKL  '
         arrknd(iiaakl) = 4
         arrdm1(iiaakl) = nohor*kmax
         arrdm2(iiaakl) = nosys
         arrdm3(iiaakl) = 1

         arrnam(iibbkl) = 'BBKL  '
         arrknd(iibbkl) = 4
         arrdm1(iibbkl) = nohor*kmax
         arrdm2(iibbkl) = nosys
         arrdm3(iibbkl) = 1

         arrnam(iicckl) = 'CCKL  '
         arrknd(iicckl) = 4
         arrdm1(iicckl) = nohor*kmax
         arrdm2(iicckl) = nosys
         arrdm3(iicckl) = 1

         arrnam(iiddkl) = 'DDKL  '
         arrknd(iiddkl) = 4
         arrdm1(iiddkl) = nohor*kmax
         arrdm2(iiddkl) = nosys
         arrdm3(iiddkl) = 1
      endif

      arrnam(iiwdmp) = 'WSTDMP'
      arrknd(iiwdmp) = 4
      arrdm1(iiwdmp) = notot
      arrdm2(iiwdmp) = nowst
      arrdm3(iiwdmp) = 2

!     the total array length

      if ( .not. l_decl ) then
         open  ( 328, file='memory_map.out' )
         write ( 328, '(/a/a/)' ) "  ==> REAL arrays in 4-byte words <==",
     &                            "  nr array name            array size"
      endif

      itota = 0
      do i_rar = 1 , nr_rar
         arrlen(i_rar) = arrdm1(i_rar)*arrdm2(i_rar)*arrdm3(i_rar)
         if ( .not. l_decl ) write ( 328, 2040 ) i_rar, arrnam(i_rar), arrlen(i_rar)
         itota = itota + arrlen(i_rar)
         if ( itota .lt. 0 ) then
            write(lunrep,2005)
            call srstop(1)
         endif
      enddo

!     Declare memory

      if ( l_decl ) then
         do i_rar = 1 , nr_rar
            iartyp = arrtyp(i_rar)
            iarlen = arrlen(i_rar)
            namarr = arrnam(i_rar)
            if ( iarlen .gt. 0 ) then
               ip = makptr(part, namarr,iartyp ,iarlen)
               if ( ip .le. 0 ) then
                  write(lunrep,2010) namarr
                  call srstop(1)
               endif
            else
               ip = 0
            endif

!           Add one extra because of the shift between rbuf(0) and a(1)

            ip = ip + 1
            ip_rar(i_rar) = ip
            arrpoi(i_rar) = ip
         enddo
      endif

!     Reset new disp and velo pointers if array's are the same

      if ( ndspn .eq. 0 ) then
         idnew = idiff
         arrknd(iidnew) = arrknd(iidiff)
         arrdm1(iidnew) = arrdm1(iidiff)
         arrdm2(iidnew) = arrdm2(iidiff)
         arrdm3(iidnew) = arrdm3(iidiff)
         arrlen(iidnew) = arrlen(iidiff)
         arrpoi(iidnew) = arrpoi(iidiff)
      endif
      if ( nveln .eq. 0 ) then
         ivnew = ivelo
         arrknd(iivnew) = arrknd(iivelo)
         arrdm1(iivnew) = arrdm1(iivelo)
         arrdm2(iivnew) = arrdm2(iivelo)
         arrdm3(iivnew) = arrdm3(iivelo)
         arrlen(iivnew) = arrlen(iivelo)
         arrpoi(iivnew) = arrpoi(iivelo)
      endif

!     New array declarations

!     (Make sure there is no allocated memory from a possible previous run.
!     This is a problem if DELWAQ is used as a library)

      if ( l_decl ) then
          if ( allocated( surface )  ) deallocate( surface )
          if ( allocated( rhs )      ) deallocate( rhs )
          if ( allocated( arhs )     ) deallocate( arhs )
          if ( allocated( adiag )    ) deallocate( adiag )
          if ( allocated( acodia )   ) deallocate( acodia )
          if ( allocated( bcodia )   ) deallocate( bcodia )
          if ( allocated( r11 )      ) deallocate( r11 )
          if ( allocated( dfluxx )   ) deallocate( dfluxx )
          if ( allocated( dfluxy )   ) deallocate( dfluxy )
          if ( allocated( s0 )       ) deallocate( s0 )
          if ( allocated( s1 )       ) deallocate( s1 )
          if ( allocated( dps )      ) deallocate( dps )
          if ( allocated( gsqs )     ) deallocate( gsqs )
          if ( allocated( sigdif )   ) deallocate( sigdif )
          if ( allocated( sigmol )   ) deallocate( sigmol )
          if ( allocated( guu    )   ) deallocate( guu    )
          if ( allocated( gvv    )   ) deallocate( gvv    )
          if ( allocated( hu     )   ) deallocate( hu     )
          if ( allocated( hv     )   ) deallocate( hv     )
          if ( allocated( thick  )   ) deallocate( thick  )
          if ( allocated( sig    )   ) deallocate( sig    )
          if ( allocated( dicuv  )   ) deallocate( dicuv  )
          if ( allocated( dicww  )   ) deallocate( dicww  )
          if ( allocated( sour   )   ) deallocate( sour   )
          if ( allocated( sink   )   ) deallocate( sink   )
          if ( allocated( areau  )   ) deallocate( areau  )
          if ( allocated( areav  )   ) deallocate( areav  )
          if ( allocated( rscale )   ) deallocate( rscale )
          if ( allocated( cell_x )   ) deallocate( cell_x )
          if ( allocated( cell_y )   ) deallocate( cell_y )
          if ( allocated( mixlen )   ) deallocate( mixlen )
          if ( allocated( gm_rhs )   ) deallocate( gm_rhs )
          if ( allocated( gm_sol )   ) deallocate( gm_sol )
          if ( allocated( gm_work )  ) deallocate( gm_work )
          if ( allocated( gm_hess )  ) deallocate( gm_hess )
          if ( allocated( gm_amat )  ) deallocate( gm_amat )
          if ( allocated( gm_diag )  ) deallocate( gm_diag )
          if ( allocated( gm_diac )  ) deallocate( gm_diac )
          if ( allocated( gm_trid )  ) deallocate( gm_trid )
          if ( allocated( flowtot )  ) deallocate( flowtot )
          if ( allocated( disptot )  ) deallocate( disptot )
          if ( allocated( theta   )  ) deallocate( theta   )
          if ( allocated( thetaseg ) ) deallocate( thetaseg )
          if ( allocated( flux )     ) deallocate( flux )
          if ( allocated( lim )      ) deallocate( lim )
          if ( allocated( maxi )     ) deallocate( maxi )
          if ( allocated( mini )     ) deallocate( mini )
          if ( allocated( l1 )       ) deallocate( l1 )
          if ( allocated( l2 )       ) deallocate( l2 )
          if ( allocated( m1 )       ) deallocate( m1 )
          if ( allocated( m2 )       ) deallocate( m2 )
          if ( allocated( n1 )       ) deallocate( n1 )
          if ( allocated( n2 )       ) deallocate( n2 )
      endif

      ierr = 0
      itota  = itota  +   noseg+nseg2
      nr_rar = nr_rar + 1
      if ( l_decl ) allocate ( surface   ( noseg+nseg2 ), stat=ierr )
      if ( ierr .ne. 0 ) then ; write(lunrep,2010) "surface             " ; call srstop(1) ; endif
      if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "surface             ", noseg+nseg2

      itota  = itota  +   noseg+nseg2
      nr_rar = nr_rar + 1
      if ( l_decl ) allocate ( wdrawal   ( noseg+nseg2 ), stat=ierr )
      if ( ierr .ne. 0 ) then ; write(lunrep,2010) "wdrawal             " ; call srstop(1) ; endif
      if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "wdrawal             ", noseg+nseg2

      if ( delmat ) then
         itota  = itota  +   nosys*noseg
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( rhs   ( nosys,noseg ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "rhs                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "rhs                 ", nosys*noseg
      endif

      if ( intsrt .eq. 11 .or. intsrt .eq. 12 .or.
     &     intsrt .eq. 13 .or. intsrt .eq. 14 .or.
     &     intsrt .eq. 24                          ) then
         itota  = itota  +  notot*(noseg+nseg2)*2                        ! arhs
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( arhs    ( notot,noseg+nseg2), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "arhs                " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "arhs                ", notot*(noseg+nseg2)*2

         itota  = itota  +  notot*(noseg+nseg2)*2                        ! adiag
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( adiag   ( notot,noseg+nseg2), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "adiag               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "adiag               ", notot*(noseg+nseg2)*2

         itota  = itota  +  notot*max((noq3+noq4),1)*2           ! acodia
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( acodia  ( notot,max(noq3+noq4,1)), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "acodia              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "acodia              ", notot*max((noq3+noq4),1)*2

         itota  = itota  +  notot*max((noq3+noq4),1)*2           ! bcodia
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( bcodia  ( notot,max(noq3+noq4,1)), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "bcodia              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "bcodia              ", notot*max((noq3+noq4),1)*2
      endif

      if ( triadi ) then
         jstart = 1 - 2 * nmax
         nmmaxj = ( 2 + mmax ) * nmax

         itota  = itota  +  (nmmaxj-jstart+1)*kmax*nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( r11   ( jstart:nmmaxj , kmax , nosys ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "r11                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "r11                 ", (nmmaxj-jstart+1)*kmax*nosys

         itota  = itota  +  (nmmaxj-jstart+1)*kmax*nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dfluxx( jstart:nmmaxj , kmax , nosys ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dfluxx              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dfluxx              ", (nmmaxj-jstart+1)*kmax*nosys

         itota  = itota  +  (nmmaxj-jstart+1)*kmax*nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dfluxy( jstart:nmmaxj , kmax , nosys ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dfluxy              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dfluxy              ", (nmmaxj-jstart+1)*kmax*nosys

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( s0    ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "s0                  " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "s0                  ", (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( s1    ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "s1                  " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "s1                  ", (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dps   ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dps                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dps                 ", (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gsqs  ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gsqs                " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gsqs                ", (nmmaxj-jstart+1)

         itota  = itota  +  nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( sigdif( nosys         )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "sigdif              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "sigdif              ",  nosys

         itota  = itota  +  nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( sigmol( nosys         )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "sigmol              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "sigmol              ",  nosys

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( guu   ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "guu                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "guu                 ",  (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gvv   ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gvv                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gvv                 ",  (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( hu    ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "hu                  " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "hu                  ",  (nmmaxj-jstart+1)

         itota  = itota  +  (nmmaxj-jstart+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( hv    ( jstart:nmmaxj )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "hv                  " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "hv                  ",  (nmmaxj-jstart+1)

         itota  = itota  +  nolay
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( thick ( nolay         )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "thick               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "thick               ",  nolay

         itota  = itota  +  nolay
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( sig   ( nolay         )               , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "sig                 " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "sig                 ",  nolay

         itota  = itota  +  (nmmaxj-jstart+1)*kmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dicuv ( jstart:nmmaxj , kmax  )       , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dicuv               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dicuv               ", (nmmaxj-jstart+1)*kmax

         itota  = itota  +  (nmmaxj-jstart+1)*(kmax+1)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dicww ( jstart:nmmaxj , 0:kmax  )     , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dicww               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dicww               ", (nmmaxj-jstart+1)*(kmax+1)

         itota  = itota  +  (nmmaxj-jstart+1)*kmax*nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( sour  ( jstart:nmmaxj, kmax, nosys )  , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "sour                " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "sour                ", (nmmaxj-jstart+1)*kmax*nosys

         itota  = itota  +  (nmmaxj-jstart+1)*kmax*nosys
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( sink  ( jstart:nmmaxj, kmax, nosys )  , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "sink                " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "sink                ", (nmmaxj-jstart+1)*kmax*nosys

         itota  = itota  +  (nmmaxj-jstart+1)*kmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( areau ( jstart:nmmaxj, kmax )         , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "areau               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "areau               ", (nmmaxj-jstart+1)*kmax

         itota  = itota  +  (nmmaxj-jstart+1)*kmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( areav ( jstart:nmmaxj, kmax )         , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "areav               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "areav               ", (nmmaxj-jstart+1)*kmax

         itota  = itota  +  (nmmaxj-jstart+1)*kmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( rscale( jstart:nmmaxj, kmax )         , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "rscale              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "rscale              ", (nmmaxj-jstart+1)*kmax
      endif

      if ( nmax*mmax .gt. 0 ) then
         itota  = itota  +  nmax*mmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( cell_x( nmax, mmax )                  , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "cell_x              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "cell_x              ", nmax*mmax

         itota  = itota  +  nmax*mmax
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( cell_y( nmax, mmax )                  , stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "cell_y              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "cell_y              ", nmax*mmax
      endif

      if ( f_solv ) then
         itota  = itota  +  noq
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( mixlen (  noq                         ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "mixlen              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "mixlen              ",  noq

         itota  = itota  + (noseg+nobnd)          *noth*2           ! gm_rhs      real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_rhs (  noseg+nobnd           ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_rhs              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_rhs              ", (noseg+nobnd)          *noth*2

         itota  = itota  + (noseg+nobnd)          *noth*2           ! gm_sol      real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_sol (  noseg+nobnd           ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_sol              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_sol              ", (noseg+nobnd)          *noth*2

         itota  = itota  + (noseg+nobnd)*(novec+5)*noth*2           ! gm_work     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_work( (noseg+nobnd)*(novec+5),noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_work             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_work             ", (noseg+nobnd)*(novec+5)*noth*2

         itota  = itota  + (novec+1)*(novec+2)    *noth*2           ! gm_hess     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_hess( (  novec+1  )*(novec+2),noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_hess             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_hess             ", (novec+1)*(novec+2)    *noth*2

         itota  = itota  +  nomat                 *noth*2           ! gm_amat     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_amat(  nomat                 ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_amat             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_amat             ",  nomat                 *noth*2

         itota  = itota  + (noseg+nobnd)          *noth*2           ! gm_diag     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_diag(  noseg+nobnd           ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_diag             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_diag             ", (noseg+nobnd)          *noth*2

         itota  = itota  + (noseg+nobnd)          *noth*2           ! gm_diac     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_diac(  noseg+nobnd           ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_diac             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_diac             ", (noseg+nobnd)          *noth*2

         itota  = itota  + 6*nolay                *noth*2           ! gm_trid     real(8)
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( gm_trid(  6*nolay               ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "gm_trid             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "gm_trid             ", 6*nolay                *noth*2

         itota  = itota  +  noq  *noth                           ! flowtot
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( flowtot ( noq  ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "flowtot             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "flowtot             ", noq  *noth

         itota  = itota  +  noq  *noth                           ! disptot
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( disptot ( noq  ,noth ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "disptot             " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "disptot             ", noq  *noth

         if ( intsrt .eq. 21 ) then
            itota  = itota  +  noq  *noth                           ! theta
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( theta   ( noq  ,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "theta               " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "theta               ", noq  *noth

            itota  = itota  +  noseg*noth                           ! thetaseg
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( thetaseg( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "thetaseg            " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "thetaseg            ", noseg*noth

            itota  = itota  +  noq  *noth                           ! flux
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( flux    ( noq  ,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "flux                " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "flux                ", noq  *noth

            itota  = itota  +  noq  *noth                           ! lim
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( lim     ( noq  ,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "lim                 " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "lim                 ", noq  *noth

            itota  = itota  +  noseg*noth                           ! maxi
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( maxi    ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "maxi                " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "maxi                ", noseg*noth

            itota  = itota  +  noseg*noth                           ! mini
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( mini    ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "mini                " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "mini                ", noseg*noth

            itota  = itota  +  noseg*noth                           ! l1
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( l1      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "l1                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "l1                  ", noseg*noth

            itota  = itota  +  noseg*noth                           ! l2
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( l2      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "l2                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "l2                  ", noseg*noth

            itota  = itota  +  noseg*noth                           ! m1
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( m1      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "m1                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "m1                  ", noseg*noth

            itota  = itota  +  noseg*noth                           ! m2
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( m2      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "m2                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "m2                  ", noseg*noth

            itota  = itota  +  noseg*noth                           ! n1
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( n1      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "n1                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "n1                  ", noseg*noth

            itota  = itota  +  noseg*noth                           ! n2
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( n2      ( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "n2                  " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "n2                  ", noseg*noth
         endif
         if ( intsrt .eq. 22 ) then
            itota  = itota  +  noq  *noth                           ! theta
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( theta   ( noq  ,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "theta               " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "theta               ", noq  *noth

            itota  = itota  +  noseg*noth                           ! thetaseg
            nr_rar = nr_rar + 1
            if ( l_decl ) allocate ( thetaseg( noseg,noth ), stat=ierr )
            if ( ierr .ne. 0 ) then ; write(lunrep,2010) "thetaseg            " ; call srstop(1) ; endif
            if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "thetaseg            ", noseg*noth
         endif
      endif
      if ( intsrt .eq. 24 ) then
         itota  = itota  +  3*noseg*2                            ! dwork
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dwork   ( 3, noseg   ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dwork               " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dwork               ", 3*noseg*2

         itota  = itota  +  noseg*2                              ! volint
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( volint  ( noseg      ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "volint              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "volint              ", noseg*2

         itota  = itota  +  notot*(noseg+nseg2)*2                        ! dconc2
         nr_rar = nr_rar + 1
         if ( l_decl ) allocate ( dconc2  ( notot, noseg+nseg2 ), stat=ierr )
         if ( ierr .ne. 0 ) then ; write(lunrep,2010) "dconc2              " ; call srstop(1) ; endif
         if ( .not. l_decl ) write ( 328, 2040 ) nr_rar, "dconc2              ", notot*(noseg+nseg2)*2
      endif
      if ( .not. l_decl ) write ( 328, '(/5x,a20,i12)' ) "Total (4 byte words)",itota

      return

 2005 format ( ' ERROR  : real array is too big. Unable to create pointer. ' )
 2010 format ( ' ERROR  : allocating real array. Name   : ',A )
 2020 format (/' Parallel processing with ',i3,' processor(s)'/)
 2030 format ('  Parallel processing with ',i3,' processor(s)')
 2040 format (   i4,1x,a20,i12 )

      end subroutine
      end module

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      program type04

C     Goal: check the extension of the semantics analysis to logical
C     scalar variables

      logical l1, l2

      read *, l1

      l2 = l1 .AND. .NOT. l1

      read *, l1

      print *, l1, l2

      end

ext: f


text:
*--#[ log:
*	$Id: ffdcxs.f,v 1.7 1996/03/22 08:13:30 gj Exp $
*	$Log: ffdcxs.f,v $
c Revision 1.7  1996/03/22  08:13:30  gj
c Fixed bug in bugfix of ffdcxs.f
c
c Revision 1.6  1996/03/14  15:53:13  gj
c Fixed bug in ffcb0: cp in C, cma=cmb=0 was computed incorrectly.
c
c Revision 1.5  1996/03/13  15:43:36  gj
c Fixed bug, when ieps unknown already some things were computed and not zero'd.
c Now I first check ieps, and then compute.
c
c Revision 1.4  1995/12/08  10:38:16  gj
c Fixed too long line
c
*--#] log:
*###[ ffdcxs:
	subroutine ffdcxs(cs3,ipi12,y,z,dyz,d2yzz,dy2z,dyzzy,xpi,piDpj,
     +						ii,ns,isoort,ier)
***#[*comment:***********************************************************
*									*
*	calculates the the difference of two S's with y(3,4),z(3,4) and *
*	y(4)z(3)-y(3)z(4) given.  Note the difference with ffdcxs4, in	*
*	which the y's are the same and only the z's different.  Here	*
*	both can be different.	Also we skip an intermediate level.	*
*	Note also that this routine is much less conservative than	*
*	ffcxs3 in its expectations of the order of the roots: it knows	*
*	that it is (z-,z+,1-z-,1-z+)!					*
*									*
*	input:	y(4,3:4)	(real)	y,1-y in S with s3,s4		*
*		z(4,3:4)	(real)	z,1-z in S with s3,s4		*
*		dyz(2,2,3:4)	(real)	y - z				*
*		d2yzz(3:4)	(real)	2*y - z+ - z-			*
*		dy2z(4,3:4)	(real)	y - 2*z				*
*		dyzzy(4)	(real)	y(i,4)*z(i,4)-y(i,3)*z(i,4)	*
*		xpi(6,3:4)	(real)	usual				*
*		piDpj(6,3:4)	(real)	usual				*
*		cs3(40)	(complex)	assumed zero.			*
*									*
*	output: cs3(40)	(complex)	mod factors pi^2/12, in array	*
*		ipi12(6)(integer)	these factors			*
*		isoort(6)(integer)	returns kind of action taken	*
*		ier	(integer)	0=ok 1=inaccurate 2=error	*
*									*
*	calls:	ffcrr,ffcxr,real/dble,DCMPLX,log,ffadd1,ffadd2,ffadd3	*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	DOUBLE COMPLEX cs3(100)
	DOUBLE PRECISION y(4,3:4),z(4,3:4),dyz(2,2,3:4),d2yzz(3:4),
     +		dy2z(4,3:4),dyzzy(4),xpi(6,3:4),piDpj(6,6,3:4)
	integer ipi12(10),ii,ns,isoort(10),ier
*
*	local variables
*
	integer i,j,k,l,m,ier0,iepsi(4),iepsj(2,2),ipi12p(4),ipitot,
     +		ipitop
	logical normal
	DOUBLE COMPLEX cs1,cs2,cs1p,cs2p,cs3p(40),c
	DOUBLE PRECISION yy,zz,yy1,zz1,dyyzz,hulp3,hulp4,absc,xhck,xmax,
     +		rloss,xm1,xm2,xm1p,xm2p,x00(3)
	save iepsi
*
*	common blocks
*
	include 'ff.h'
*
*	data
*
	data iepsi /-2,+2,+2,-2/
*
*	statement function
*
	absc(c) = abs(DBLE(c)) +abs(DIMAG(c))
*
*	check constants
	if ( ltest .and. ns .ne. 6 ) print *,'ffdcxs: error: ns <> 6'
*  #] declarations:
*  #[ check input:
	if ( ltest ) then
	    rloss = xloss**2*DBLE(10)**(-mod(ier,50))
	    if ( lwrite ) print *,'rloss = ',rloss
	    do 2 k=3,4
		do 1 i=1,2
		    xhck = y(i,k) + y(i+2,k) - 1
		    xmax = max(abs(y(i,k)),x1)
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: 1 - y(',i,k,') <> 1-y(',i,k,'): ',
     +			y(i,k),y(i+2,k),xhck
		    xhck = z(i,k) + z(i+2,k) - 1
		    xmax = max(abs(z(i,k)),x1)
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: 1 - z(',i,k,') <> 1-z(',i,k,'): ',
     +			z(i,k),z(i+2,k),xhck
		    xhck = dyz(2,i,k) - y(2,k) + z(i,k)
		    xmax = max(abs(y(2,k)),abs(z(i,k)))
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: dyz(2',i,k,')<>y(2',k,')-z(',i,k,'): ',
     +			dyz(2,i,k),y(2,k),z(i,k),xhck
		    xhck = dy2z(i,k) - y(2,k) + 2*z(i,k)
		    xmax = max(abs(y(2,k)),2*abs(z(i,k)))
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: dy2z(',i,k,')<>y(2',k,')-2*z(',i,k,
     +			'): ',dy2z(i,k),y(2,k),2*z(i,k),xhck
		    xhck = dy2z(i+2,k) - y(4,k) + 2*z(i+2,k)
		    xmax = max(abs(y(4,k)),2*abs(z(i+2,k)))
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: dy2z(',i+2,k,')<>y(4',k,')-2z(',i+2,k,
     +			'): ',dy2z(i+2,k),y(4,k),2*z(i+2,k),xhck
		    l = 2*k+i - 6
		    m = 2*(k/2)
		    xhck = dyzzy(l) - y(m,4)*z(m+i-2,3) +
     +			y(m,3)*z(m+i-2,4)
		    xmax = max(abs(dyzzy(l)),abs(y(m,4)*z(m+i-2,3)))
		    if ( rloss*xhck .gt. precx*xmax ) print *,'ffdcxs:',
     +			' error: dyzzy(',l,') <> ...',dyzzy(l),
     +			y(m,4)*z(m+i-2,3),y(m,3)*z(m+i-2,4),xhck
    1		continue
    2	    continue
	endif
*  #] check input:
*  #[ normal case:
	normal = .FALSE.
   10	continue
	if ( normal .or. isoort(1) .ne. isoort(9) .or. isoort(1) .lt.
     +		10 ) then
	    if ( lwrite ) print *,'ffdcxs: normal case'
	    call ffcxs3(cs3( 1),ipi12(1),y(1,3),z(1,3),dyz(1,1,3),
     +		d2yzz(3),dy2z(1,3),xpi(1,3),piDpj(1,1,3),ii,6,
     +		isoort(1),ier)
	    call ffcxs3(cs3(81),ipi12(9),y(1,4),z(1,4),dyz(1,1,4),
     +		d2yzz(4),dy2z(1,4),xpi(1,4),piDpj(1,1,4),ii,6,
     +		isoort(9),ier)
	    return
	endif
*  #] normal case:
*  #[ rotate R's:
	if ( abs(y(2,3)) .lt. 1/xloss ) then
	    if ( lwrite ) print *,'ffdcxs: rotating R''s'
	    do 102 i=1,2
	    do 101 j=1,2
*		iepsi() = /-2,+2,+2,-2/
*		BUT I AM NOT YET SURE OF THE SIGNS (29/6/89)
		k = 2*(i-1)+j
		if ( y(2*i,3) .gt. 0 ) then
		    iepsj(j,i) = iepsi(k)
		else
		    iepsj(j,i) = -iepsi(k)
		endif
		if ( y(2*i,3) .gt. 0 .neqv. y(2*i,4) .gt. 0 ) then
*		    I have no clue to the ieps, take normal route
*		    iepsj(j,i) = 0
		    if ( lwrite ) print *,'ffdcxs: don''t know ieps ',i
		    normal = .TRUE.
		    goto 10
		endif
  101	    continue
  102	    continue
*	    loop over y,z , 1-y,1-z
	    do 120 i=1,2
*	    loop over z+ , z-
	    do 110 j=1,2
		if ( j .eq. 2 ) then
*		    do not calculate if not there (isoort=0, one root)
*			(this is probably not needed as this case should
*			have been dealt with in ffdxc0)
		    if ( isoort(9) .eq. 0 ) then
			if ( lwrite ) then
			    print *,'ffdcxs: isoort(9).eq.0, not needed'
			endif
			goto 110
		    endif
*		    or if not needed (isoort=2, two equal roots)
		    if ( mod(isoort(9),10) .eq. 2 ) then
			if ( lwrite ) print *,'ffdcxs: skipped next ',
     +				'R as it is equal'
*			we use that l still contains the correct value
			do 105 m=1,7
			    cs3(10*(l-1)+m) = 2*DBLE(cs3(10*(l-1)+m))
  105			continue
			ipi12(l) = 2*ipi12(l)
			goto 110
		    endif
		endif
		k = 2*(i-1)+j
		l = 8*(i-1)+j
		if ( dyzzy(k) .ne. 0 ) then
*		    minus sign wrong in thesis (2.78)
		    hulp3 = -dyz(2,j,3)/dyzzy(k)
		    hulp4 = +dyz(2,j,4)/dyzzy(k)
		    yy = y(2*i,3)*hulp4
		    yy1 = y(2*i,4)*hulp3
		    zz = z(k,3)*hulp4
		    zz1 = z(k,4)*hulp3
		    dyyzz = dyz(2,j,3)*hulp4
		    if ( i .eq. 2 ) then
			yy = -yy
			yy1 = -yy1
			zz = -zz
			zz1 = -zz1
		    endif
*		    if ( ltest ) then
*			if ( rloss*abs(yy+yy1-1) .gt. precx*max(abs(yy),
*     +				x1) ) print *,'ffdcxs: error: 1 - yy ',
*     +				'<> yy1',yy,yy1,yy+yy1-1
*			if ( rloss*abs(zz+zz1-1) .gt. precx*max(abs(zz),
*     +				x1) ) print *,'ffdcxs: error: 1 - zz ',
*     +				'<>zz1',zz,zz1,zz+zz1-1
*			if ( rloss*abs(dyyzz-yy+zz) .gt. precx*max(abs(
*     +				yy),abs(zz)) ) print *,'ffdcxs: error:',
*				' dyyzz<>yy-zz',dyyzz,yy,zz,dyyzz-yy+zz
*		    endif
		    if ( lwrite ) then
			do 109 m=3,4
			    print *,'arg1',m,' was ',+y(2,m)/dyz(2,k,m)
			    print *,'arg2',m,' was ',-y(4,m)/dyz(2,k,m)
  109			continue
			print *,'arg1',m,'  is ',+yy/dyyzz
			print *,'arg2',m,'  is ',-yy1/dyyzz
		    endif
		    call ffcxr(cs3(10*l-9),ipi12(l),yy,yy1,zz,zz1,dyyzz,
     +			.FALSE.,x0,x0,x0,.FALSE.,x00,iepsj(j,i),ier)
		else
		    if ( lwrite ) print *,'  y(4)z(3)-y(3)z(4)=0 -> S=0'
		endif
  110	    continue
  120	    continue
	    goto 800
	endif
*  #] rotate R's:
*  #[ other cases (not ready):
	if ( lwrite ) print *,'ffdcxs: warning: special case not',
     +		' yet implemented, trying normal route'
	call ffcxs3(cs3( 1),ipi12(1),y(1,3),z(1,3),dyz(1,1,3),
     +		d2yzz(3),dy2z(1,3),xpi(1,3),piDpj(1,1,3),ii,ns,
     +		isoort(1),ier)
	call ffcxs3(cs3(81),ipi12(9),y(1,4),z(1,4),dyz(1,1,4),
     +		d2yzz(4),dy2z(1,4),xpi(1,4),piDpj(1,1,4),ii,ns,
     +		isoort(9),ier)
	return
*  #] other cases (not ready):
*  #[ debug:
  800	if ( lwrite ) then
	    ier0 = 0
	    do 805 i=1,40
		cs3p(i) = 0
  805	    continue
	    print '(a)','   #[ compare: '
	    call ffcxs3(cs3p( 1),ipi12p(1),y(1,3),z(1,3),dyz(1,1,3),
     +		d2yzz(3),dy2z(1,3),xpi(1,3),piDpj(1,1,3),ii,ns,
     +		isoort(1),ier0)
	    call ffcxs3(cs3p(21),ipi12p(3),y(1,4),z(1,4),dyz(1,1,4),
     +		d2yzz(4),dy2z(1,4),xpi(1,4),piDpj(1,1,4),ii,ns,
     +		isoort(9),ier0)
	    print '(a)','   #] compare: '
	    cs1 = 0
	    cs2 = 0
	    cs1p = 0
	    cs2p = 0
	    xm1 = 0
	    xm2 = 0
	    xm1p = 0
	    xm2p = 0
	    do 810 i=1,20
		cs1 = cs1 + cs3(i)
		xm1 = max(xm1,absc(cs1))
		cs2 = cs2 + cs3(i+80)
		xm2 = max(xm2,absc(cs2))
		cs1p = cs1p + cs3p(i)
		xm1p = max(xm1p,absc(cs1p))
		cs2p = cs2p + cs3p(i+20)
		xm2p = max(xm2p,absc(cs2p))
  810	    continue
	    ipitot = ipi12(1) + ipi12(2) - ipi12(9) - ipi12(10)
	    ipitop = ipi12p(1) + ipi12p(2) - ipi12p(3) - ipi12p(4)
 1000	    format(2g24.16,g12.4)
	    print *,'ffdcxs: compare:'
	    print *,'    Originally:'
	    print 1000,cs1p,xm1p
	    print 1000,-cs2p,xm2p
	    if ( ipitot .ne. 0 ) print 1000,ipitot*pi12,0.
	    print *,'+ ------------'
	    print 1000,cs1p-cs2p+ipitot*DBLE(pi12),max(xm1p,xm2p)
	    print *,'    Now:'
	    print 1000,cs1,xm1
	    print 1000,-cs2,xm2
	    if ( ipitop .ne. 0 ) print 1000,ipitop*pi12,0.
	    print *,'+ ------------'
	    print 1000,cs1-cs2+ipitop*DBLE(pi12),max(xm1,xm2)
	endif
*  #] debug:
*###] ffdcxs:
	end
*###[ ffdcs:
	subroutine ffdcs(cs3,ipi12,cy,cz,cdyz,cd2yzz,cdyzzy,cdyyzz,
     +		cpi,cpiDpj,ii,ns,isoort,ier)
***#[*comment:***********************************************************
*									*
*	calculates the the difference of two S's with cy(3,4),cz(3,4),	*
*	cy(4)cz(3)-cy(3)cz(4) given.  Note the difference with ffdcs4,	*
*	in which the cy's are the same and only the cz's different.	*
*	Here both can be different.	Also we skip an intermediat	*
*	level.								*
*									*
*	input:	cy(4,3:4)    (complex)	cy,1-cy in S with s3,s4		*
*		cz(4,3:4)    (complex)	cz,1-cz in S with s3,s4		*
*		cdyz(2,2,3:4)(complex)	cy - cz				*
*		cd2yzz(3:4)  (complex)	2*cy - cz+ - cz-		*
*		cdyzzy(4)    (complex)	cy(i,4)*cz(i,4)-cy(i,3)*cz(i,4)	*
*		cdyyzz(2)    (complex)	cy(i,4)-cz(i,4)-cy(i,3)+cz(i,4)	*
*		cpi(6,3:4)   (complex)	usual				*
*		cpiDpj(6,3:4)(complex)	usual				*
*		cs3(40)	     (complex)	assumed zero.			*
*									*
*	output: cs3(40)	     (complex)	mod factors pi^2/12, in array	*
*		ipi12(6)     (integer)	these factors			*
*		isoort(6)    (integer)	returns kind of action taken	*
*		ier	     (integer)	number of digits lost		*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	DOUBLE COMPLEX cs3(100)
	DOUBLE COMPLEX cy(4,3:4),cz(4,3:4),cdyz(2,2,3:4),cd2yzz(3:4),
     +		cdyzzy(4),cdyyzz(2),cpi(6,3:4),cpiDpj(6,6,3:4)
	integer ipi12(10),ii,ns,isoort(10),ier
*
*	local variables
*
	integer i,j,k,l,m,n,ier0,ieps,ni(4,3:4),ntot(3:4),
     +		n1a,n1b,ii1,nffeta,nffet1,i2pi,n2a,ip,ipi12p(4),ipitot,
     +		ipitop
	DOUBLE COMPLEX cs1,cs2,cs1p,cs2p,cs3p(40),c,cc,clogy,zfflog,
     +		zfflo1,cmip,yy,zz,yy1,zz1,dyyzz,hulp3,hulp4,xhck
	DOUBLE PRECISION rloss,xm1,xm2,xm1p,xm2p,absc,xmax,s1,s2,s3,s4,
     +		y1m,y1m1,y1p,y1p1
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) +abs(DIMAG(c))
*
*	check constants
	if ( ltest .and. ns .ne. 6 ) print *,'ffdcs: error: ns <> 6'
*  #] declarations:
*  #[ check input:
	if ( ltest ) then
	    rloss = xloss*DBLE(10)**(-mod(ier,50)-2)
	    if ( lwrite ) print *,'rloss = ',rloss
	    do 20 k=3,4
		do 10 i=1,2
		    xhck = cy(i,k) + cy(i+2,k) - 1
		    xmax = max(absc(cy(i,k)),x1)
		    if ( rloss*absc(xhck) .gt. precc*xmax )
     +			print *,'ffdcs: error: 1 - cy(',i,k,') <> 1-cy('
     +			,i,k,'): ',cy(i,k),cy(i+2,k),xhck
		    xhck = cz(i,k) + cz(i+2,k) - 1
		    xmax = max(absc(cz(i,k)),x1)
		    if ( rloss*absc(xhck) .gt. precc*xmax )
     +			print *,'ffdcs: error: 1 - cz(',i,k,') <> 1-cz('
     +			,i,k,'): ',cz(i,k),cz(i+2,k),xhck
		    xhck = cdyz(2,i,k) - cy(2,k) + cz(i,k)
		    xmax = max(absc(cy(2,k)),absc(cz(i,k)))
		    if ( rloss*absc(xhck) .gt. precc*xmax )
     +			print *,'ffdcs: error: cdyz(2',i,k,')<>cy(2',k,
     +			')-cz(',i,k,'): ',cdyz(2,i,k),cy(2,k),cz(i,k),
     +			xhck
		    l = 2*k+i - 6
		    m = 2*(k/2)
		    xhck = cdyzzy(l) - cy(m,4)*cz(m+i-2,3) +
     +			cy(m,3)*cz(m+i-2,4)
		    xmax =max(absc(cdyzzy(l)),absc(cy(m,4)*cz(m+i-2,3)))
		    if ( rloss*absc(xhck) .gt. precc*xmax )
     +			print *,'ffdcs: error: cdyzzy(',l,') <> ...',
     +			cdyzzy(l),cy(m,4)*cz(m+i-2,3),cy(m,3)*cz(m+i-2,4
     +			),xhck,ier
   10		continue
   20	    continue
	endif
*  #] check input:
*  #[ normal case:
	if ( mod(isoort(1),5).ne.mod(isoort(9),5) .or. isoort(1).gt.-5
     +								) then
	    if ( lwrite ) print *,'ffdcs: normal case'
	    if ( ltest .and. isoort(1) .le. -100 ) then
		print *,'ffdcs: error: wrong value for isoort'
	    endif
	    call ffcs3(cs3( 1),ipi12(1),cy(1,3),cz(1,3),cdyz(1,1,3),
     +		cd2yzz(3),cpi(1,3),cpiDpj(1,1,3),ii,6,isoort(1),ier)
	    call ffcs3(cs3(81),ipi12(9),cy(1,4),cz(1,4),cdyz(1,1,4),
     +		cd2yzz(4),cpi(1,4),cpiDpj(1,1,4),ii,6,isoort(9),ier)
	    return
	endif
*  #] normal case:
*  #[ rotate R's:
	if ( absc(cy(2,3)) .lt. 1/xloss .or. isoort(1) .le. -100 ) then
	    if ( lwrite ) print *,'ffdcs: rotated R''s'
*
*	    loop over cy,cz , 1-cy,1-cz
	    do 190 i=1,2

	    if ( isoort(1).le.-100 .and. i.eq.2 ) then
*
*		special case del2s=0, a limit has been taken
*
		if ( ii .eq. 2  ) then
*
*		    we took the wrong sign for the dilogs...
*
		    do 110 j=1,20
			cs3(j) = -cs3(j)
  110		    continue
		    ipi12(1) = -ipi12(1)
		    ipi12(2) = -ipi12(2)
		endif
*
*		now the remaining logs.  take care to get the ieps
*		correct!
*
		if ( i.eq.1 .eqv. DBLE(cy(2*i,3)).gt.0 ) then
		    ieps = -3
		else
		    ieps = +3
		endif
		call ffclg2(cs3(81),ipi12(9),cy(2,3),cz(1,3),
     +			cdyz(2,1,3),cy(2,4),cz(1,4),cdyz(2,1,4),
     +			cdyyzz(1),isoort(1),isoort(9),ii,ieps,ier)
		if ( ii .eq. 2 ) then
*		    we have the wrong sign
		    do 120 j=81,83
			cs3(j) = -cs3(j)
  120		    continue
		    ipi12(9) = -ipi12(9)
		endif
		if ( mod(isoort(1),5).eq.0 .and. mod(isoort(9),5).eq.0
     +								) then
		    if ( lwrite ) print *,'ffdcs: skipped other logs ',
     +			'as they are the complex conjugate'
		    do 130 j=81,83
			cs3(j) = 2*DBLE(cs3(j))
  130		    continue
		    ipi12(9) = 2*ipi12(9)
		else
		    print *,'ffdcs: error: not yet tested'
		    call ffclg2(cs3(91),ipi12(10),cy(2,3),cz(2,3),
     +			cdyz(2,2,3),cy(2,4),cz(2,4),cdyz(2,2,4),
     +			cdyyzz(2),isoort(1),isoort(9),ii,-ieps,ier)
		    if ( ii .eq. 2 ) then
*			we have the wrong sign
			do 140 j=91,93
			    cs3(j) = -cs3(j)
  140			continue
			ipi12(10) = -ipi12(10)
		    endif
		endif
		goto 190
	    endif
*
*	    loop over cz- , cz+
	    do 180 j=1,2
		if ( j .eq. 2 ) then
		    if ( isoort(9) .eq. 0 .or. isoort(1) .eq. 0 ) then
*
*			(this is not correct as this case should
*			have been dealt with in ffdxc0,ffdcc0)
*
			call fferr(79,ier)
			goto 180
		    elseif ( mod(isoort(9),5) .eq. 0 .and.
     +			     mod(isoort(1),5) .eq. 0 ) then
*
*			or if not needed (isoort=-10, two conjugate roots)
*
			if ( lwrite ) print *,'ffdcs: skipped next ',
     +				'R as it is the conjugate'
*			we use that l still contains the correct value
			do 150 m=1,9
			    cs3(10*(l-1)+m) = 2*DBLE(cs3(10*(l-1)+m))
  150			continue
			ipi12(l) = 2*ipi12(l)
			goto 180
		    elseif ( mod(isoort(9),10) .eq. 2 ) then
			if ( lwrite ) print *,'ffdcs: skipped next ',
     +				'R as it is equal'
*			we use that l still contains the correct value
			do 160 m=1,9
			    cs3(10*(l-1)+m) = 2*cs3(10*(l-1)+m)
  160			continue
			ipi12(l) = 2*ipi12(l)
			goto 180
		    endif
		endif
		k = 2*(i-1)+j
		l = 8*(i-1)+j
		if ( cdyzzy(k) .ne. 0 ) then
		    hulp3 = -cdyz(2,j,3)/cdyzzy(k)
		    hulp4 = cdyz(2,j,4)/cdyzzy(k)
		    yy = cy(2*i,3)*hulp4
		    yy1 = cy(2*i,4)*hulp3
		    zz = cz(k,3)*hulp4
		    zz1 = cz(k,4)*hulp3
		    dyyzz = cdyz(2,j,3)*hulp4
		    if ( i .eq. 2 ) then
			yy = -yy
			yy1 = -yy1
			zz = -zz
			zz1 = -zz1
		    endif
*
*		    ieps = 3 means: dear ffcrr, do not use eta terms,
*		    they are calculated here.  The sign gives the sign
*		    of the imag. part of the argument of the dilog, not
*		    y-z.
*
		    if ( i.eq.1 .eqv. j.eq.1 .eqv. DBLE(cy(2*i,3)).gt.0
     +								) then
			ieps = -3
		    else
			ieps = +3
		    endif
		    call ffcrr(cs3(10*l-9),ipi12(l),yy,yy1,zz,zz1,dyyzz,
     +			.FALSE.,c0,c0,c0,isoort(j),ieps,ier)
*
*		    eta terms of the R's (eta(.)*log(c1)-eta(.)*log(c2))
*
		    do 170 m=3,4
*			no eta terms in the real case
			if ( DIMAG(cz(k,m)) .eq. 0 .and.
     +			     DIMAG(cdyz(2,j,m)) .eq. 0 ) then
			    ni(k,m) = 0
			elseif ( i .eq. 1 ) then
			    ni(k,m) = nffeta(-cz(k,m),1/cdyz(2,j,m),ier)
			else
			    ni(k,m) = nffeta(cz(k,m),1/cdyz(2,j,m),ier)
			endif
  170		    continue
		    if ( ni(k,3) .ne. 0 .or. ni(k,4) .ne. 0 ) then
			if ( lwrite ) print *,'n3,n4: ',ni(k,3),ni(k,4)
			if ( ni(k,3) .ne. ni(k,4) ) then
			    do 175 m=3,4
				c = cy(2*i,m)/cdyz(2,j,m)
				if ( i .eq. 2 ) c = -c
				cc = c-1
				if ( absc(cc) .lt. xloss ) then
				    if ( lwrite ) print *,'c = ',c
				    c = cz(k,m)/cdyz(2,j,m)
				    if ( lwrite ) print *,'c+= ',1-c
				    clogy = zfflo1(c,ier)
				else
				    clogy = zfflog(c,0,c0,ier)
				endif
				n = 10*l + (m-3) - 2
				if ( ltest .and. cs3(n) .ne. 0 ) then
				    print *,'ffdcs: error: cs3(',n,
     +					') != 0'
				endif
				if ( m .eq. 3 ) then
				   cs3(n) = + ni(k,m)*c2ipi*clogy
				else
				   cs3(n) = - ni(k,m)*c2ipi*clogy
				endif
				if ( lwrite ) then
				print *,'eta',n,'= ',ni(k,m)*c2ipi*clogy
				if ( m .eq. 4 ) print *,'som = ',cs3(n)
     +							+ cs3(n-1)
				endif
 175			    continue
			else
			    if ( i .eq. 1 ) then
				n1a = nffeta(cy(k,3)/cdyz(2,j,3),
     +					cdyz(2,j,4)/cy(k,4),ier)
			    else
				n1a = nffeta(-cy(k,3)/cdyz(2,j,3),
     +					-cdyz(2,j,4)/cy(k,4),ier)
			    endif
			    if ( n1a .ne. 0 ) then
				call fferr(80,ier)
			    endif
			    c =cy(k,3)*cdyz(2,j,4)/(cdyz(2,j,3)*cy(k,4))
			    cc = c-1
			    if ( absc(cc) .lt. xloss ) then
				if ( lwrite ) print *,'1-c = ',1-c
				c = -cdyzzy(k)/(cdyz(2,j,3)*cy(k,4))
				if ( lwrite ) print *,'1-c+= ',c
				clogy = zfflo1(c,ier)
			    else
				clogy = zfflog(c,0,c0,ier)
			    endif
			    n = 10*l - 2
			    if ( ltest .and. cs3(n) .ne. 0 ) then
				print *,'ffdcs: error: cs3(',n,') not 0'
			    endif
			    if ( i .eq. 1 ) then
				cs3(n) = +ni(k,3)*c2ipi*clogy
			    else
				cs3(n) = -ni(k,3)*c2ipi*clogy
			    endif
			    if ( lwrite ) print *,'both etas ',cs3(n)
			endif
		    endif
		else
		    if ( lwrite ) print *,'  cy(4)cz(3)-cy(3)cz(4)=0',
     +			' -> S=0'
		endif
  180	    continue
  190	    continue
	    goto 700
	endif
*  #] rotate R's:
*  #[ other cases (not ready):
	if ( lwrite ) print *,'ffdcs: warning: special case not',
     +		' yet implemented, trying normal route'
	call ffcs3(cs3( 1),ipi12(1),cy(1,3),cz(1,3),cdyz(1,1,3),
     +		cd2yzz(3),cpi(1,3),cpiDpj(1,1,3),ii,ns,isoort(1),ier)
	call ffcs3(cs3(81),ipi12(9),cy(1,4),cz(1,4),cdyz(1,1,4),
     +		cd2yzz(4),cpi(1,4),cpiDpj(1,1,4),ii,ns,isoort(9),ier)
	return
*  #] other cases (not ready):
*  #[ get eta's:
  700	continue
	ip = ii+3
	do 740 k=3,4
	l = 8*(k-3) + 1
	if ( DIMAG(cpi(ip,k)) .eq. 0 ) then
*
*	complex because of a complex root in y or z
*
	if ( (mod(isoort(l),10).eq.-1 .or. mod(isoort(l),10).eq.-3)
     +		.and. isoort(l+1) .ne. 0 ) then
*
*	    isoort = -1:    y is complex, possibly z as well
*	    isoort = -3:    y,z complex, but (y-z-)(y-z+) real
*	    isoort = 0:     y is complex, one z root only
*	    isoort = -10:   y is real, z is complex
*	    isoort = -5,-6: y,z both real
*
	    cmip = DCMPLX(DBLE(x0),-DBLE(cpi(ip,k)))
	    if ( DIMAG(cz(1,k)) .eq. 0 ) then
		ni(1,k) = 0
	    else
		ni(1,k) = nffet1(-cz(1,k),-cz(2,k),cmip,ier)
		i = nffet1(cz(3,k),cz(4,k),cmip,ier)
		if ( i .ne. ni(1,k) ) call fferr(53,ier)
	    endif
	    ni(2,k) = 0
	    if ( DBLE(cd2yzz(k)).eq.0 .and. ( DIMAG(cz(1,k)).eq.0 .and.
     +		 DIMAG(cz(2,k)).eq.0 .or. DBLE(cdyz(2,1,k)).eq.0 .and.
     +		 DBLE(cdyz(2,2,k)) .eq. 0 ) ) then
*		follow the i*epsilon prescription as (y-z-)(y-z+) real
		if ( DBLE(cpi(ip,k)) .lt. 0 ) then
		    ni(3,k) = -1
		else
		    ni(3,k) = 0
		endif
		ni(4,k) = -nffet1(cdyz(2,1,k),cdyz(2,2,k),cmip,ier)
	    else
		if ( DBLE(cpi(ip,k)) .lt. 0 .and. DIMAG(cdyz(2,1,k)*
     +			cdyz(2,2,k)) .lt. 0 ) then
		    ni(3,k) = -1
		else
		    ni(3,k) = 0
		endif
		ni(4,k) = -nffeta(cdyz(2,1,k),cdyz(2,2,k),ier)
	    endif
	elseif ( (mod(isoort(l),10).eq.-1 .or. mod(isoort(l),10).eq.-3)
     +		.and. isoort(l+1).eq.0 ) then
	    ni(1,k) = 0
	    if ( DIMAG(cz(1,k)) .ne. 0 ) then
		ni(2,k) = nffet1(-cpiDpj(ii,ip,k),-cz(1,k),DCMPLX(DBLE(0
     +			),DBLE(-1)),ier)
	    else
		ni(2,k) = nffet1(-cpiDpj(ii,ip,k),DCMPLX(DBLE(0),
     +			DBLE(1)),DCMPLX(DBLE(0),DBLE(-1)),ier)
	    endif
	    ni(3,k) = 0
	    ni(4,k) = -nffeta(-cpiDpj(ii,ip,k),cdyz(2,1,k),ier)
	else
	    if ( mod(isoort(l),5).ne.0 .and. mod(isoort(l),5).ne.-1
     +			.and. mod(isoort(l),5).ne.-3 ) then
		call fferr(81,ier)
		print *,'isoort(',l,') = ',isoort(l)
	    endif
	    ni(1,k) = 0
	    ni(2,k) = 0
	    ni(3,k) = 0
	    ni(4,k) = 0
	endif
	else
	    print *,'ffdcs: error: cpi complex should not occur'
	    stop
	endif
  740	continue
	if ( lwrite ) then
	    print *,'ffdcs: eta''s are: '
	    print *,'s3: ',(ni(i,3),i=1,4)
	    print *,'s4: ',(ni(i,4),i=1,4)
	endif
*  #] get eta's:
*  #[ add eta's:
	do 750 k=3,4
	    ntot(k) = ni(1,k)+ni(2,k)+ni(3,k)+ni(4,k)
  750	continue
	if ( ntot(3) .ne. 0 .and. ntot(3) .eq. ntot(4) ) then
	    if ( lwrite ) print *,'ffdcs: warning: could be smarter...'
	endif
	do 760 k=3,4
	    if ( ntot(k) .ne. 0 ) call ffclgy(cs3(20+80*(k-3)),
     +		ipi12(2+8*(k-3)),ni(1,k),cy(1,k),cz(1,k),cd2yzz(k),ier)
  760	continue
*  #] add eta's:
*  #[ debug:
  800	if ( lwrite ) then
	    ier0 = 0
	    do 805 i=1,40
		cs3p(i) = 0
  805	    continue
	    do 806 i=1,4
		ipi12p(i) = 0
  806		continue
	    if ( isoort(1) .gt. -100 ) then
		print '(a)','   #[ compare: '
		call ffcs3(cs3p( 1),ipi12p(1),cy(1,3),cz(1,3),cdyz(1,1,3
     +			),cd2yzz(3),cpi(1,3),cpiDpj(1,1,3),ii,ns,
     +			isoort(1),ier0)
		call ffcs3(cs3p(21),ipi12p(3),cy(1,4),cz(1,4),cdyz(1,1,4
     +			),cd2yzz(4),cpi(1,4),cpiDpj(1,1,4),ii,ns,
     +			isoort(9),ier0)
		print '(a)','   #] compare: '
	    endif
	    cs1 = 0
	    cs2 = 0
	    cs1p = 0
	    cs2p = 0
	    xm1 = 0
	    xm2 = 0
	    xm1p = 0
	    xm2p = 0
	    do 810 i=1,20
		cs1 = cs1 + cs3(i)
		xm1 = max(xm1,absc(cs1))
		cs2 = cs2 + cs3(i+80)
		xm2 = max(xm2,absc(cs2))
		cs1p = cs1p + cs3p(i)
		xm1p = max(xm1p,absc(cs1p))
		cs2p = cs2p + cs3p(i+20)
		xm2p = max(xm2p,absc(cs2p))
  810	    continue
	    ipitot = ipi12(1) + ipi12(2) - ipi12(9) - ipi12(10)
	    ipitop = ipi12p(1) + ipi12p(2) - ipi12p(3) - ipi12p(4)
 1000	    format(2g24.16,g12.4)
	    print *,'ffdcs: compare:'
	    print *,'    Originally:'
	    print 1000,cs1p,xm1p
	    print 1000,-cs2p,xm2p
	    if ( ipitop .ne. 0 ) print 1000,ipitop*DBLE(pi12),0.
	    print *,'+ ------------'
	    print 1000,cs1p-cs2p+ipitop*DBLE(pi12),max(xm1p,xm2p)
	    print *,'    Now:'
	    print 1000,cs1,xm1
	    print 1000,-cs2,xm2
	    if ( ipitot .ne. 0 ) print 1000,ipitot*pi12,0.
	    print *,'+ ------------'
	    print 1000,cs1-cs2+ipitot*DBLE(pi12),max(xm1,xm2)
	endif
*  #] debug:
*###] ffdcs:
	end
*###[ ffclg2:
	subroutine ffclg2(cs3,ipi12,cy3,cz3,cdyz3,cy4,cz4,cdyz4,cdyyzz,
     +		isort3,isort4,ii,ieps,ier)
***#[*comment:***********************************************************
*									*
*	Calculate the finite part of the divergent dilogs in case	*
*	del2s=0.  These are given by					*
*									*
*		log^2(-cdyz3)/2 - log^2(-cdyz4)/2			*
*									*
*	Note that often we only need the imaginary part, which may be	*
*	very unstable even if the total is not.				*
*									*
*									*
*	Input:	cy3,cz3,cdyz3	(complex)	y,z,diff in C with s3	*
*		cy4,cz4,cdyz4	(complex)	y,z,diff in C with s4	*
*		cdyyzz		(complex)	y4 - z4 - y3 + z3	*
*		isort3,4	(integer)				*
*									*
*	Output	cs3(4)		(complex)	output			*
*		ipi12		(integer)	terms pi^2/12		*
*		ier		(integer)	error flag		*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	DOUBLE COMPLEX cs3(3),cy3(3),cz3(3),cdyz3,cy4(3),cz4(3),cdyz4,
     +		cdyyzz
	integer ipi12,ieps,ier,isort3,isort4,ii
*
*	local variables
*
	integer n1,nffeta,nffet1,ipi3,ipi4
	DOUBLE COMPLEX c,cc,chck,clog3,clog4,clog1,zfflo1,cipi
	DOUBLE PRECISION absc,rloss
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ check input:
	if ( ltest ) then
	    if ( abs(ieps) .ne. 3 ) print *,'ffclg2: error: |ieps| <> 3'
	    rloss = xloss**2*DBLE(10)**(-mod(ier,50))
	    if ( cs3(1) .ne. c0 .or. cs3(2) .ne. c0 .or. cs3(3) .ne. c0)
     +		print *,'ffclg2: error: cs3 <> 0',cs3
	    chck = cz3(1) + cz3(3) - 1
	    if ( rloss*absc(chck) .gt. precc*max(x1,absc(cz3(1))) )
     +		print *,'ffclg2: error: 1 - cz3 <> (1-cz3)',cz3(1),cz3(3
     +		),chck
	    chck = cz4(1) + cz4(3) - 1
	    if ( rloss*absc(chck) .gt. precc*max(x1,absc(cz4(1))) )
     +		print *,'ffclg2: error: 1 - cz4 <> (1-cz4)',cz4(1),cz4(3
     +		),chck
	    chck = cdyz3 - cy3(1) + cz3(1)
	    if ( rloss*absc(chck) .gt. precc*max(absc(cdyz3),absc(cz3(1)
     +		)) ) print *,'ffclg2: error: cdyz3 <> cy3-cz3',
     +		cdyz3,cy3(1),cz3(1),chck
	    chck = cdyz4 - cy4(1) + cz4(1)
	    if ( rloss*absc(chck) .gt. precc*max(absc(cdyz4),absc(cz4(1)
     +		)) ) print *,'ffclg2: error: cdyz4 <> cy4-cz4',
     +		cdyz4,cy4(1),cz4(1),chck
	    chck = cdyyzz - cy4(1) + cz4(1) + cy3(1) - cz3(1)
	    if ( rloss*absc(chck) .gt. precc*max(absc(cy4(1)),absc(cz4(1
     +		)),absc(cy3(1))) ) print *,'ffclg2: error: cdyyzz <> ',
     +		'terms',cdyyzz,cy4(1),cz4(1),cy3(1),cz3(1),chck
	endif
*  #] check input:
*  #[ calculations:
	cipi = DCMPLX(DBLE(x0),DBLE(pi))
	if ( DBLE(cdyz3) .lt. 0 ) then
	    clog3 = log(-cdyz3)
	    ipi3 = 0
	else
	    clog3 = log(cdyz3)
	    if ( DIMAG(cdyz3) .gt. 0 ) then
		ipi3 = -1
	    elseif ( DIMAG(cdyz3) .lt. 0 ) then
		ipi3 = +1
	    else
		ipi3 = sign(1,-ieps)
	    endif
	endif
	if ( DBLE(cdyz4) .lt. 0 ) then
	    clog4 = log(-cdyz4)
	    ipi4 = 0
	else
	    clog4 = log(cdyz4)
	    if ( DIMAG(cdyz4) .gt. 0 ) then
		ipi4 = -1
	    elseif ( DIMAG(cdyz4) .lt. 0 ) then
		ipi4 = +1
	    else
		ipi4 = sign(1,-ieps)
	    endif
	endif
	cc = clog3-clog4
	if ( absc(cc) .ge. xloss*absc(clog3) ) then
	    cs3(1) = -(clog3+ipi3*cipi)**2/2
	    cs3(2) = +(clog4+ipi4*cipi)**2/2
	    if ( lwrite ) clog1 = -123
	else
	    c = cdyyzz/cdyz4
	    clog1 = zfflo1(c,ier)
*
*	    notice that zfflog return log(a-ieps) (for compatibility
*	    with the dilog)		   ^
*
	    if ( DIMAG(cdyz3) .eq. 0 ) then
		n1 = nffet1(DCMPLX(DBLE(0),DBLE(-ieps)),-1/cdyz4,-c,
     +			ier)
	    elseif ( DIMAG(cdyz3) .eq. 0 ) then
		n1 = nffet1(-cdyz3,DCMPLX(DBLE(0),DBLE(ieps)),-c,ier)
	    else
		n1 = nffeta(-cdyz3,-1/cdyz4,ier)
	    endif
	    if ( n1 .ne. 0 ) then
		clog1 = clog1 - n1*c2ipi
	    endif
	    cs3(1) = -clog3*clog1/2
	    cs3(2) = -clog4*clog1/2
	    cs3(3) = -(ipi3+ipi4)*cipi*clog1/2
*	    we could split off a factor 2*pi^2 if needed
	endif
*	ATTENTION: now (23-jul-1989) ffdcs assumes that only *3* cs are
*	set. Change ffdcs as well if this is no longer true!
*  #] calculations:
*  #[ debug:
	if ( lwrite ) then
	    if ( clog1 .ne. -123 ) then
		print *,'ffclg2: originally:'
		print '(a,2g24.15)','S3: ',+(clog3+ipi3*cipi)**2/2
		print '(a,2g24.15)','S4: ',-(clog4+ipi4*cipi)**2/2
		print '(a,2g24.15,2i6)','sum:',+(clog3+ipi3*cipi)**2/2
     +			-(clog4+ipi4*cipi)**2/2
	    endif
	    print *,'ffclg2: now:'
	    print '(a,2g24.15)','S3: ',-cs3(1)
	    print '(a,2g24.15)','S4: ',-cs3(2)
	    print '(a,2g24.15)','Spi:',-cs3(3)
	    print '(a,2g24.15,2i6)','sum:',-cs3(1)-cs3(2)-cs3(3),-ipi12,
     +		ier
	endif
*  #] debug:
*###] ffclg2:
	end


text:
module Mod_NSFractionalStep
   use typre
   use Mod_Mesh
   use Mod_Memor

   use Mod_NavierStokes
   use Mod_ParallelSystemInterface
   private
   public NSFractionalStepProblem,NSFractionalStepProblem_Const

   type, extends(NavierStokesProblem) :: NSFractionalStepProblem

      !Linear System
      class(ParallelSystemInterface), pointer   :: LinearSystemP => NULL()
      real(rp), allocatable                     :: unknoP(:,:)
      integer(ip),  allocatable :: kfl_fixpr(:)       ! Nodal fixity for the pressure (NSF)
      real(rp), allocatable :: bvpress(:)

      integer(ip) :: kfl_restar = 0_ip                ! I need to know if this was a restart problem
      !to decide wether to initialize hydrostatic pressure field or not
      integer(ip) :: FractionalExternalIterations = 1

contains

      procedure :: SetExmod                           => nsf_SetExmod
      procedure :: LinearSystemMemall                 => nsf_LinearSystemMemall
      procedure :: LinearSystemTurnof                 => nsf_LinearSystemTurnof
      procedure :: SpecificReanut                     => nsf_reanut
      procedure :: SpecificReaMPI                     => nsf_reampi
      procedure :: SpecificCrankNicolsonEndste        => NULLSUB
      !procedure :: SpecificReabcsMPI                  => nsf_reabcsMPI
      procedure :: Solite                             => nsf_solite
      procedure :: Endite                             => nsf_endite
      procedure :: SpecificTurnof                     => nsf_turnof
      procedure :: InnerResiduals                     => nsf_InnerResiduals
      procedure ::  SpecificReabcs                    => nsf_Reabcs

      !For Initializing the hydrostatic pressure field
      !Need to note down if this is a initial restart problem
      procedure :: SpecificRestart                    => nsi_restar
      procedure :: SpecificBegste                     => nsf_begste
      procedure :: SpecificOuterr                     => nsf_outerr
      procedure :: SpecificRefine                     => nsf_refine
   end type

   interface
      subroutine nsf_solite(a)
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a

      end subroutine

      subroutine nsf_endite(a,itask)
         use typre
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a
         integer(ip)                    :: itask
      end subroutine

      subroutine nsf_InnerResiduals(a,rinsi,rprnsi)
         use typre
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a
         real(rp) :: rinsi,rprnsi
      end subroutine

      subroutine nsi_restar(a,itask)
         use typre
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a
         integer(ip), intent(in)                   :: itask
      end subroutine

      subroutine nsf_begste(a)
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a
      end subroutine

      subroutine nsf_refine(a,itask)
         use typre
         import NSFractionalStepProblem
         implicit none
         class(NSFractionalStepProblem) :: a
         character(6) :: itask
      end subroutine

   end interface

  interface NSFractionalStepProblem_Const
      procedure constructor
  end interface NSFractionalStepProblem_Const

  contains

    function constructor()
        class(NSFractionalStepProblem), pointer :: constructor

        allocate(constructor)

    end function constructor

   subroutine nsf_LinearSystemMemall(a)
      implicit none
      class(NSFractionalStepProblem) :: a

      integer(ip) :: npoin,ndime
      character(150) :: exmod,exfile
      character(150) :: auxstring

      exfile = 'nsi'

      !Velocity linear system
      exmod = 'nsf'
      call a%Mesh%GetNdime(ndime)
      auxstring = trim(a%InputFolder)//'/'//adjustl(trim(a%namda))//'.'//adjustl(trim(exfile))//'.sol '

      !Create ParallelSystem
      call a%ParallelLibrary%CreateSystem(a%LinearSystem, a%Memor)
      call a%LinearSystem%SetFlush(a%kfl_flush)
      call a%LinearSystem%SetMPI(a%MPIcomm,a%MPIsize,a%MPIroot,a%MPIrank)
      call a%Mesh%InitializeSystem(ndime,0,a%LinearSystem,a%Memor,auxstring,exmod,a%lun_solve)

      call a%Mesh%GetNpoin(npoin)
      call a%Memor%alloc(ndime,npoin,a%unkno,'unkno_nsf',trim(exmod)//'_LinearSystemMemall')

      !Pressure Linear System
      exmod = 'nsp'
      auxstring = trim(a%InputFolder)//'/'//adjustl(trim(a%namda))//'.'//adjustl(trim(exfile))//'.sol '
      call a%ParallelLibrary%CreateSystem(a%LinearSystemP, a%Memor)
      call a%LinearSystemP%SetFlush(a%kfl_flush)
      call a%LinearSystemP%SetMPI(a%MPIcomm,a%MPIsize,a%MPIroot,a%MPIrank)
      call a%Mesh%InitializeSystem(1,0,a%LinearSystemP,a%Memor,auxstring,exmod,a%lun_solve)

      call a%Mesh%GetNpoin(npoin)
      call a%Memor%alloc(1,npoin,a%unknoP,'unknoP',trim(exmod)//'_LinearSystemMemall')

   end subroutine

   subroutine nsf_LinearSystemTurnof(a)
      use typre
      use Mod_PhysicalProblem
      implicit none
      class(NSFractionalStepProblem) :: a

      integer(ip) :: npoin,ndime

      !Velocity Linear System
      call a%Mesh%GetNpoin(npoin)
      call a%Mesh%GetNdime(ndime)
      call a%Memor%dealloc(ndime,size(a%unkno,2),a%unkno,'unkno_nsf','nsf_LinearSystemTurnof')
      call a%LinearSystem%Deallocate
      call a%ParallelLibrary%DeallocateSystem(a%LinearSystem, a%Memor)
      !call a%Mesh%DeallocateSystem(a%LinearSystem,a%Memor)

      !Pressure Linear System
      call a%Memor%dealloc(1,size(a%unknoP,2),a%unknoP,'unknoP','nsf_LinearSystemTurnof')
      call a%LinearSystemP%Deallocate
      call a%ParallelLibrary%DeallocateSystem(a%LinearSystemP, a%Memor)
      !call a%Mesh%DeallocateSystem(a%LinearSystemP,a%Memor)

end subroutine

   subroutine nsf_SetExmod(a)
      use typre
      implicit none

      class(NSFractionalStepProblem) :: a

      a%exmod = 'nsi'
      a%namod = 'NSTINC_FractionalStep'

   end subroutine

   subroutine nsf_reanut(a,itask)
      use typre
      use MPI
      implicit none
      class(NSFractionalStepProblem) :: a
      integer(ip) :: itask
      interface
         subroutine nsi_reanut(a,itask)
            use typre
            use Mod_NavierStokes
            implicit none
            class(NavierStokesProblem) :: a
            integer(ip) :: itask
         end subroutine
      end interface

      !For Listener
      real(rp), pointer     :: param(:) => NULL()
      character(5), pointer :: words(:) => NULL()
      integer(ip), pointer  :: nnpar => NULL(),nnwor => NULL()
      integer(ip) :: ierr

      call nsi_reanut(a,itask)
      if (itask == 0) then
         a%kfl_repro_SkipFE = 1    !Default is skip Finite Element part of the residual
         a%FractionalExternalIterations = 1

      elseif (itask == 1) then
         call a%Listener%getarrs(words,param,nnpar,nnwor)
         if(words(1)=='EXTIT') then
            a%FractionalExternalIterations = param(1)
         endif

      elseif (itask == 100) then
         if (a%kfl_repro == 0) 	a%kfl_repro = 1
      endif
   end subroutine

   subroutine nsf_reampi(a)
      use MPI
      implicit none
      class(NSFractionalStepProblem) :: a
      interface
         subroutine nsi_reampi(a)
            use typre
            use Mod_NavierStokes
            implicit none
            class(NavierStokesProblem) :: a
         end subroutine
      end interface

      integer(ip) :: ierr

      call nsi_reampi(a)
      CALL MPI_BCAST(a%FractionalExternalIterations, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)
   end subroutine

   subroutine nsf_outerr(a)
      use typre
      implicit none
      class(NSFractionalStepProblem) :: a
      integer(ip) :: itask
      interface
         subroutine nsi_outerr(a)
            use typre
            use Mod_NavierStokes
            implicit none
            class(NavierStokesProblem) :: a
         end subroutine
      end interface

      call nsi_outerr(a)
      if(a%kfl_timei==0)then
         call runend('FRACTIONAL STEP CAN NOT BE USED IN NON TRANSIENT CASES')
      end if
      if(a%kfl_inist==1)then
         call runend('INISTOKES NOT READY IN FRACTIONAL STEP')
      end if
   end subroutine

   subroutine nsf_Reabcs(a,itask,kflag)
      implicit none
      class(NSFractionalStepProblem) :: a
      integer(ip) :: itask
      integer(ip), optional :: kflag

      integer(ip) :: npoin
      interface
         subroutine nsi_Reabcs(a,itask,kflag)
            use typre
            use Mod_NavierStokes
            implicit none
            class(NavierStokesProblem) :: a
            integer(ip) :: itask
            integer(ip), optional :: kflag
         end subroutine

         subroutine nsf_pbcpre(a)
            use typre
            import NSFractionalStepProblem
            implicit none
            class(NSFractionalStepProblem) :: a
         end subroutine
      end interface

      call nsi_Reabcs(a,itask,kflag)
      if (itask == 100) then
         call a%Mesh%GetNpoin(npoin)
         call a%Memor%alloc(npoin,a%kfl_fixpr,'kfl_fixpr','nsf_reabcsMPI')
         call a%Memor%alloc(npoin,a%bvpress,'bvpress','nsf_reabcsMPI')
         call nsf_pbcpre(a)
      endif

   end subroutine

   subroutine nsf_turnof(a)
      implicit none
      class (NSFractionalStepProblem) :: a

      integer(ip) :: npoin
      interface
         subroutine nsi_turnof(a)
            use typre
            use Mod_NavierStokes
            implicit none
            class(NavierStokesProblem) :: a
         end subroutine
      end interface

      call nsi_turnof(a)
      call a%Mesh%GetNpoin(npoin)
      call a%Memor%dealloc(size(a%kfl_fixpr),a%kfl_fixpr,'kfl_fixpr','nsf_reabcsMPI')
      call a%Memor%dealloc(size(a%bvpress),a%bvpress,'bvpress','nsf_reabcsMPI')

   end subroutine

   subroutine NULLSUB(a)
      use typre
      implicit none
      class(NSFractionalStepProblem) :: a
   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
! Subroutine getomi finds poles and weights for representing epsInv.
! Written by J. Kas, 2006
      SUBROUTINE getomi(wgrid, epsInv, NPts, NPoles, omi, gi, Deltai,
     &     eps0)
      IMPLICIT NONE
!     Input variables
!     wGrid(NPts)  - grid of omega
!     epsInv(NPts) - loss function
!     NPoles       - Number of poles to represent epsInv

!     Output
!     omi          - omega position of poles
!     gi           - pole strengths
!     Deltai       - energry region about omi where
      INTEGER NPoles, NPts
      DOUBLE PRECISION wGrid(NPts), epsInv(NPts), omi(NPoles),
     &     gi(NPoles), Deltai(NPoles), eps0

!     Local variables
      INTEGER maxPts, i1, i2
      PARAMETER (maxPts = 50000)
      DOUBLE PRECISION M0, M0i, M1i, MM1i, M1Old, M1New, omNew, omOld,
     $       pi, Tmp1(maxPts), Tmp2(maxPts), IntGrl, a, b, c, frac

c     INTEGER MaxPoles, iMax,
c     DOUBLE PRECISION egap
      PARAMETER (pi = 3.14159 26535 89793 23846 26433d0)
      EXTERNAL IntGrl
!     Find zeroth moment
!      frac = (1.d0-0.01d0)
      frac = 0.999d0
      M0 = 0.5d0*epsInv(1)*wGrid(1)
      DO i1 = 2, NPts
         M0 = M0 + 0.5d0*(epsInv(i1)+epsInv(i1-1))*
     &        (wGrid(i1)-wGrid(i1-1))
         Tmp1(i1) = epsInv(i1)/wGrid(i1)
         Tmp2(i1) = epsInv(i1)*wGrid(i1)
      END DO
!     Find gap energy
!      DO i1 = 1, NPts

!     Now find energy regions such that the integral of epsinv in each region is
!     equal.
      omOld = 0.d0
      omNew = 0.d0
      M1Old = 0.d0
      M1New = 0.d0
      b     = 0.d0
      ! Save the first 10 poles for low energies.
      DO i1 = 1, NPoles
         M0i = 0.5d0*epsInv(1)*wGrid(1)
         M1New = 0.5d0*epsInv(1)*wGrid(1)**2
         DO i2 = 2, NPts
            IF(M0i.ge.(DBLE(i1)*frac*M0/DBLE(NPoles))) THEN
               omOld = omNew
               omNew = wGrid(i2-1)
               IF(i2.eq.1) THEN
                  M0i = 0.d0
               ELSEIF(i2.eq.2) THEN
                  M0i = 0.5d0*epsInv(1)*wGrid(1)
               ELSE
                  M0i = M0i - 0.5d0*(epsInv(i2-1)+ epsInv(i2-2))*
     $              (wGrid(i2-1)-wGrid(i2-2))
               END IF
               GOTO 5
            ELSE
               M0i = M0i + 0.5d0*(epsInv(i2)+epsInv(i2-1))*
     $              (wGrid(i2)-wGrid(i2-1))
               M1New = M1New + 0.5d0*(epsInv(i2)*wGrid(i2)+ epsInv(i2-1)
     $              *wGrid(i2-1))*(wGrid(i2)-wGrid(i2-1))
            END IF
         END DO
 5       CONTINUE

!        The region has been found up to the wGrid accuracy. Now interpolate
!        to find omi
         IF(i2.eq.2) THEN
            a = epsInv(i2-1)/omNew
         ELSE
            a = (epsInv(i2-1) - epsInv(i2-2))/(omNew - wGrid(i2-2))
            b = epsInv(i2-2)
         END IF
         a = MAX(a, 1.d-6)
         c = 2.d0/a*(M0i - frac*M0*(DBLE(i1)/DBLE(NPoles)))

         omNew = wGrid(i2-2) - b/a + SQRT((b/a)**2 - c)
         IF((omNew.lt.wGrid(i2-2)).or.(omNew.gt.wGrid(i2-1))) THEN
            omNew = wGrid(i2-2) - b/a - SQRT((b/a)**2 - c)
         END IF

         IF(.FALSE.) THEN
            omOld = DBLE(i1-1)*500.d0/DBLE(NPoles)
            omNew = DBLE(i1)*500.d0/DBLE(NPoles)
         END IF
!        Now set omi and g s.t. the inverse and first moments are preserved for
!        this region.
         MM1i = 0.d0
         M1i  = 0.d0
         IF(omOld.eq.0.d0) THEN
            MM1i = 0.5d0*epsInv(1)
            M1i  = 0.5d0*epsInv(1)*wGrid(1)**2
            omOld = wGrid(1)
         END IF
         MM1i = MM1i + IntGrl(omOld, omNew, wGrid, Tmp1, NPts)
         M1i  = M1i + IntGrl(omOld, omNew, wGrid, Tmp2, NPts)
         gi(i1)  = 2.d0/pi*MM1i
         omi(i1) = SQRT(M1i/MM1i)
         Deltai(i1) = omNew-omOld
         M1Old  = M1New
         !PRINT '(a27,i6,20f20.10)', 'i1, omi, omOld, omNew, a, b', i1,
         !&    omi(i1), omOld, omNew, a, b
      END DO

!     Shift and scale the poles to set inverse moment to 1-1/eps0, without changing
!     first moment. (eps0 is input by user)
      CALL getdom(gi,omi,NPoles,eps0)

      RETURN
      END

      DOUBLE PRECISION FUNCTION IntGrl(a, b, x, y, n)
      INTEGER n, i, iStart, iLast
      DOUBLE PRECISION x(n), y(n), dx, a, b,
     &     ya, yb
c      DOUBLE PRECISION aEnd, bEnd,

      iStart = -1
      iLast = -1
      IntGrl = 0.d0
      dx = x(2) - x(1)

      DO i = 1, n
         IF(x(i).ge.a.and.x(i).lt.b) THEN
            iStart = i
            goto 5
         ELSE
            iStart = -1
         END IF
      END DO
 5    CONTINUE

      DO i = n, 1, -1
         IF(x(i).le.b.and.x(i).gt.a) THEN
            iLast = i
            goto 10
         ELSE
            iLast = -1
         END IF
      END DO
 10   CONTINUE

      IF (iStart.lt.0.and.iLast.lt.0) THEN
         CALL terp(x, y, n, 1, a, ya)
         CALL terp(x, y, n, 1, b, yb)
         IntGrl = IntGrl + 0.5d0*(yb+ya)*(b-a)
      ELSE
         DO i = iStart, iLast-1
            dx=x(i+1)-x(i)
            IntGrl = IntGrl + 0.5d0*(y(i)+y(i+1))*dx
         END DO

         CALL terp(x, y, n, 1, a, ya)
         CALL terp(x, y, n, 1, b, yb)

         IntGrl = IntGrl + 0.5d0*(ya+y(iStart))*(x(iStart) - a)
         IntGrl = IntGrl + 0.5d0*(yb+y(iLast))*(b-x(iLast))
      END IF
      RETURN
      END

! This will change the pole representation of the loss function without
! changing the f-sum rules. Written by J. Kas 2006
      SUBROUTINE getdom(gi,omi,NPoles,eps0)
      IMPLICIT NONE

!     Input/output
      INTEGER NPoles
      DOUBLE PRECISION gi(NPoles), omi(NPoles), eps0

!     Local
      INTEGER i1, i2, MxIter
      DOUBLE PRECISION MM1, tol, dom, dom1, dom2, sumg1, sumg2
      PARAMETER(MxIter = 10000, tol = 1.d-6)
      IF(eps0.lt.-1.5d0) RETURN
      IF(eps0.lt.0.d0) THEN
         MM1 = 1.d0
      ELSE
         MM1   = 1.d0 - 1.d0/eps0
      END IF
      dom1  = -omi(1)*0.999d0
      dom2  = omi(1)*0.999d0
!     Bracket dom between dom1 and dom2
      DO i1 = 1, MxIter
         sumg1 = 0.d0
         sumg2 = 0.d0
         DO i2 = 1, NPoles
            sumg1 = sumg1 + gi(i2)*(omi(i2)/(omi(i2)-dom1))**2
            sumg2 = sumg2 + gi(i2)*(omi(i2)/(omi(i2)-dom2))**2
         END DO
         IF((sumg1.lt.MM1).and.(sumg2.gt.MM1)) GOTO 5
         !IF(sumg1.gt.MM1) dom1 = dom1*2.d0
         !IF(sumg2.lt.MM1) dom2 = dom2/2.d0
         PRINT*, 'Error in getdom.f.'
         PRINT*, 'sumg1, sumg2', i1, sumg1, sumg2, MM1, dom1, dom2
         STOP
      END DO
 5    CONTINUE                  ! dom bracketed

      dom = 0.5d0*(dom1+dom2)
      DO i1 = 1, MxIter
         sumg1 = 0.d0
         DO i2 = 1, NPoles
            sumg1 = sumg1 + gi(i2)*(omi(i2)/(omi(i2)-dom))**2
         END DO
         IF(abs((sumg1-MM1)/MM1).lt.tol) GOTO 10
         IF((sumg1-MM1).gt.0.d0) THEN
            dom2 = dom
         ELSE
            dom1 = dom
         END IF
         dom = 0.5d0*(dom1+dom2)
      END DO
 10   CONTINUE

      DO i1 = 1, NPoles
         gi(i1)  = gi(i1)*(omi(i1)/(omi(i1)-dom))**2
         omi(i1) = omi(i1)-dom
      END DO
      RETURN
      END