text:
      SUBROUTINE GFLUXV(DTDEL,TDEL,TAUCUMIN,WDEL,CDEL,UBAR0,F0PI,RSF,  &
                        BTOP,BSURF,FMIDP,FMIDM,DIFFV,FLUXUP,FLUXDN,    &
                        detau)

!     Legacy Mars GCM v24
!     Mars Climate Modeling Center
!     NASA Ames Research Center

!  THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS
!  FOR THE VISIBLE  FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
!  THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
!  MEASURED FROM THE TOP OF EACH LAYER.  (DTAU) TOP OF EACH LAYER HAS
!  OPTICAL DEPTH TAU(N).IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
!  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
!  FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
! THIS SUBROUTINE DIFFERS FROM ITS IR COUNTERPART IN THAT HERE WE SOLVE FOR
! THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR
! J.A.S., 37, 630-642, 1980.
! THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
! VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
!
! NLL           = NUMBER OF LEVELS (NAYER + 1) THAT WILL BE SOLVED
! NAYER         = NUMBER OF LAYERS (NOTE DIFFERENT SPELLING HERE)
! WAVEN         = WAVELENGTH FOR THE COMPUTATION
! DTDEL(NLAYER) = ARRAY OPTICAL DEPTH OF THE LAYERS
! TDEL(NLL)     = ARRAY COLUMN OPTICAL DEPTH AT THE LEVELS
! WDEL(NLEVEL)  = SINGLE SCATTERING ALBEDO
! CDEL(NLL)     = ASYMMETRY FACTORS, 0=ISOTROPIC
! UBARV         = AVERAGE ANGLE,
! UBAR0         = SOLAR ZENITH ANGLE
! F0PI          = INCIDENT SOLAR DIRECT BEAM FLUX
! RSF           = SURFACE REFLECTANCE
! BTOP          = UPPER BOUNDARY CONDITION ON DIFFUSE FLUX
! BSURF         = REFLECTED DIRECT BEAM = (1-RSFI)*F0PI*EDP-TAU/U
! FP(NLEVEL)    = UPWARD FLUX AT LEVELS
! FM(NLEVEL)    = DOWNWARD FLUX AT LEVELS
! FMIDP(NLAYER) = UPWARD FLUX AT LAYER MIDPOINTS
! FMIDM(NLAYER) = DOWNWARD FLUX AT LAYER MIDPOINTS
! added Dec 2002
! DIFFV         = downward diffuse solar flux at the surface
!
!======================================================================!

      use grid_h
      use radinc_h

      implicit none

      integer, parameter :: NLP = 101     ! MUST BE LARGER THAN NLEVEL

      REAL*8 EM, EP
      REAL*8 W0(L_NLAYRAD), COSBAR(L_NLAYRAD), DTAU(L_NLAYRAD)
      REAL*8 TAU(L_NLEVRAD), WDEL(L_NLAYRAD), CDEL(L_NLAYRAD)
      REAL*8 DTDEL(L_NLAYRAD), TDEL(L_NLEVRAD)
      REAL*8 FMIDP(L_NLAYRAD), FMIDM(L_NLAYRAD)
      REAL*8 LAMDA(NLP), ALPHA(NLP), XK1(NLP), XK2(NLP)
      REAL*8 G1(NLP), G2(NLP), G3(NLP), GAMA(NLP), CP(NLP), CM(NLP)
      REAL*8 CPM1(NLP)
      REAL*8 CMM1(NLP), E1(NLP), E2(NLP), E3(NLP), E4(NLP), EXPTRM(NLP)
      REAL*8 FLUXUP, FLUXDN
      REAL*8 FACTOR, TAUCUMIN(L_LEVELS), TAUCUM(L_LEVELS)
      real*8  :: detau

      integer NAYER, L, K
      real*8  ubar0, f0pi, rsf, btop, bsurf, g4, denom, am, ap
      real*8  taumax, taumid, cpmid, cmmid
      real*8  diffv

!======================================================================C

      NAYER  = L_NLAYRAD
      TAUMAX = L_TAUMAX    !Default is 35.0

!  Delta-Eddington Scaling

      FACTOR    = 1.0D0 - WDEL(1)*CDEL(1)**2

      TAU(1)    = TDEL(1)*FACTOR
      TAUCUM(1) = 0.0D0
      TAUCUM(2) = TAUCUMIN(2)*FACTOR
      TAUCUM(3) = TAUCUM(2) +(TAUCUMIN(3)-TAUCUMIN(2))*FACTOR

      DO L=1,L_NLAYRAD-1
        FACTOR      = 1.0D0 - WDEL(L)*CDEL(L)**2
        W0(L)       = WDEL(L)*(1.0D0-CDEL(L)**2)/FACTOR
        COSBAR(L)   = CDEL(L)/(1.0D0+CDEL(L))
        DTAU(L)     = DTDEL(L)*FACTOR
        TAU(L+1)    = TAU(L)+DTAU(L)
        K           = 2*(L+1)
        TAUCUM(K)   = TAU(L+1)
        TAUCUM(K+1) = TAUCUM(K) + (TAUCUMIN(K+1)-TAUCUMIN(K))*FACTOR
      END DO

!  Bottom layer

      L             = L_NLAYRAD
      FACTOR        = 1.0D0 - WDEL(L)*CDEL(L)**2
      W0(L)         = WDEL(L)*(1.0D0-CDEL(L)**2)/FACTOR
      COSBAR(L)     = CDEL(L)/(1.0D0+CDEL(L))
      DTAU(L)       = DTDEL(L)*FACTOR
      TAU(L+1)      = TAU(L)+DTAU(L)
      TAUCUM(2*L+1) = TAU(L+1)
      detau         = TAUCUM(2*L+1)

!     WE GO WITH THE QUADRATURE APPROACH HERE.  THE "SQRT(3)" factors
!     ARE THE UBARV TERM.

      DO L=1,L_NLAYRAD

        ALPHA(L)=SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )

!       SET OF CONSTANTS DETERMINED BY DOM

        G1(L)    = (SQRT(3.0)*0.5)*(2.0- W0(L)*(1.0+COSBAR(L)))
        G2(L)    = (SQRT(3.0)*W0(L)*0.5)*(1.0-COSBAR(L))
        G3(L)    = 0.5*(1.0-SQRT(3.0)*COSBAR(L)*UBAR0)
        LAMDA(L) = SQRT(G1(L)**2 - G2(L)**2)
        GAMA(L)  = (G1(L)-LAMDA(L))/G2(L)
      END DO

      DO L=1,L_NLAYRAD
        G4    = 1.0-G3(L)
        DENOM = LAMDA(L)**2 - 1./UBAR0**2

!       THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
!       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
!       THE SCATTERING GOES TO ZERO
!       PREVENT THIS WITH AN IF STATEMENT

        IF ( DENOM .EQ. 0.) THEN
          DENOM=1.E-10
        END IF

        AM = F0PI*W0(L)*(G4   *(G1(L)+1./UBAR0) +G2(L)*G3(L) )/DENOM
        AP = F0PI*W0(L)*(G3(L)*(G1(L)-1./UBAR0) +G2(L)*G4    )/DENOM

!       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
!       AT THE TOP OF THE LAYER, THAT IS LOWER   OPTICAL DEPTH TAU(L)

        CPM1(L) = AP*EXP(-MIN(TAU(L)/UBAR0,MAXEXP))
        CMM1(L) = AM*EXP(-MIN(TAU(L)/UBAR0,MAXEXP))

!       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
!       BOTTOM OF THE LAYER.  THAT IS AT HIGHER OPTICAL DEPTH TAU(L+1)

        CP(L) = AP*EXP(-MIN(TAU(L+1)/UBAR0,MAXEXP))
        CM(L) = AM*EXP(-MIN(TAU(L+1)/UBAR0,MAXEXP))

      END DO

!     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
!     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD

      DO L=1,L_NLAYRAD
        EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*DTAU(L))  ! CLIPPED EXPONENTIAL
        EP = EXP(EXPTRM(L))

        EM        = 1.0/EP
        E1(L)     = EP+GAMA(L)*EM
        E2(L)     = EP-GAMA(L)*EM
        E3(L)     = GAMA(L)*EP+EM
        E4(L)     = GAMA(L)*EP-EM
      END DO

      CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,        &
                   BSURF,RSF,XK1,XK2)

!     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.

      DO L=1,L_NLAYRAD-1
        EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*(TAUCUM(2*L+1)-TAUCUM(2*L)))

        EP = EXP(EXPTRM(L))

        EM    = 1.0/EP
        G4    = 1.0-G3(L)
        DENOM = LAMDA(L)**2 - 1./UBAR0**2

!       THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
!       THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
!       THE SCATTERING GOES TO ZERO
!       PREVENT THIS WITH A IF STATEMENT

        IF ( DENOM .EQ. 0.) THEN
          DENOM=1.E-10
        END IF

        AM = F0PI*W0(L)*(G4   *(G1(L)+1./UBAR0) +G2(L)*G3(L) )/DENOM
        AP = F0PI*W0(L)*(G3(L)*(G1(L)-1./UBAR0) +G2(L)*G4    )/DENOM

!       CPMID AND CMMID  ARE THE CPLUS AND CMINUS TERMS EVALUATED
!       AT THE MIDDLE OF THE LAYER.

        TAUMID   = TAUCUM(2*L+1)

        CPMID = AP*EXP(-MIN(TAUMID/UBAR0,MAXEXP))
        CMMID = AM*EXP(-MIN(TAUMID/UBAR0,MAXEXP))

        FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
        FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID

!       ADD THE DIRECT FLUX TO THE DOWNWELLING TERM

        FMIDM(L)= FMIDM(L)+UBAR0*F0PI*EXP(-MIN(TAUMID/UBAR0,MAXEXP))

      END DO

!     FLUX AT THE top layer

      EP    = 1.0
      EM    = 1.0
      G4    = 1.0-G3(1)
      DENOM = LAMDA(1)**2 - 1./UBAR0**2

!     THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
!     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
!     THE SCATTERING GOES TO ZERO
!     PREVENT THIS WITH A IF STATEMENT

      IF ( DENOM .EQ. 0.) THEN
        DENOM=1.E-10
      END IF

      AM = F0PI*W0(1)*(G4   *(G1(1)+1./UBAR0) +G2(1)*G3(1) )/DENOM
      AP = F0PI*W0(1)*(G3(1)*(G1(1)-1./UBAR0) +G2(1)*G4    )/DENOM

!     CPMID AND CMMID  ARE THE CPLUS AND CMINUS TERMS EVALUATED
!     AT THE MIDDLE OF THE LAYER.

      CPMID  = AP
      CMMID  = AM

      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
      FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID

!     ADD THE DIRECT FLUX TO THE DOWNWELLING TERM

      fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1)/UBAR0,MAXEXP))

!     This is for the "special" bottom layer, where we take
!     DTAU instead of DTAU/2.

      L     = L_NLAYRAD
      EXPTRM(L) = MIN(TAUMAX,LAMDA(L)*(TAUCUM(L_LEVELS)-               &
                                       TAUCUM(L_LEVELS-1)))

      EP    = EXP(EXPTRM(L))
      EM    = 1.0/EP
      G4    = 1.0-G3(L)
      DENOM = LAMDA(L)**2 - 1./UBAR0**2

!     THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
!     THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
!     THE SCATTERING GOES TO ZERO
!     PREVENT THIS WITH A IF STATEMENT

      IF ( DENOM .EQ. 0.) THEN
        DENOM=1.E-10
      END IF

      AM = F0PI*W0(L)*(G4   *(G1(L)+1./UBAR0) +G2(L)*G3(L) )/DENOM
      AP = F0PI*W0(L)*(G3(L)*(G1(L)-1./UBAR0) +G2(L)*G4    )/DENOM

!     CPMID AND CMMID  ARE THE CPLUS AND CMINUS TERMS EVALUATED
!     AT THE MIDDLE OF THE LAYER.

      TAUMID   = MIN(TAUCUM(L_LEVELS),TAUMAX)
      CPMID    = AP*EXP(-MIN(TAUMID/UBAR0,MAXEXP))
      CMMID    = AM*EXP(-MIN(TAUMID/UBAR0,MAXEXP))

      FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
      FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID

!  Save the diffuse downward flux for TEMPGR calculations

      DIFFV = FMIDM(L)

!     ADD THE DIRECT FLUX TO THE DOWNWELLING TERM

      FMIDM(L)= FMIDM(L)+UBAR0*F0PI*EXP(-MIN(TAUMID/UBAR0,MAXEXP))

      RETURN
      end subroutine gfluxv

max_stars_repo_name: nasa/legacy-mars-global-climate-model

max_issues_repo_name: nasa/legacy-mars-global-climate-model

max_forks_repo_name: nasa/legacy-mars-global-climate-model


text:
      subroutine sloss(idout,dslost,wmelt,nowcrp)
c*******************************************************************
c                                                                  *
c   This subroutine called from SR CONTIN calculates sediment      *
c   concentration and sediment yield on a storm-by-storm basis,    *
c   and prints out soil loss and sediment load information by      *
c   input segment and by erosion /deposition section. It also      *
c   prints out an abbreviated hydrology output for the event       *
c   under consideration, by calling SR HYDOUT.                     *
c                                                                  *
c*******************************************************************
c*******************************************************************
c                                                                  *
c   Arguments                                                      *
c     idout - flag for hydrology output                            *
c     dslost - net soil loss/gain at each point on hillslope for   *
c              a storm event (time*dG/dx) - (kg/m**2)              *
c                                                                  *
c*******************************************************************
c                                                                  *
      include 'pmxcrp.inc'
      include 'pmxnsl.inc'
      include 'pmxpln.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
      include 'pmxprt.inc'
      include 'pmxelm.inc'
      include 'pmxpts.inc'
      include 'pmxseg.inc'
      include 'pmxpnd.inc'
c
c*******************************************************************
c                                                                  *
c    Common Blocks                                                 *
c                                                                  *
c*******************************************************************
c
      include 'cavloss.inc'
c
c*******************************************************************
c                                                                  *
c    avloss variables updated                                      *
c       dsmon(mxplan,100),dsyear(mxplan,100),dsavg(mxplan,100)     *
c       avsols,avsole,avsolm,avsoly,avsolf                         *
c                                                                  *
c*******************************************************************
c
      include 'ccntour.inc'
c
      include 'cconsta.inc'
c
      include 'ccover.inc'
c
      include 'ccrpout.inc'
c
      include 'cdist.inc'
c
      include 'cefflen.inc'
c     read: efflen(pxplan)
      include 'cends.inc'
c
c*******************************************************************
c                                                                  *
c    ends variables updated                                        *
c       qsout                                                      *
c                                                                  *
c*******************************************************************
c
      include 'cenrpas.inc'
c
c*******************************************************************
c                                                                  *
c    enrpas variables updated                                      *
c       enrmm1,enrmm2,enryy1,enryy2                                *
c       enrff1,enrff2                                              *
c       frcmm1(10),frcmm2(10)                                      *
c       frcyy1(10),frcyy2(10),frcff1(10),frcff2(10)                *
c                                                                  *
c*******************************************************************
c
      include 'cerdval.inc'
c
c*******************************************************************
c                                                                  *
c  erdval variables updated                                        *
c      load(101)                                                   *
c                                                                  *
c*******************************************************************
c
      include 'chydrol.inc'
c
      include 'ciravlo.inc'
c
c**********************************************************************
c                                                                     *
c  iravlo variables updated                                           *
c     irsold,irsolm, irsoly, irsolt                               *
c                                                                     *
c**********************************************************************
c
      include 'cirriga.inc'
c
      include 'cparame.inc'
c
      include 'cpart.inc'
c
      include 'cparval.inc'
c
      include 'cseddet.inc'
      include 'cslpopt.inc'
      include 'cstore.inc'
      include 'cstruc.inc'
      include 'cupdate.inc'
c
      include 'cwater.inc'
c
      integer idout, nowcrp
      real wmelt
c
c*******************************************************************
c                                                                  *
c    local variables                                               *
c     conc   - sediment concentration at OFE end - (kg/kg)         *
c     dslend - dimensional sediment load at OFE end - (kg/m)       *
c     dslod1 - dimensional sediment load at point - (kg/m)         *
c     dslod2 - dimensional sediment load at point - (kg/m)         *
c     dslost - net soil loss/gain at each point on hillslope for   *
c              a storm event (time*dG/dx) - (kg/m**2)              *
c     csedls - sediment loss from contoured OFEs in (kg)           *
c                                                                  *
c*******************************************************************
c
      save
c     real conc(mxplan)
      integer i, ipart, j, jun
      real dslend, dslod1, dslod2, sumslc
      real dslost(mxplan,100),csedls(mxplan)
      data csedls /mxplan * 0.0/
c
c*******************************************************************
c
      jun = 31
      if (imodel.eq.2) jun = 32
c
c     initialize load variables
c
      dslend = 0.0
      dslod1 = 0.0
      dslod2 = 0.0
c     conc(iplane)=0.0
      do 10 j = 2, 101
        dslost(iplane,j-1) = 0.0
   10 continue

c
c     changed 1-11-2010 - always initialize csedls, don't use save from above - jrf
c
      csedls(iplane) = 0.0
c
      if (ioutpt.eq.1.and.imodel.eq.1)
     1     call hydout(jun,idout,wmelt,nowcrp)
c
c     Calculate dimensional sediment load at beginning of OFE (x=0)
c
      dslod1 = load(1) * effdrn(iplane) * tcend * width(iplane) /
     1    rspace(iplane)
c
c     If contours present, then compute an average detachment rate on
c     the OFE using the sediment discharge at the end of the contours.
c     contours are assumed not to fail if present.
c
c     START OF CONTOUR ELSE-IF
c
      ofelod(iplane)=0.0
      if (conseq(nowcrp,iplane).ne.0) then
        dslend = load(101) * effdrn(iplane) * tcend * width(iplane) /
     1      rspace(iplane)
        avsols = dslend / cntlen(iplane)
c
        do 20 j = 2, 101
          dslost(iplane,j-1) = avsols
   20   continue

c
c       following line added 2/4/91 to calculate interrill contribution
c       for an OFE on which contours have held
c
        irdgdx(iplane) = (theta*tcend*effdrn(iplane)*width(iplane)) / (
     1      rspace(iplane)*cntlen(iplane))
c
c       6-19-2002 Need to report the sediment leaving the sides
c       of each contoured OFE in the output file, as well as the
c       accumulated sediment leaving the sides of the OFEs. Add
c       new variable "avsolc(iplane)" to track the storm sediment
c       loss from the contoured OFE sides.   dcf

        avsolc(iplane) = dslend

c       Compute the total sediment leaving the contoured OFE to the
c       side, by finding an equivalent width (OFE area/contour lgth)
        csedls(iplane) = avsolc(iplane)*
     1                   (slplen(iplane)*fwidth(1)/cntlen(iplane))

c       Set the average sediment discharge off the plane equal to zero
c       since if contours hold no sediment reaches the next OFE from
c       above.
c
        avsols = 0.0
c
c     ELSE if no contours are present on the OFE
c     calculate the net sediment loss or gain at a point by taking
c     the difference between the sediment load at the points divided
c     by the distance between the points.
c
      else
c
        do 30 j = 2, 101
          dslod2 = load(j) * effdrn(iplane) * tcend * width(iplane) /
     1        rspace(iplane)
          dslost(iplane,j-1) = (dslod2-dslod1) / (slplen(iplane)*0.01)
          dslod1 = dslod2
   30   continue
        ofelod(iplane)=dslod1
c
c       ADDED 11/8/90 by dcf - to allow estimation of interrill
c       contribution to sediment loss.
c
        irdgdx(iplane) = (theta*tcend*effdrn(iplane)*width(iplane)) / (
     1      rspace(iplane)*slplen(iplane))
c
c
        if (cntlen(iplane).lt.0.01) cntlen(iplane) = slplen(iplane)
        avsols = dslod2 / cntlen(iplane)
c
c     END OF CONTOUR ELSE-IF
c
      end if
c
c     summing for monthly and annual dslost
c
      do 40 j = 2, 101
        dsmon(iplane,j-1) = dsmon(iplane,j-1) + dslost(iplane,j-1)
        dsyear(iplane,j-1) = dsyear(iplane,j-1) + dslost(iplane,j-1)
        dsavg(iplane,j-1) = dsavg(iplane,j-1) + dslost(iplane,j-1)
   40 continue
c
c     START OF IPLANE=NPLANE ELSE-IF  (this the last
c     OFE on the current hillslope)
c
      if (iplane.eq.nplane) then
c
c       If contours are not present the sediment reaching the end of
c       the hillslope is equal to the sediment load at the last point
c       of the last OFE on the current hillslope.
c       else - if contours are present then they are assumed to not
c       fail and all sediment is assumed routed to the side of the
c       hillslope and none reaches the end of the OFE.
c       6/26/2002 - a new calculation for "avsole" is made for the
c       case of a contoured hillslope - that includes all sediment
c       leaving the profile in any direction.
c
        sumslc = 0.0
        do 45 j = 1,nplane
          sumslc = sumslc + csedls(j)
  45    continue
        if (qout.gt.0.0 .and. conseq(nowcrp,nplane).eq.0)then
          avsole = (sumslc + dslod2*fwidth(1))/fwidth(1)
        else
          avsole = sumslc/fwidth(1)
        endif
c
        if (avsole.lt.0.0) avsole = 0.0
c
        avsolm = avsolm + avsole
        avsoly = avsoly + avsole
        avsolf = avsolf + avsole
c
        irsold = avsole * irfrac
        irsolm = irsolm + irsold
        irsoly = irsoly + irsold
        irsolt = irsolt + irsold
c
c       summing for weighted enrichment ratio
c       sets single storm and sums for monthly and annual
c
        enrmm1 = enrmm1 + (enrato(nplane)*avsole)
        enrmm2 = enrmm2 + avsole
        enryy1 = enryy1 + (enrato(nplane)*avsole)
        enryy2 = enryy2 + avsole
        enrff1 = enrff1 + (enrato(nplane)*avsole)
        enrff2 = enrff2 + avsole
c
        do 50 i = 1, npart
c
          frcmm1(i) = frcmm1(i) + (frcflw(i,nplane)*avsole)
          frcmm2(i) = frcmm2(i) + avsole
          frcyy1(i) = frcyy1(i) + (frcflw(i,nplane)*avsole)
          frcyy2(i) = frcyy2(i) + avsole
          frcff1(i) = frcff1(i) + (frcflw(i,nplane)*avsole)
          frcff2(i) = frcff2(i) + avsole
c
   50   continue
c
c       compute sediment concentration (kg/m**3),
c       for each particle size leaving the hillslope,
c       avsole is in (kg/m) and runvol is in (m**3)
c
        if (iplane.eq.nplane) then
          do 60 ipart = 1, npart
            if (peakro(nplane).le.0.0) then
              sedcon(ipart,nplane) = 0.0
            else
              sedcon(ipart,nplane) = avsole / (runoff(nplane)*
     1            efflen(nplane)) * frcflw(ipart,nplane)
            end if
   60     continue
        end if
c
c
c     END OF IPLANE=NPLANE IF-ELSE
c
      end if
c
c
c     If contours hold, qsout equals zero since no sediment reaches
c     the next OFE.  If contours fail, qsout equals the sediment load
c     at the end of the current OFE.
c     if contours exist they are assumed not to fail
c
      if (conseq(nowcrp,iplane).eq.0) then
        if (qout.gt.0.0) then
          qsout = dslod2 / effdrn(iplane)
        else
          qsout = 0.0
        end if
      else
        qsout = 0.0
      end if
c
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
!
! File: amg_d_gs_solver_mod.f90
!
! Module: amg_d_gs_solver_mod
!
!  This module defines:
!  - the amg_d_gs_solver_type data structure containing the ingredients
!    for a local Gauss-Seidel iteration. We provide Forward GS (FWGS) and
!    backward GS (BWGS). The iterations are local to a process (they operate
!    on the block diagonal). Combined with a Jacobi smoother will generate a
!    hybrid-Gauss-Seidel solver, i.e. Gauss-Seidel within each process, Jacobi
!    among the processes.
!    With two objects as pre- and post-smoothers it is possible to build a
!    Forward-Backward smoother, suitable for symmetric iterations.
!
module amg_d_gs_solver

  use amg_d_base_solver_mod

  type, extends(amg_d_base_solver_type) :: amg_d_gs_solver_type
    type(psb_dspmat_type)      :: l, u
    integer(psb_ipk_)          :: sweeps
    real(psb_dpk_)             :: eps
  contains
    procedure, pass(sv) :: dump    => amg_d_gs_solver_dmp
    procedure, pass(sv) :: check   => d_gs_solver_check
    procedure, pass(sv) :: clone   => amg_d_gs_solver_clone
    procedure, pass(sv) :: clone_settings => amg_d_gs_solver_clone_settings
    procedure, pass(sv) :: clear_data     => amg_d_gs_solver_clear_data
    procedure, pass(sv) :: build   => amg_d_gs_solver_bld
    procedure, pass(sv) :: cnv     => amg_d_gs_solver_cnv
    procedure, pass(sv) :: apply_v => amg_d_gs_solver_apply_vect
    procedure, pass(sv) :: apply_a => amg_d_gs_solver_apply
    procedure, pass(sv) :: free    => d_gs_solver_free
    procedure, pass(sv) :: cseti   => d_gs_solver_cseti
    procedure, pass(sv) :: csetc   => d_gs_solver_csetc
    procedure, pass(sv) :: csetr   => d_gs_solver_csetr
    procedure, pass(sv) :: descr   => d_gs_solver_descr
    procedure, pass(sv) :: default => d_gs_solver_default
    procedure, pass(sv) :: sizeof  => d_gs_solver_sizeof
    procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros
    procedure, nopass   :: get_wrksz => d_gs_solver_get_wrksize
    procedure, nopass   :: get_fmt    => d_gs_solver_get_fmt
    procedure, nopass   :: get_id    => d_gs_solver_get_id
    procedure, nopass   :: is_iterative => d_gs_solver_is_iterative
  end type amg_d_gs_solver_type

  type, extends(amg_d_gs_solver_type) :: amg_d_bwgs_solver_type
  contains
    procedure, pass(sv) :: build    => amg_d_bwgs_solver_bld
    procedure, pass(sv) :: apply_v  => amg_d_bwgs_solver_apply_vect
    procedure, pass(sv) :: apply_a  => amg_d_bwgs_solver_apply
    procedure, nopass   :: get_fmt  => d_bwgs_solver_get_fmt
    procedure, nopass   :: get_id   => d_bwgs_solver_get_id
    procedure, pass(sv) :: descr    => d_bwgs_solver_descr
  end type amg_d_bwgs_solver_type

  private :: d_gs_solver_bld, d_gs_solver_apply, &
       &  d_gs_solver_free,   &
       &  d_gs_solver_descr,  d_gs_solver_sizeof, &
       &  d_gs_solver_default, d_gs_solver_dmp, &
       &  d_gs_solver_apply_vect, d_gs_solver_get_nzeros, &
       &  d_gs_solver_get_fmt, d_gs_solver_check,&
       &  d_gs_solver_is_iterative, &
       &  d_bwgs_solver_get_fmt, d_bwgs_solver_descr, &
       &  d_gs_solver_get_id, d_bwgs_solver_get_id, d_gs_solver_get_wrksize

  interface
    subroutine amg_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         & trans,work,wv,info,init,initu)
      import :: psb_desc_type, amg_d_gs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_
      implicit none
      type(psb_desc_type), intent(in)             :: desc_data
      class(amg_d_gs_solver_type), intent(inout) :: sv
      type(psb_d_vect_type),intent(inout)         :: x
      type(psb_d_vect_type),intent(inout)         :: y
      real(psb_dpk_),intent(in)                    :: alpha,beta
      character(len=1),intent(in)                   :: trans
      real(psb_dpk_),target, intent(inout)         :: work(:)
      type(psb_d_vect_type),intent(inout)         :: wv(:)
      integer(psb_ipk_), intent(out)                :: info
      character, intent(in), optional                :: init
      type(psb_d_vect_type),intent(inout), optional   :: initu
    end subroutine amg_d_gs_solver_apply_vect
    subroutine amg_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         & trans,work,wv,info,init,initu)
      import :: psb_desc_type, amg_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_
      implicit none
      type(psb_desc_type), intent(in)             :: desc_data
      class(amg_d_bwgs_solver_type), intent(inout) :: sv
      type(psb_d_vect_type),intent(inout)         :: x
      type(psb_d_vect_type),intent(inout)         :: y
      real(psb_dpk_),intent(in)                    :: alpha,beta
      character(len=1),intent(in)                   :: trans
      real(psb_dpk_),target, intent(inout)         :: work(:)
      type(psb_d_vect_type),intent(inout)         :: wv(:)
      integer(psb_ipk_), intent(out)                :: info
      character, intent(in), optional                :: init
      type(psb_d_vect_type),intent(inout), optional   :: initu
    end subroutine amg_d_bwgs_solver_apply_vect
  end interface

  interface
    subroutine amg_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info,init,initu)
      import :: psb_desc_type, amg_d_gs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_
      implicit none
      type(psb_desc_type), intent(in)      :: desc_data
      class(amg_d_gs_solver_type), intent(inout) :: sv
      real(psb_dpk_),intent(inout)         :: x(:)
      real(psb_dpk_),intent(inout)         :: y(:)
      real(psb_dpk_),intent(in)            :: alpha,beta
      character(len=1),intent(in)           :: trans
      real(psb_dpk_),target, intent(inout) :: work(:)
      integer(psb_ipk_), intent(out)        :: info
      character, intent(in), optional       :: init
      real(psb_dpk_),intent(inout), optional :: initu(:)
    end subroutine amg_d_gs_solver_apply
    subroutine amg_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
         & trans,work,info,init,initu)
      import :: psb_desc_type, amg_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_
      implicit none
      type(psb_desc_type), intent(in)      :: desc_data
      class(amg_d_bwgs_solver_type), intent(inout) :: sv
      real(psb_dpk_),intent(inout)         :: x(:)
      real(psb_dpk_),intent(inout)         :: y(:)
      real(psb_dpk_),intent(in)            :: alpha,beta
      character(len=1),intent(in)           :: trans
      real(psb_dpk_),target, intent(inout) :: work(:)
      integer(psb_ipk_), intent(out)        :: info
      character, intent(in), optional       :: init
      real(psb_dpk_),intent(inout), optional :: initu(:)
    end subroutine amg_d_bwgs_solver_apply
  end interface

  interface
    subroutine amg_d_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
      import :: psb_desc_type, amg_d_gs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
      type(psb_dspmat_type), intent(in), target, optional :: b
      class(psb_d_base_sparse_mat), intent(in), optional  :: amold
      class(psb_d_base_vect_type), intent(in), optional   :: vmold
      class(psb_i_base_vect_type), intent(in), optional   :: imold
    end subroutine amg_d_gs_solver_bld
    subroutine amg_d_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
      import :: psb_desc_type, amg_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_bwgs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
      type(psb_dspmat_type), intent(in), target, optional :: b
      class(psb_d_base_sparse_mat), intent(in), optional  :: amold
      class(psb_d_base_vect_type), intent(in), optional   :: vmold
      class(psb_i_base_vect_type), intent(in), optional   :: imold
    end subroutine amg_d_bwgs_solver_bld
  end interface

  interface
    subroutine amg_d_gs_solver_cnv(sv,info,amold,vmold,imold)
      import :: amg_d_gs_solver_type, psb_dpk_, &
           & psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
      class(amg_d_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
      class(psb_d_base_sparse_mat), intent(in), optional  :: amold
      class(psb_d_base_vect_type), intent(in), optional   :: vmold
      class(psb_i_base_vect_type), intent(in), optional   :: imold
    end subroutine amg_d_gs_solver_cnv
  end interface

  interface
    subroutine amg_d_gs_solver_dmp(sv,desc,level,info,prefix,head,solver,global_num)
      import :: psb_desc_type, amg_d_gs_solver_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, &
           & psb_ipk_
      implicit none
      class(amg_d_gs_solver_type), intent(in) :: sv
      type(psb_desc_type), intent(in)             :: desc
      integer(psb_ipk_), intent(in)              :: level
      integer(psb_ipk_), intent(out)             :: info
      character(len=*), intent(in), optional     :: prefix, head
      logical, optional, intent(in)              :: solver, global_num
    end subroutine amg_d_gs_solver_dmp
  end interface

  interface
    subroutine amg_d_gs_solver_clone(sv,svout,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_solver_type, amg_d_gs_solver_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_gs_solver_type), intent(inout)               :: sv
      class(amg_d_base_solver_type), allocatable, intent(inout) :: svout
      integer(psb_ipk_), intent(out)              :: info
    end subroutine amg_d_gs_solver_clone
  end interface

  interface
    subroutine amg_d_gs_solver_clone_settings(sv,svout,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_solver_type, amg_d_gs_solver_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_gs_solver_type), intent(inout) :: sv
      class(amg_d_base_solver_type), intent(inout) :: svout
      integer(psb_ipk_), intent(out)               :: info
    end subroutine amg_d_gs_solver_clone_settings
  end interface

  interface
    subroutine amg_d_gs_solver_clear_data(sv,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_gs_solver_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_gs_solver_type), intent(inout) :: sv
      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_d_gs_solver_clear_data
  end interface

contains

  subroutine d_gs_solver_default(sv)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv

    sv%sweeps = ione
    sv%eps    = dzero

    return
  end subroutine d_gs_solver_default

  subroutine d_gs_solver_check(sv,info)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='d_gs_solver_check'

    call psb_erractionsave(err_act)
    info = psb_success_

    call amg_check_def(sv%sweeps,&
         & 'GS Sweeps',ione,is_int_positive)

    if (info /= psb_success_) goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return

  end subroutine d_gs_solver_check

  subroutine d_gs_solver_cseti(sv,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv
    character(len=*), intent(in)                  :: what
    integer(psb_ipk_), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    integer(psb_ipk_)  :: err_act
    character(len=20)  :: name='d_gs_solver_cseti'

    info = psb_success_
    call psb_erractionsave(err_act)

    select case(psb_toupper(what))
    case('SOLVER_SWEEPS')
      sv%sweeps = val
    case default
      call sv%amg_d_base_solver_type%set(what,val,info,idx=idx)
    end select

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_gs_solver_cseti

  subroutine d_gs_solver_csetc(sv,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv
    character(len=*), intent(in)                  :: what
    character(len=*), intent(in)                  :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    integer(psb_ipk_)  :: err_act, ival
    character(len=20)  :: name='d_gs_solver_csetc'

    info = psb_success_
    call psb_erractionsave(err_act)

    call sv%amg_d_base_solver_type%set(what,val,info,idx=idx)

    if (info /= psb_success_) then
      info = psb_err_from_subroutine_
      call psb_errpush(info, name)
      goto 9999
    end if

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_gs_solver_csetc

  subroutine d_gs_solver_csetr(sv,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv
    character(len=*), intent(in)                  :: what
    real(psb_dpk_), intent(in)                     :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    integer(psb_ipk_)  :: err_act
    character(len=20)  :: name='d_gs_solver_csetr'

    call psb_erractionsave(err_act)
    info = psb_success_

    select case(psb_toupper(what))
    case('SOLVER_EPS')
      sv%eps = val
    case default
      call sv%amg_d_base_solver_type%set(what,val,info,idx=idx)
    end select

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_gs_solver_csetr

  subroutine d_gs_solver_free(sv,info)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(inout) :: sv
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_)  :: err_act
    character(len=20)  :: name='d_gs_solver_free'

    call psb_erractionsave(err_act)
    info = psb_success_

    call sv%l%free()
    call sv%u%free()

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_gs_solver_free

  subroutine d_gs_solver_descr(sv,info,iout,coarse)

    Implicit None

    ! Arguments
    class(amg_d_gs_solver_type), intent(in) :: sv
    integer(psb_ipk_), intent(out)             :: info
    integer(psb_ipk_), intent(in), optional    :: iout
    logical, intent(in), optional       :: coarse

    ! Local variables
    integer(psb_ipk_)      :: err_act
    character(len=20), parameter :: name='amg_d_gs_solver_descr'
    integer(psb_ipk_) :: iout_

    call psb_erractionsave(err_act)
    info = psb_success_
    if (present(iout)) then
      iout_ = iout
    else
      iout_ = psb_out_unit
    endif

    if (sv%eps<=dzero) then
      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
           &  sv%sweeps,' sweeps'
    else
      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
           &  sv%eps,' and maxit', sv%sweeps
    end if

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_gs_solver_descr

  function d_gs_solver_get_nzeros(sv) result(val)

    implicit none
    ! Arguments
    class(amg_d_gs_solver_type), intent(in) :: sv
    integer(psb_epk_) :: val
    integer(psb_ipk_)        :: i

    val = 0
    val = val + sv%l%get_nzeros()
    val = val + sv%u%get_nzeros()

    return
  end function d_gs_solver_get_nzeros

  function d_gs_solver_sizeof(sv) result(val)

    implicit none
    ! Arguments
    class(amg_d_gs_solver_type), intent(in) :: sv
    integer(psb_epk_) :: val
    integer(psb_ipk_)        :: i

    val = psb_sizeof_ip
    val = val + sv%l%sizeof()
    val = val + sv%u%sizeof()

    return
  end function d_gs_solver_sizeof

  function d_gs_solver_get_fmt() result(val)
    implicit none
    character(len=32)  :: val

    val = "Forward Gauss-Seidel solver"
  end function d_gs_solver_get_fmt

  function d_gs_solver_get_id() result(val)
    implicit none
    integer(psb_ipk_)  :: val

    val = amg_gs_
  end function d_gs_solver_get_id

  !
  ! If this is true, then the solver needs a starting
  ! guess. Currently only handled in JAC smoother.
  !
  function d_gs_solver_is_iterative() result(val)
    implicit none
    logical  :: val

    val = .true.
  end function d_gs_solver_is_iterative

  subroutine d_bwgs_solver_descr(sv,info,iout,coarse)

    Implicit None

    ! Arguments
    class(amg_d_bwgs_solver_type), intent(in) :: sv
    integer(psb_ipk_), intent(out)             :: info
    integer(psb_ipk_), intent(in), optional    :: iout
    logical, intent(in), optional       :: coarse

    ! Local variables
    integer(psb_ipk_)      :: err_act
    character(len=20), parameter :: name='amg_d_bwgs_solver_descr'
    integer(psb_ipk_) :: iout_

    call psb_erractionsave(err_act)
    info = psb_success_
    if (present(iout)) then
      iout_ = iout
    else
      iout_ = psb_out_unit
    endif

    if (sv%eps<=dzero) then
      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
           &  sv%sweeps,' sweeps'
    else
      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
           &  sv%eps,' and maxit', sv%sweeps
    end if

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)
    return
  end subroutine d_bwgs_solver_descr

  function d_bwgs_solver_get_fmt() result(val)
    implicit none
    character(len=32)  :: val

    val = "Backward Gauss-Seidel solver"
  end function d_bwgs_solver_get_fmt

  function d_bwgs_solver_get_id() result(val)
    implicit none
    integer(psb_ipk_)  :: val

    val = amg_bwgs_
  end function d_bwgs_solver_get_id

  function d_gs_solver_get_wrksize() result(val)
    implicit none
    integer(psb_ipk_)  :: val

    val = 2
  end function d_gs_solver_get_wrksize

end module amg_d_gs_solver


text:
C/MEMBR ADD NAME=WAXPY,SSI=0
c     Copyright INRIA
      subroutine waxpy(n,sr,si,xr,xi,incx,yr,yi,incy)
c!but
c
c     cette subroutine multiplie la constante complexe double
c     precision s (dont la partie reelle est sr et la partie
c     imaginaire si) par le vecteur complexe double precision
c     (dont les parties reelles de ses composantes sont dans
c     le vecteur double precision xr). le produit ainsi
c     obtenu est additionne au vecteur complexe y (dont les
c     parties reelles de ses composantes sont rangees dans le
c     vecteur double precision yr et les parties imaginaires
c     dans le vecteur double precision yr). le resultat de
c     l'addition reste dans y.
c
c!liste d'appel
c
c      subroutine waxpy(n,sr,si,xr,xi,incx,yr,yi,incy)
c
c     n: entier, taille des vecteurs traites
c
c     sr, si: double precision, parties reel et imaginaire de s
c
c     xr, xi: vecteurs double precision, parties rellees et
c     imaginaires, respectivement du vecteur complexe x.
c
c     yr, yi: vecteurs double precision, parties rellees et
c     imaginaires, respectivement du vecteur complexe y.
c
c     incx, incy: entiers, increments entre deux composantes
c     successives des vecteurs x et y.
c
c!auteur
c
c     cleve moler.- mathlab.
c
c!
      double precision sr,si,xr(*),xi(*),yr(*),yi(*)
      if (n .le. 0) return
      if (sr .eq. 0.0d+0 .and. si .eq. 0.0d+0) return
      ix = 1
      iy = 1
      if (incx.lt.0) ix = (-n+1)*incx + 1
      if (incy.lt.0) iy = (-n+1)*incy + 1
      do 10 i = 1, n
         yr(iy) = yr(iy) + sr*xr(ix) - si*xi(ix)
         yi(iy) = yi(iy) + sr*xi(ix) + si*xr(ix)
         ix = ix + incx
         iy = iy + incy
   10 continue
      return
      end

ext: f

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


! { dg-do compile }
! { dg-options "-Werror -O3" }
!
! Tests the fix for PR78619, in which the recursive use of 'f' at line 13
! caused an ICE.
!
! Contributed by Gerhard Steinmetz  <gerhard.steinmetz.fortran@t-online.de>
!
  print *, g(1.0) ! 'g' is OK
contains
  function f(x) result(z)
    real :: x, z
    z = sign(1.0, f) ! { dg-error "calling itself recursively|must be the same type" }
  end
  real function g(x)
    real :: x
    g = -1
    g = -sign(1.0, g) ! This is OK.
  end
end
! { dg-message "all warnings being treated as errors" "" { target *-*-* } 0 }


text:
program test
    use arithmetic_mod  ! importing arithmetic module
    implicit none

    ! declaring and initializing variables
    integer :: i_num_1 = 2, i_num_2 = 3
    real :: r_num_1 = 25.0, r_num_2 = 35.0
    complex :: com_num_1 = (7, 8), com_num_2 = (5, -7)

    ! Testing adder interface
    print "(a21)", "Test adder interface:"
    print "(a8,i1)", "2 + 3 = ", adder(i_num_1, i_num_2)
    print "(a14,f4.1)", "25.0 + 35.0 = ", adder(r_num_1, r_num_2)
    print "(a19,(f6.3,f6.3),(f6.3, f6.3))", "(7, 8) + (5, -7) = ", &
    adder(com_num_1, com_num_2)

    ! Testing subtractor interface
    print "(/,a26)", "Test subtractor interface:"
    print "(a8,i2)", "2 - 3 = ", subtractor(i_num_1, i_num_2)
    print "(a14,f5.1)", "25.0 - 35.0 = ", subtractor(r_num_1, r_num_2)
    print "(a19,(f5.3,f7.3))", "(7, 8) - (5, -7) = ", &
    subtractor(com_num_1, com_num_2)

    ! Testing multiplier interface
    print "(/,a26)", "Test multiplier interface:"
    print "(a8,i2)", "2 * 3 = ", multiplier(i_num_1, i_num_2)
    print "(a14,f5.1)", "25.0 * 35.0 = ", multiplier(r_num_1, r_num_2)
    print "(a19,(f6.3,f7.3))", "(7, 8) * (5, -7) = ", &
    multiplier(com_num_1, com_num_2)

    ! Testing divider interface
    print "(/,a23)", "Test divider interface:"
    print "(a8,i2)", "2 / 3 = ", divider(i_num_1, i_num_2)
    print "(a14,f8.6)", "25.0 / 35.0 = ", divider(r_num_1, r_num_2)
    print "(a19,(f9.6,f9.6))", "(7, 8) / (5, -7) = ", &
    divider(com_num_1, com_num_2)
    ! devided by zero tests
    print "(/,a22)", "Divided by zero tests:"
    print "(a8,i1)", "2 / 0 = ", divider(i_num_1, 0)
    print "(a13,f8.6)", "25.0 / 0.0 = ", divider(r_num_1, 0.0)
    print "(a18,(f9.6,f9.6))", "(7, 8) / (0, 0) = ", &
    divider(com_num_1, (0.0, 0.0))

end program test

max_stars_repo_name: EverLookNeverSee/FCS

max_issues_repo_name: EverLookNeverSee/FCS

max_forks_repo_name: EverLookNeverSee/FCS


text:
      program wrtdat
!
!--------------------------------------------------------------------
!  Tim Marchok, ????
!     This program reads in data from an ascii text file and writes
!     the data back out in GrADS format (direct access).
!     The input for this program is an ascii file containing
!     mean track/intensity .
!     The data in the ascii file should be in f6.1 format to be read in.
!--------------------------------------------------------------------
!  Fanglin Yang, January 2016
!    Include standard deviation and extend to 168 hours.
!--------------------------------------------------------------------
!
      parameter(npts=15)
      real cdat(npts,100),cdat1(npts),cdat2(npts)
      character  INFILE*65,OUTFILE*65,cmodname*8
!
      namelist/namin/INFILE

!     Read name of input file from the namelist....
      read (5,namin,end=1000)
 1000 continue

      print *,'After namelist read, infile= ',infile

      open (unit=11,file=INFILE
     &  ,access='sequential',action='read')

      OUTFILE = INFILE(1:LASTDOT(INFILE))//'.gr'
      open (unit=51,file=OUTFILE,form='unformatted',status='unknown')

      print *,'Name of input ascii track file is  ',infile
      print *,'Name of output GrADS track file is ',outfile
      print *,' '

      n = 0
      do while (.true.)
        n = n + 1
        read (11,31,end=99,err=100) cmodname,(cdat(i,n),i=1,npts)
        write (6,31) cmodname,(cdat(i,n),i=1,npts)
      enddo
  31  format (1x,a4,3x,15(1x,f6.1))
  99  continue

      mm=n/2
      print *,'Number of models = ',mm
      do k=1,mm
        do i=1,npts
         cdat1(i)=cdat(i,k)
         cdat2(i)=cdat(i,mm+k)
        enddo
        write (51) cdat1   !mean track error
        write (51) cdat2   !track error standard deviation
      enddo
      goto 200
 100  print *,'ERROR writing output file'
 200  continue
      stop
      end
c
c--------------------------------------------
c
c--------------------------------------------
      FUNCTION LASTDOT (STRING)
C
C**   RETURNS THE POSITION OF THE LAST CHARACTER of a
c     OF A STRING before the last dot in the string.  For
c     example, in the string trkdat.ep.ascii, it would
c     return "9", for the position of the "p".
C
      CHARACTER*(*) STRING
C
      LAST = LEN(STRING)
C
      DO 10 I = LAST,1,-1
        IF (STRING(I:I).EQ.'.') then
          itmp = i - 1
          GO TO 20
        endif
   10 CONTINUE
C
      LASTDOT = 0
      RETURN
C
   20 LASTDOT = itmp
      RETURN
C
      END

max_stars_repo_name: yangfanglin/hurtrack

max_issues_repo_name: yangfanglin/hurtrack

max_forks_repo_name: yangfanglin/hurtrack


text:
      PROGRAM ZAMTST
C----------------------------------------------------------------------C
C                                                                      C
C     ZAMBRETTI FORECASTER TEST HARNESS                                C
C                                                                      C
C     AUTHOR: TJH 30-08-2019                                           C
C                                                                      C
C----------------------------------------------------------------------C
      INTEGER IMONTH, ILOOP
      REAL PRESS, PDIFF
      CHARACTER*2 CIND
      CHARACTER*20 CFORE
C
C     MONTH = APRIL (SPRING / AUTUMN CONDITIONS)
C
      IMONTH = 4
C
C     FALLING PRESSURE
C
      PDIFF = -2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"FALLING PRESSURE, SPRING"
      WRITE(*,*)"------------------------"
      DO 100 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  100 CONTINUE
C
C     STEADY PRESSURE
C
      PDIFF = 0.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"STEADY PRESSURE, SPRING"
      WRITE(*,*)"------------------------"
      DO 200 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  200 CONTINUE
C
C     RISING PRESSURE
C
      PDIFF = 2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"RISING PRESSURE, SPRING"
      WRITE(*,*)"------------------------"
      DO 300 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  300 CONTINUE
C
C     MONTH = JANUARY (WINTER CONDITIONS)
C
      IMONTH = 1
C
C     FALLING PRESSURE
C
      PDIFF = -2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"FALLING PRESSURE, WINTER"
      WRITE(*,*)"------------------------"
      DO 400 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  400 CONTINUE
C
C     STEADY PRESSURE
C
      PDIFF = 0.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"STEADY PRESSURE, WINTER"
      WRITE(*,*)"------------------------"
      DO 500 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  500 CONTINUE
C
C     RISING PRESSURE
C
      PDIFF = 2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"RISING PRESSURE, WINTER"
      WRITE(*,*)"------------------------"
      DO 600 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  600 CONTINUE
C
C     MONTH = JULY (SUMMER CONDITIONS)
C
      IMONTH = 7
C
C     FALLING PRESSURE
C
      PDIFF = -2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"FALLING PRESSURE, SUMMER"
      WRITE(*,*)"------------------------"
      DO 700 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  700 CONTINUE
C
C     STEADY PRESSURE
C
      PDIFF = 0.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"STEADY PRESSURE, SUMMER"
      WRITE(*,*)"------------------------"
      DO 800 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  800 CONTINUE
C
C     RISING PRESSURE
C
      PDIFF = 2.0
      WRITE(*,*)"------------------------"
      WRITE(*,*)"RISING PRESSURE, SUMMER"
      WRITE(*,*)"------------------------"
      DO 900 ILOOP=940,1060
         PRESS=REAL(ILOOP)
         CALL ZFCAST(IMONTH, PRESS, PDIFF, CIND, CFORE)
         WRITE(*,*) PRESS, " ", CIND, " ", CFORE
  900 CONTINUE
      STOP
      END

max_stars_repo_name: psychotimmy/weather

max_issues_repo_name: psychotimmy/weather

max_forks_repo_name: psychotimmy/weather


text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK3.0                       .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c
c
c ... file divgs.f
c
c     this file includes documentation and code for
c     subroutine divgs          i
c
c ... files which must be loaded with divgs.f
c
c     sphcom.f, hrfft.f, vhags.f, shsgs.f, gaqd.f
c
c
c     subroutine divgs(nlat,nlon,isym,nt,divg,idiv,jdiv,br,bi,mdb,ndb,
c    +                 wshsgs,lshsgs,work,lwork,ierror)
c
c     given the vector spherical harmonic coefficients br and bi, precomputed
c     by subroutine vhags for a vector field (v,w), subroutine divgs
c     computes the divergence of the vector field in the scalar array divg.
c     divg(i,j) is the divergence at the gaussian colatitude point theta(i)
c     (see nlat as input parameter) and east longitude
c
c            lambda(j) = (j-1)*2*pi/nlon
c
c     on the sphere.  i.e.
c
c            dv(i,j) = 1/sint*[ d(sint*v(i,j))/dtheta + d(w(i,j))/dlambda ]
c
c     where sint = sin(theta(i)).  w is the east longitudinal and v
c     is the colatitudinal component of the vector field from which
c     br,bi were precomputed
c
c
c     input parameters
c
c     nlat   the number of points in the gaussian colatitude grid on the
c            full sphere. these lie in the interval (0,pi) and are computed
c            in radians in theta(1) <...< theta(nlat) by subroutine gaqd.
c            if nlat is odd the equator will be included as the grid point
c            theta((nlat+1)/2).  if nlat is even the equator will be
c            excluded as a grid point and will lie half way between
c            theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3.
c            note: on the half sphere, the number of grid points in the
c            colatitudinal direction is nlat/2 if nlat is even or
c            (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than zero. the axisymmetric case corresponds to nlon=1.
c            the efficiency of the computation is improved when nlon
c            is a product of small prime numbers.
c
c
c     isym   a parameter which determines whether the divergence is
c            computed on the full or half sphere as follows:
c
c      = 0
c
c            the symmetries/antsymmetries described in isym=1,2 below
c            do not exist in (v,w) about the equator.  in this case the
c            divergence is neither symmetric nor antisymmetric about
c            the equator.  the divergence is computed on the entire
c            sphere.  i.e., in the array divg(i,j) for i=1,...,nlat and
c            j=1,...,nlon.
c
c      = 1
c
c            w is antisymmetric and v is symmetric about the equator.
c            in this case the divergence is antisymmetyric about
c            the equator and is computed for the northern hemisphere
c            only.  i.e., if nlat is odd the divergence is computed
c            in the array divg(i,j) for i=1,...,(nlat+1)/2 and for
c            j=1,...,nlon.  if nlat is even the divergence is computed
c            in the array divg(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c      = 2
c            w is symmetric and v is antisymmetric about the equator
c            in this case the divergence is symmetyric about the
c            equator and is computed for the northern hemisphere
c            only.  i.e., if nlat is odd the divergence is computed
c            in the array divg(i,j) for i=1,...,(nlat+1)/2 and for
c            j=1,...,nlon.  if nlat is even the divergence is computed
c            in the array divg(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c
c     nt     nt is the number of scalar and vector fields.  some
c            computational efficiency is obtained for multiple fields.
c            in the program that calls divgs, the arrays br,bi, and divg
c            can be three dimensional corresponding to an indexed multiple
c            vector field.  in this case multiple scalar synthesis will
c            be performed to compute the divergence for each field.  the
c            third index is the synthesis index which assumes the values
c            k=1,...,nt.  for a single synthesis set nt = 1.  the
c            description of the remaining parameters is simplified by
c            assuming that nt=1 or that all the arrays are two dimensional.
c
c     idiv   the first dimension of the array divg as it appears in
c            the program that calls divgs. if isym = 0 then idiv
c            must be at least nlat.  if isym = 1 or 2 and nlat is
c            even then idiv must be at least nlat/2. if isym = 1 or 2
c            and nlat is odd then idiv must be at least (nlat+1)/2.
c
c     jdiv   the second dimension of the array divg as it appears in
c            the program that calls divgs. jdiv must be at least nlon.
c
c     br,bi  two or three dimensional arrays (see input parameter nt)
c            that contain vector spherical harmonic coefficients
c            of the vector field (v,w) as computed by subroutine vhags.
c     ***    br and bi must be computed by vhags prior to calling
c            divgs.
c
c     mdb    the first dimension of the arrays br and bi as it
c            appears in the program that calls divgs. mdb must be at
c            least min0(nlat,nlon/2) if nlon is even or at least
c            min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     ndb    the second dimension of the arrays br and bi as it
c            appears in the program that calls divgs. ndb must be at
c            least nlat.
c
c
c     wshsgs an array which must be intialized by subroutine shsgsi.
c            once initialized,
c            wshsgs can be used repeatedly by divgs as long as nlon
c            and nlat remain unchanged.  wshsgs must not be altered
c            between calls of divgs.
c
c
c     lshsgs the dimension of the array wshsgs as it appears in the
c            program that calls divgs. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsgs must be at least
c
c               nlat*(3*(l1+l2)-2)+(l1-1)*(l2*(2*nlat-l1)-3*l1)/2+nlon+15
c
c
c     work   a work array that does not have to be saved.
c
c     lwork  the dimension of the array work as it appears in the
c            program that calls divgs. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2                    if nlat is even or
c               l2 = (nlat+1)/2                if nlat is odd
c
c            if isym = 0 then lwork must be at least
c
c               nlat*((nt+1)*nlon+2*nt*l1+1)
c
c            if isym > 0 then lwork must be at least
c
c               (nt+1)*l2*nlon+nlat*(2*nt*l1+1)
c
c     **************************************************************
c
c     output parameters
c
c
c    divg   a two or three dimensional array (see input parameter nt)
c           that contains the divergence of the vector field (v,w)
c           whose coefficients br,bi where computed by subroutine
c           vhags.  divg(i,j) is the divergence at the gaussian colatitude
c           point theta(i) and longitude point lambda(j) = (j-1)*2*pi/nlon.
c           the index ranges are defined above at the input parameter
c           isym.
c
c
c    ierror = 0  no errors
c           = 1  error in the specification of nlat
c           = 2  error in the specification of nlon
c           = 3  error in the specification of isym
c           = 4  error in the specification of nt
c           = 5  error in the specification of idiv
c           = 6  error in the specification of jdiv
c           = 7  error in the specification of mdb
c           = 8  error in the specification of ndb
c           = 9  error in the specification of lshsgs
c           = 10 error in the specification of lwork
c **********************************************************************
c
c
      subroutine divgs(nlat,nlon,isym,nt,divg,idiv,jdiv,br,bi,mdb,ndb,
     +                 wshsgs,lshsgs,work,lwork,ierror)

      dimension divg(idiv,jdiv,nt),br(mdb,ndb,nt),bi(mdb,ndb,nt)
      dimension wshsgs(lshsgs),work(lwork)
c
c     check input parameters
c
      ierror = 1
      if(nlat .lt. 3) return
      ierror = 2
      if(nlon .lt. 4) return
      ierror = 3
      if (isym.lt.0 .or. isym.gt.2) return
      ierror = 4
      if(nt .lt. 0) return
      ierror = 5
      imid = (nlat+1)/2
      if((isym.eq.0 .and. idiv.lt.nlat) .or.
     1   (isym.gt.0 .and. idiv.lt.imid)) return
      ierror = 6
      if(jdiv .lt. nlon) return
      ierror = 7
      if(mdb .lt. min0(nlat,(nlon+1)/2)) return
      mmax = min0(nlat,(nlon+2)/2)
      ierror = 8
      if(ndb .lt. nlat) return
      ierror = 9
      imid = (nlat+1)/2
      lpimn = (imid*mmax*(nlat+nlat-mmax+1))/2
c     check permanent work space length
      l2 = (nlat+1)/2
      l1 = min0((nlon+2)/2,nlat)
      lp=nlat*(3*(l1+l2)-2)+(l1-1)*(l2*(2*nlat-l1)-3*l1)/2+nlon+15
      if(lshsgs.lt.lp) return
      ierror = 10
c
c     verify unsaved work space (add to what shses requires, file f3)
c
      ls = nlat
      if(isym .gt. 0) ls = imid
      nln = nt*ls*nlon
c
c     set first dimension for a,b (as requried by shses)
c
      mab = min0(nlat,nlon/2+1)
      mn = mab*nlat*nt
      if(lwork.lt. nln+ls*nlon+2*mn+nlat) return
      ierror = 0
c
c     set work space pointers
c
      ia = 1
      ib = ia+mn
      is = ib+mn
      iwk = is+nlat
      lwk = lwork-2*mn-nlat
      call divgs1(nlat,nlon,isym,nt,divg,idiv,jdiv,br,bi,mdb,ndb,
     +work(ia),work(ib),mab,work(is),wshsgs,lshsgs,work(iwk),lwk,
     +ierror)
      return
      end

      subroutine divgs1(nlat,nlon,isym,nt,divg,idiv,jdiv,br,bi,mdb,ndb,
     +                  a,b,mab,sqnn,wshsgs,lshsgs,wk,lwk,ierror)
      dimension divg(idiv,jdiv,nt),br(mdb,ndb,nt),bi(mdb,ndb,nt)
      dimension a(mab,nlat,nt),b(mab,nlat,nt),sqnn(nlat)
      dimension wshsgs(lshsgs),wk(lwk)
c
c     set coefficient multiplyers
c
      do 1 n=2,nlat
      fn = float(n-1)
      sqnn(n) = sqrt(fn*(fn+1.))
    1 continue
c
c     compute divergence scalar coefficients for each vector field
c
      do 2 k=1,nt
      do 3 n=1,nlat
      do 4 m=1,mab
      a(m,n,k) = 0.0
      b(m,n,k) = 0.0
    4 continue
    3 continue
c
c     compute m=0 coefficients
c
      do 5 n=2,nlat
      a(1,n,k) = -sqnn(n)*br(1,n,k)
      b(1,n,k) = -sqnn(n)*bi(1,n,k)
    5 continue
c
c     compute m>0 coefficients using vector spherepack value for mmax
c
      mmax = min0(nlat,(nlon+1)/2)
      do 6 m=2,mmax
      do 7 n=m,nlat
      a(m,n,k) = -sqnn(n)*br(m,n,k)
      b(m,n,k) = -sqnn(n)*bi(m,n,k)
    7 continue
    6 continue
    2 continue
c
c     synthesize a,b into divg
c
      call shsgs(nlat,nlon,isym,nt,divg,idiv,jdiv,a,b,
     +           mab,nlat,wshsgs,lshsgs,wk,lwk,ierror)
      return
      end

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
    SUBROUTINE testWrite(ifltab1, messageUnit, status)
!
!    Test basic writing and reading functions
!    Internal function
!
      implicit none
!
      integer messageUnit, status
      integer(8) ifltab1(600)
!
      real data1(1000), data2(1000)
      integer userHeader(1000)
      character cpath1*100, cpath2*100
      character cunits*12, ctype*12
      character ctemp*10
      integer istat, ierr
      integer i, npath1, nuserHeader, nvals
      logical lfound

      common /lchk/ lcheck
      logical lcheck
!
      integer iplan, jplan
!     IPLAN = 0   Always write
!     IPLAN = 1   Only write if new record
!     IPLAN = 2   Only write if old record
!
!     JPLAN = 0:  Use limits NHEAD and NDATA, but do not update NHEAD,
!                 and NDATA
!     JPLAN = 1:  Use limits, but update NHEAD and NDATA
!     JPLAN = 2:  Same as 1, but do not write a warning message if data
!                 is not found
!
!
      WRITE(messageUnit, *)'Begin testWrite'
1     Format('Basic zwrite test FAILED')
!
!     This is a fundemental check on the write/check/read capability of DSS
      cpath1 = '/Test Write/1/c/d/e/f/'
      call chrlnb(cpath1, npath1)
      !  Simple zwrite check
      nvals = 1000
      do  i=1,nvals
        data1(i) = float(i)
      end do

    !call zset('mlvl', ' ', 17)

      iplan = 0
      call zwritea(ifltab1, cpath1, npath1, userHeader, 0, data1, nvals, iplan, lfound)
      call zinqir(ifltab1, 'error', ctemp, status)
      if (status.ne.0) go to 900
!
!     Now check that the record is there
      call zcheck(ifltab1, cpath1, npath1, nuserHeader, nvals, lfound)
      call zinqir(ifltab1, 'error', ctemp, status)
      if (status.ne.0) go to 900
      if (.not.lfound) Then
        write(messageUnit, 1)
        write(messageUnit, 100) cpath1
 100    format('First record written in DB NOT found: ', A)
        go to 910
     endif
     if (nuserHeader .ne. 0) Then
        write(messageUnit, 1)
        write(messageUnit, 101) nuserHeader
 101    format('Invalid length for User Header: ', i8)
        go to 910
     endif
     if (nvals .ne.1000) Then
        write(messageUnit, 1)
        write(messageUnit, 102) nvals
 102    format('Invalid length for Data Area: ', i8)
        go to 910
     endif
!
!
!   Now read the record and verify that the values are correct
      jplan = 1
      nvals = 1000
      call zreada(ifltab1, cpath1, npath1, userHeader, nuserHeader, data2, nvals, jplan, lfound)
      call zinqir(ifltab1, 'error', ctemp, status)
      if (status.ne.0) go to 900
      if (.not.lfound) Then
        write(messageUnit, 1)
        write(messageUnit, 200) cpath1
 200    format('First record written in DB NOT found: ', A)
        go to 910
     endif
     if (nuserHeader .ne. 0) Then
        write(messageUnit, 1)
        write(messageUnit, 201) nuserHeader
 201    format('Invalid length for User Header: ', i8)
        go to 910
     endif
     if (nvals .ne.1000) Then
        write(messageUnit, 1)
        write(messageUnit, 202) nvals
 202    format('Invalid length for Data Area: ', i8)
        go to 910
     endif
!
!     Check data values
      do i=1,nvals
        if (data2(i).ne.data1(i)) then
            write(messageUnit, 1)
            write(messageUnit, 245)
            write(messageUnit, 246) i, data1(i), data2(i)
 245         Format('Values read do not match those written')
 246         Format('At ordinate: ', I8, ' Values: ', 2F8.2)
            go to 900
        endif
      end do

!
!
    !call zset('mlvl', ' ', 17)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!  Tests where the hash codes are hardwired
    !!!  To use, set the testing code on in zhash and
    !! boolDoubleCheckPath in zcheck and zupdateBin
    !!
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 50, .true., status)
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 100, .true., status)
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 200, .true., status)
    !call writeARec(ifltab1, '/g/h/i/j/k/l/', 1000, .true., status)
!    do 30 i=1,800
!    write (cpath1, 29) i
29  format('/a/b/c/d/e/', i3.3, '/')
!    call chrlnb(cpath1, npath1)
!    call writeARec(ifltab1, cpath1(1:npath1), i, .false., status)
!    if (status.ne.0) then
!        write (*,*)'ISSSUES!'
!        return
!     endif
 30 continue
!    if (lcheck) call zcheckFile(ifltab1, status)
!    call zclose(ifltab1)
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 50, .true., status)
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 100, .true., status)
!    call writeARec(ifltab1, '/a/b/c/d/e/f/', 200, .true., status)
    !call writeARec(ifltab1, '/g/h/i/j/k/l/', 1000, .true., status)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
 800   Continue
       if (lcheck) call zcheckFile(ifltab1, status)
       if (status.eq.0) then
            Write(messageUnit, 810)
 810        Format('Basic zwrite Passed Successfully')
        else
            Write(messageUnit, 820)
 820        Format('File Check FAILED in zwrite')
        endif
       return
!
 900    Continue
        write(messageUnit, 1)
 910    Continue
        if (status.eq.0) status = -1
        return
        end

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


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.

      subroutine rwfunc ( iopt   , nitem  , nvals  , item   , nrec   ,
     &                    nhtot  , ifact  , dtflg  , dtflg3 , lununf ,
     &                    iwidth , ioutpt , ierr   )

!       Deltares Software Centre

!>\file
!>           Reads function values (Fourier and Harmonic components)
!>
!>           Harmonic components for iopt = 3 and Fouriers for iopt = 4/n
!>           Reads:
!>           - number of frequencies exclusive of average value
!>           - matrix(nvals*nitem) of average values
!>           - the base frequency (Fourrier series only)
!>           Per frequency:
!>           - period (only for Harmonics)
!>           - a phase of the component
!>           - matrix(nvals*nitem) of values for the component
!>           Notes:
!>           - because there is an average value, there are nharm+1 matrices
!>           - because there is a phase, there are 1+nvals*nitem values per component
!>           - the periods are filled in for the Fouriers, so they look further more the same
!>           - the zero'th period contains the value nvals*nitem
!>           - the zero'th phase contains the number of harmoncs
!>           The routine is entered with ierr = -1  for boundaries and -2 for wastes

!     Created            : March 1988  by M.E. Sileon

!     Modified           : April 1997  by R. Bruinsma : Tokenized input data file reading added
!                        : May   2011  by Leo Postma  : Fortran90 look and feel,
!                                                       no external memory needed

!     Subroutines called : conver
!                          cnvtim

!     Functions called   : gettok tokenized input data file reading

!     Logical units      : lunut  = unit formatted output file
!                          lununf = unit unformatted output file

      use rd_token       ! for the reading of tokens
      use timers       !   performance timers

      implicit none

!     Parameters

!     kind        function         name            Descriptipon

      integer(4), intent(in   ) :: iopt          !< 3 Harmonics, 4 Fourier
      integer(4), intent(in   ) :: nitem         !< number of input items
      integer(4), intent(in   ) :: nvals         !< number of values per item
      integer(4), intent(in   ) :: item (nitem)  !< item numbers
      integer(4), intent(inout) :: nrec          !< number of harmonic records
      integer(4), intent(inout) :: nhtot         !< total harmonic array space
      integer(4), intent(in   ) :: ifact         !< factor between clocks
      logical   , intent(in   ) :: dtflg         !< "date"-format
      logical   , intent(in   ) :: dtflg3        !< 'date'-format (F;ddmmhhss,T;yydddhh)
      integer(4), intent(in   ) :: lununf        !< unit nr unformatted file
      integer(4), intent(in   ) :: iwidth        !< width of theoutput file
      integer(4), intent(in   ) :: ioutpt        !< how extensive output ?
      integer(4), intent(inout) :: ierr          !< error count

      integer(4)                   ndim          ! total size of the matrix
      integer(4)                   nhar          ! number of harmonics
      integer(4)                   ibase         ! base period of Fouriers
      integer(4)                   ierr2         ! error hlp variable
      integer(4)                   i, k          ! loop variables
      integer(4)                   ib, ie        ! limits for printed output
      integer(4)                   i1, i2        ! print loop counters
      logical                      bound         ! true if boundary processing
      logical                      waste         ! true if waste processing
      integer(4), allocatable   :: iperio(:)     ! workspace for frequencies
      real   (4), allocatable   :: value (:,:)   ! workspace for values
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "rwfunc", ithndl )

      bound = .false.
      waste = .false.
      if ( ierr .eq. -1 ) bound = .true.
      if ( ierr .eq. -2 ) waste = .true.
      ierr  = 0
      ndim  = nitem*nvals

      if ( gettoken( nhar, ierr2 ) .gt. 0 ) goto 100
      allocate ( iperio(nhar+1), value(ndim+1,nhar+1) )

      select case ( iopt )

         case ( 3 )        !      read values if IOPT = 3 ( harmonic function )

            do k = 2, ndim+1
               if ( gettoken( value(k,1), ierr2 ) .gt. 0 ) goto 100
            enddo
            do i = 2, nhar+1
               if ( gettoken( iperio(i), ierr2 ) .gt. 0 ) goto 100
               do k = 1, ndim+1
                  if ( gettoken( value(k,i), ierr2 ) .gt. 0 ) goto 100
               enddo
            enddo
            call conver ( iperio(2) , nhar  , ifact  , dtflg , dtflg3 )
            value(1,1) = float(nhar)

         case ( 4 )        !      read values if IOPT = 4 ( fourier function )

            if ( gettoken( ibase, ierr2 ) .gt. 0 ) goto 100
            do k = 2, ndim+1
               if ( gettoken( value(k,1), ierr2 ) .gt. 0 ) goto 100
            enddo
            do i = 2, nhar+1
               do k = 1, ndim+1
                  if ( gettoken( value(k,i), ierr2 ) .gt. 0 ) goto 100
              enddo
            enddo
            call cnvtim ( ibase , ifact , dtflg , dtflg3 )
            value(1,1) = float(nhar)
            do i = 2,nhar+1
               iperio(i)  = ibase/(i-1)
            enddo

      end select

!        control writing

      if ( ioutpt .lt. 4 ) then
         write ( lunut , 2070 )
      else
         if ( iopt .eq. 3 ) then
            write ( lunut, 2000 ) nhar
         else
            write ( lunut, 2010 ) nhar, ibase
         endif
         write ( lunut, 2020 )
         do i1 = 1, nvals, iwidth
            write ( lunut, 2030 ) (k,k=i1,min(i1+iwidth-1,nvals))
            do i2 = 1, nitem
               ib = (i2-1)*nvals+1+i1
               ie = (i2-1)*nvals+1+min(i1+iwidth-1,nvals)
               write ( lunut, 2040 ) item(i2),(value(k,1),k=ib,ie)
            enddo
         enddo
         do i  = 2,nhar+1
            write ( lunut, 2050 ) iperio(i),value(1,i)
            if ( iperio(i) .le. 0 ) then
               write ( lunut, 2060 )
               ierr = ierr+1
            endif
            do i1 = 1, nvals, iwidth
               write ( lunut, 2030 ) (k,k=i1,min(i1+iwidth-1,nvals))
               do i2 = 1, nitem
                  ib = (i2-1)*nvals+1+i1
                  ie = (i2-1)*nvals+1+min(i1+iwidth-1,nvals)
                  write ( lunut, 2040 ) item(i2),(value(k,i),k=ib,ie)
               enddo
            enddo
         enddo
      endif

!        calculate new settings

      iperio(1)  = ndim
      if ( bound .or. waste ) write ( lununf ) nhar+1
      do i = 1, nhar+1
         write ( lununf ) iperio(i) , (value(k,i),k=1,ndim+1)
      enddo
      nrec  = nrec  +  nhar+1
      nhtot = nhtot + (ndim+1)*(nhar+1)
      if (timon) call timstop( ithndl )
      return

  100 ierr = ierr + 1
      if (timon) call timstop( ithndl )
      return

!        output formats

 2000 format( /,' Number of harmonics:',I4)
 2010 format( /,' Number of Fouriers:',I4,' base period:',I10)
 2020 format(   ' Mean values :')
 2030 format(   '      Item',I8,9I12 )
 2040 format(           I10,2X,1P,10E12.4)
 2050 format(   ' Period:',I10,'   Phase:',1P,E12.4)
 2060 format(   ' ERROR, PERIOD is less or equal to zero!' )
 2070 format(   ' Printed output for output option 4 and higher !' )

      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
c=======================================================================
c
c     subroutine OPVOL
c
c     Optimization utility for ALLOCVOL
c
c-----------------------------------------------------------------------
      SUBROUTINE opvol ( n, cov, rho, mu, neq, nin, ccst, bcst,
     &                  cinf, csup,
     &                  iwork, dwork, wopt, lagr, info )
c-----------------------------------------------------------------------
c
c     INPUT
c            n      : size of portfolio                          integer
c            cov    : covariance matrix (n*n)                     double
c            rho    : mean returns vector (n)                     double
c            mu     : performance target                          double
c            neq    : number of initial equality constraints     integer
c            nin    : number of initial inequality constraints   integer
c            ccst   : matrix of constraints (n*(neq + nin))       double
c            bcst   : vector of constraints (neq + nin)           double
c            cinf   : assets inferior limit (n)                   double
c            csup   : assets superior limit (n)                   double
c
c     WORKSPACE
c            iwork  : 3*n + 2*nin + neq + 3                      integer
c            dwork  : n*(n+neq+nin+8) + neq  + 3*nin + 3          double

c     OUTPUT
c            wopt   : optimal portfolio (n)                       double
c            lagr   : Lagrange (n + neq + 1 + nin)                double
c            info   : = 0 successful exit                        integer
c
c     CALL
c            CTMV : utility function (cf. utalloc.f)
c            QP   : quadratic solver (vectorized version)
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin
      DOUBLE PRECISION mu
      DOUBLE PRECISION cov(*), rho(*), cinf(*), csup(*), wopt(*),
     &                 ccst(n,*), bcst(*), lagr(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER neqtot, nintot, piw, pdlin, pdw, pdbcs, pdccs
c
c     external subroutines
      EXTERNAL ctmv, qp
c
c-----------------------------------------------------------------------
c
c     initializations
      info = 0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piw  = 1
c     piw  : pointer for QP internal workspace, who needs
c            ( 3*n + 2*nintot + neqtot + 1 )
c             with neqtot = neq
c                  nintot = nin + 1
c
c     Total size of iwork array = ( 3*n + 2*nin + neq + 3 )
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pdlin  = 1
c     pdlin  : pointer for linear part vector, so ( n ) more
      pdbcs  = pdlin +  ( n )
c     pdbcs  : pointer for the constraints vector ( nintot + neqtot ),
c              so ( nin + 1 + neq ) more
      pdccs  = pdbcs +  ( nin + 1 + neq )
c     pdccs  : pointer for the constraints matrix
c              n*( nintot + neqtot ),
c           so n*( nin + 1 + neq ) more
      pdw    = pdccs + ( n*( nin + 1 + neq ) )
c     pdw    : pointer for QP internal workspace who needs
c              (n*n + 6*n + 2*nintot)
c
c     Total size of dwork array = n
c                               + nin + 1 + neq
c                               + n*( nin + 1 + neq )
c                               + n*(n + 6) + 2*(nin + 1)
c                               = n*(n + neq + nin + 8) + neq  + 3*nin + 3
c
c      = n*(n + neq + nin + 8) + neq  + 3*nin + 3
c
c-----------------------------------------------------------------------
c
c     construction of the constraints matrix and vector, and linear part
      CALL ctmv ( n, rho, mu, neq, nin, ccst, bcst, neqtot, nintot,
     &            dwork(pdccs), dwork(pdbcs), dwork(pdlin))
c
c     quadratic solver
      CALL qp ( n, cov, dwork(pdlin), dwork(pdccs),
     &          dwork(pdbcs), cinf, csup, neqtot, nintot,
     &          iwork(piw), dwork(pdw), lagr, wopt,
     &          info )
c
      RETURN
      END


! { dg-do run }
! { dg-add-options ieee }
! PR48589 Invalid G0/G0.d editing for NaN/infinity
! Test case by Thomas Henlich
program test_g0_special

    call check_all("(g10.3)", "(f10.3)")
    call check_all("(g10.3e3)", "(f10.3)")
    call check_all("(spg10.3)", "(spf10.3)")
    call check_all("(spg10.3e3)", "(spf10.3)")
    !print *, "-----------------------------------"
    call check_all("(g0)", "(f0.0)")
    call check_all("(g0.15)", "(f0.0)")
    call check_all("(spg0)", "(spf0.0)")
    call check_all("(spg0.15)", "(spf0.0)")
contains
    subroutine check_all(fmt1, fmt2)
        character(len=*), intent(in) :: fmt1, fmt2
        real(8) :: one = 1.0D0, zero = 0.0D0, nan, pinf, minf

        nan = zero / zero
        pinf = one / zero
        minf = -one / zero
        call check_equal(fmt1, fmt2, nan)
        call check_equal(fmt1, fmt2, pinf)
        call check_equal(fmt1, fmt2, minf)
    end subroutine check_all
    subroutine check_equal(fmt1, fmt2, r)
        real(8), intent(in) :: r
        character(len=*), intent(in) :: fmt1, fmt2
        character(len=80) :: s1, s2

        write(s1, fmt1) r
        write(s2, fmt2) r
        if (s1 /= s2) STOP 1
        !if (s1 /= s2) print "(6a)", trim(fmt1), ": '", trim(s1), "' /= '", trim(s2), "'"
	!print "(6a)", trim(fmt1), ": '", trim(s1), "' /= '", trim(s2), "'"
    end subroutine check_equal
end program test_g0_special