text:
CCC   'viscoop.f' COMPUTES INTERNAL REGION STRESSES FOR VISCOELASTIC MATERIALS

      subroutine dstrq(coords,nxt,nyt,nzt,dh,dt,nz)

c     4th order finite-difference of stress components

c     nxt   nodal points in x dir          (integer)(sent)
c     nyt   nodal points in y dir          (integer)(sent)
c     nzt   nodal points in z dir          (integer)(sent)
c     dh    spatial discretization         (real)   (sent)
c     dt    temporal discretization        (real)   (sent)

      use parstat
      implicit none

      integer, dimension(3) :: coords
      integer :: nxt,nyt,nzt,nz
      real :: dh,dt

      call allq(xi(18),xf(18),yi(18),yf(18),zi(18),zf(18),
     +          coords,nxt,nyt,nzt,dh,dt,nz)

      return
      end

      subroutine dstrqc(coords,nxt,nyt,nzt,dh,dt,nz,rank)

c     4th order finite-difference of stress components

c     nxt   nodal points in x dir          (integer)(sent)
c     nyt   nodal points in y dir          (integer)(sent)
c     nzt   nodal points in z dir          (integer)(sent)
c     dh    spatial discretization         (real)   (sent)
c     dt    temporal discretization        (real)   (sent)

      use parstat
      implicit none

      integer, dimension(3) :: coords
      integer :: nxt,nyt,nzt,nz,rank
      real :: dh,dt

      call allq(1,nxt,1,nyt,1,nzt,coords,nxt,nyt,nzt,dh,dt,nz)

      return
      end

      subroutine allq(nxb,nxe,nyb,nye,nzb,nze,coords,nxt,nyt,nzt,dh,dt,nz)

c     4th order finite-difference of normal stresses at t+1

c     nxb   starting point for FD in x dir (integer)(sent)
c     nxe   ending point for FD in x dir   (integer)(sent)
c     nyb   starting point for FD in y dir (integer)(sent)
c     nye   ending point for FD in y dir   (integer)(sent)
c     nzb   starting point for FD in z dir (integer)(sent)
c     nze   ending point for FD in z dir   (integer)(sent)
c     dh    spatial discretization         (real)   (sent)
c     dt    temporal discretization        (real)   (sent)

      use parstat
      implicit none

      integer, dimension(3) :: coords
      integer :: nxt,nyt,nzt,nz
      integer :: kk,jj,i,j,k,itx,ity,itz,iii,jjj,kkk,num
      integer :: nxb,nxe,nyb,nye,nzb,nze
      real    :: dh1,dth,dh,dt,c1,c2,dt1,tauu
      real, dimension(nxb:nxe) :: vxl,vxm,vqpa,vqsa
      real, dimension(nxb:nxe) :: vx1,vx2,vt1,vt2,vt3,vt4,vt5,vt6,vt7,va1,vh1,vh2_r4,vh2_r5,vh2_r6,vxmu,vaas1,vtmp,vaas2,vaas3
      real, dimension(nxb:nxe) :: ww,weights_p,weights_s,weights_s1,weights_s2,weights_s3
      real, dimension (2,2,2) :: tau1,tau2

      dh1 = 1./dh
      dth = dt/dh
      c1 = 9./8.
      c2 = -1./24.

      dt1 = 1./dt
      do k=1,2
        do j=1,2
           do i=1,2
              tauu=tau(i,j,k)
c     old way: no exponnetials
c              tau1(i,j,k)=1/((tauu*dt1)+(1./2.))
c              tau2(i,j,k)=(tauu*dt1)-(1./2.)

              tau2(i,j,k)=exp(-dt/tauu)
              tau1(i,j,k)=0.5*(1.-tau2(i,j,k))
           enddo
        enddo
      enddo

      do kk= nzb,nze,kblock
        do jj= nyb,nye,jblock
          do 50 k= kk,min(kk+kblock-1,nze)
            itz = mod(((nz+1)-(coords(3)*nzt+k)),2)+1
              do 50 j= jj,min(jj+jblock-1,nye)
                ity = mod((coords(2)*nyt+j),2)+1
                itx = mod((coords(1)*nxt+nxb),2)+1
c                do i=nxb,nxe
c                  vx1(i)=tau1(itx,ity,itz)
c                  vx2(i)=tau2(itx,ity,itz)
c                  ww(i)=8.*weights(itx,ity,itz)
c                  itx = 3 - itx
c                end do

C == ADDED BY LUIS FOR SGSN DYNAMIC FAULT MODEL
                vxl(nxb:nxe)=8./(lam(nxb:nxe,j,k)+lam(nxb+1:nxe+1,j,k) +
     +               lam(nxb:nxe,j-1,k)+lam(nxb+1:nxe+1,j-1,k) +
     +               lam(nxb:nxe,j,k-1)+lam(nxb+1:nxe+1,j,k-1) +
     +               lam(nxb:nxe,j-1,k-1)+lam(nxb+1:nxe+1,j-1,k-1))

                vxm(nxb:nxe)=8./(mu(nxb:nxe,j,k)+mu(nxb+1:nxe+1,j,k) +
     +               mu(nxb:nxe,j-1,k)+mu(nxb+1:nxe+1,j-1,k) +
     +               mu(nxb:nxe,j,k-1)+mu(nxb+1:nxe+1,j,k-1) +
     +               mu(nxb:nxe,j-1,k-1)+mu(nxb+1:nxe+1,j-1,k-1))

                vxl(nxb:nxe)=vxl(nxb:nxe)+2.*vxm(nxb:nxe)

                vqpa(nxb:nxe)=0.125*(qp(nxb:nxe,j,k)+qp(nxb+1:nxe+1,j,k) +
     +                 qp(nxb:nxe,j-1,k)+qp(nxb+1:nxe+1,j-1,k) +
     +                 qp(nxb:nxe,j,k-1)+qp(nxb+1:nxe+1,j,k-1) +
     +                 qp(nxb:nxe,j-1,k-1)+qp(nxb+1:nxe+1,j-1,k-1))

                vqsa(nxb:nxe)=0.125*(qs(nxb:nxe,j,k)+qs(nxb+1:nxe+1,j,k) +
     +                 qs(nxb:nxe,j-1,k)+qs(nxb+1:nxe+1,j-1,k) +
     +                 qs(nxb:nxe,j,k-1)+qs(nxb+1:nxe+1,j,k-1) +
     +               qs(nxb:nxe,j-1,k-1)+qs(nxb+1:nxe+1,j-1,k-1))
                vaas1(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j,k-1))
                vaas2(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j-1,k))
                vaas3(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb+1:nxe+1,j,k))

                do i=nxb,nxe
                  vx1(i)=tau1(itx,ity,itz)
                  vx2(i)=tau2(itx,ity,itz)
c                  ww(i)=8.*weights(itx,ity,itz)

              if (1./vqsa(i) .lt. 200.0) then
c                 print *, 9999999999999
              num=0
              do iii=1,2
                 do jjj=1,2
                    do kkk=1,2
                       weights_los(iii,jjj,kkk)=coeff(1+num)*(vqsa(i)**2)+coeff(2+num)*vqsa(i)
c                    weights_los(iii,jjj,kkk)=coeff(1+num)*(1./vqsa(i))**4.+coeff(2+num)*(1./vqsa(i))**3.+
c     +                     coeff(3+num)*(1./vqsa(i))**2.+coeff(4+num)*(1./vqsa(i))+coeff(5+num);
                     num=num+2
                     end do
                     end do
                     end do

                     weights_s(i)=weights_los(itx,ity,itz)
                  else
                     weights_s(i)=8.*weights(itx,ity,itz)*vqsa(i)
                     endif

              if (1./vqpa(i) .lt. 200.0) then
              num=0
              do iii=1,2
                 do jjj=1,2
                    do kkk=1,2
                       weights_lop(iii,jjj,kkk)=coeff(1+num)*(vqpa(i)**2)+coeff(2+num)*vqpa(i)
                     num=num+2
                     end do
                     end do
                     end do

                     weights_p(i)=weights_lop(itx,ity,itz)
                  else
                     weights_p(i)=8.*weights(itx,ity,itz)*vqpa(i)
                     endif
cccccccccccccccccccccccccccccccc
              if (1./vaas1(i) .lt. 200.0) then
              num=0
              do iii=1,2
                 do jjj=1,2
                    do kkk=1,2
                       weights_los(iii,jjj,kkk)=coeff(1+num)*(vaas1(i)**2)+coeff(2+num)*vaas1(i)
                     num=num+2
                     end do
                     end do
                     end do

c                     weights_s1(i)=weights_los(itx,ity,itz)/vaas(i)
                     weights_s1(i)=weights_los(itx,ity,itz)
                  else
                     weights_s1(i)=8.*weights(itx,ity,itz)*vaas1(i)
                     endif
ccccccccccccccccccccccccccccccccccccccccc
              if (1./vaas2(i) .lt. 200.0) then
              num=0
              do iii=1,2
                 do jjj=1,2
                    do kkk=1,2
                       weights_los(iii,jjj,kkk)=coeff(1+num)*(vaas2(i)**2)+coeff(2+num)*vaas2(i)
                     num=num+2
                     end do
                     end do
                     end do

c                     weights_s1(i)=weights_los(itx,ity,itz)/vaas(i)
                     weights_s2(i)=weights_los(itx,ity,itz)
                  else
                     weights_s2(i)=8.*weights(itx,ity,itz)*vaas2(i)
                     endif
ccccccccccccccccccccccccccccccccccccccccc
              if (1./vaas3(i) .lt. 200.0) then
              num=0
              do iii=1,2
                 do jjj=1,2
                    do kkk=1,2
                       weights_los(iii,jjj,kkk)=coeff(1+num)*(vaas3(i)**2)+coeff(2+num)*vaas3(i)
                     num=num+2
                     end do
                     end do
                     end do

c                     weights_s1(i)=weights_los(itx,ity,itz)/vaas(i)
                     weights_s3(i)=weights_los(itx,ity,itz)
                  else
                     weights_s3(i)=8.*weights(itx,ity,itz)*vaas3(i)
                     endif
ccccccccccccccccccccccccccccccccccccccccc

                  itx = 3 - itx
                end do

                vt1(nxb:nxe)=c1*(u1(nxb+1:nxe+1,j,k)-u1(nxb:nxe,j,k))+c2*(u1(nxb+2:nxe+2,j,k)-u1(nxb-1:nxe-1,j,k))
                vt2(nxb:nxe)=c1*(v1(nxb:nxe,j,k)-v1(nxb:nxe,j-1,k))+c2*(v1(nxb:nxe,j+1,k)-v1(nxb:nxe,j-2,k))
                vt3(nxb:nxe)=c1*(w1(nxb:nxe,j,k)-w1(nxb:nxe,j,k-1))+c2*(w1(nxb:nxe,j,k+1)-w1(nxb:nxe,j,k-2))
                vt4(nxb:nxe)=c1*(u1(nxb:nxe,j+1,k)-u1(nxb:nxe,j,k))+c2*(u1(nxb:nxe,j+2,k)-u1(nxb:nxe,j-1,k)) +
     $                 c1*(v1(nxb:nxe,j,k)-v1(nxb-1:nxe-1,j,k))+c2*(v1(nxb+1:nxe+1,j,k)-v1(nxb-2:nxe-2,j,k))
                vt5(nxb:nxe)=c1*(u1(nxb:nxe,j,k+1)-u1(nxb:nxe,j,k))+c2*(u1(nxb:nxe,j,k+2)-u1(nxb:nxe,j,k-1)) +
     $                 c1*(w1(nxb:nxe,j,k)-w1(nxb-1:nxe-1,j,k))+c2*(w1(nxb+1:nxe+1,j,k)-w1(nxb-2:nxe-2,j,k))
                vt6(nxb:nxe)=c1*(v1(nxb:nxe,j,k+1)-v1(nxb:nxe,j,k))+c2*(v1(nxb:nxe,j,k+2)-v1(nxb:nxe,j,k-1))
                vt7(nxb:nxe)=c1*(w1(nxb:nxe,j+1,k)-w1(nxb:nxe,j,k))+c2*(w1(nxb:nxe,j+2,k)-w1(nxb:nxe,j-1,k))

c kbwapril15                vh1(nxb:nxe)=-vxm(nxb:nxe)*vqsa(nxb:nxe)
                vh1(nxb:nxe)=-vxm(nxb:nxe)
                vtmp(nxb:nxe)=vxl(nxb:nxe)*(vt1(nxb:nxe)+vt2(nxb:nxe)+vt3(nxb:nxe))
                !va1(nxb:nxe)=-vqpa(nxb:nxe)*dh1*vtmp(nxb:nxe)/2.
c     KBWapril15                va1(nxb:nxe)=-vqpa(nxb:nxe)*vxl(nxb:nxe)*dh1*(vt1(nxb:nxe)+vt2(nxb:nxe)+vt3(nxb:nxe))/2.
                va1(nxb:nxe)=-vxl(nxb:nxe)*dh1*(vt1(nxb:nxe)+vt2(nxb:nxe)+vt3(nxb:nxe))/2.

                !vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb:nxe,j,k-1))
                !vh2(nxb:nxe)=-vxmu(nxb:nxe,j,k)*vaas(nxb:nxe,j,k)/2.
                !vaas(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j,k-1))

C   NORMAL STRESS XX, YY AND ZZ

                xx(nxb:nxe,j,k)=xx(nxb:nxe,j,k)+dth*(vtmp(nxb:nxe)-2.*vxm*(vt2(nxb:nxe)+vt3(nxb:nxe)))+dt*r1(nxb:nxe,j,k)
                yy(nxb:nxe,j,k)=yy(nxb:nxe,j,k)+dth*(vtmp(nxb:nxe)-2.*vxm*(vt1(nxb:nxe)+vt3(nxb:nxe)))+dt*r2(nxb:nxe,j,k)
                zz(nxb:nxe,j,k)=zz(nxb:nxe,j,k)+dth*(vtmp(nxb:nxe)-2.*vxm*(vt1(nxb:nxe)+vt2(nxb:nxe)))+dt*r3(nxb:nxe,j,k)

                !vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb:nxe,j,k-1))
                !vh2_r4(nxb:nxe)=-0.25*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j,k-1))*vxmu(nxb:nxe)
                vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb:nxe,j,k-1))
c1                vaas(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j,k-1))
c                vh2_r4(nxb:nxe)=-vaas(nxb:nxe)*vxmu(nxb:nxe)/2.
                 vh2_r4(nxb:nxe)=-vxmu(nxb:nxe)/2.
                xy(nxb:nxe,j,k)=xy(nxb:nxe,j,k)+vxmu(nxb:nxe)*dth*vt4(nxb:nxe)+dt*r4(nxb:nxe,j,k)
                !vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb:nxe,j-1,k))
                !vh2_r5(nxb:nxe)=-0.25*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j-1,k))*vxmu(nxb:nxe)
                vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb:nxe,j-1,k))
c2                vaas(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb:nxe,j-1,k))
c                vh2_r5(nxb:nxe)=-vaas(nxb:nxe)*vxmu(nxb:nxe)/2.
                vh2_r5(nxb:nxe)=-vxmu(nxb:nxe)/2.
                xz(nxb:nxe,j,k)=xz(nxb:nxe,j,k)+vxmu(nxb:nxe)*dth*vt5(nxb:nxe)+dt*r5(nxb:nxe,j,k)
                !vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb+1:nxe+1,j,k))
                !vh2_r6(nxb:nxe)=-0.25*(qs(nxb:nxe,j,k)+qs(nxb+1:nxe+1,j,k))*vxmu(nxb:nxe)
                vxmu(nxb:nxe) = 2./(mu(nxb:nxe,j,k)+mu(nxb+1:nxe+1,j,k))
c3                vaas(nxb:nxe) = 0.5*(qs(nxb:nxe,j,k)+qs(nxb+1:nxe+1,j,k))
c                vh2_r6(nxb:nxe)=-vaas(nxb:nxe)*vxmu(nxb:nxe)/2.
                vh2_r6(nxb:nxe)=-vxmu(nxb:nxe)/2.
                yz(nxb:nxe,j,k)=yz(nxb:nxe,j,k)+vxmu(nxb:nxe)*(dth*vt6(nxb:nxe)+dth*vt7(nxb:nxe))+dt*r6(nxb:nxe,j,k)

                r1(nxb:nxe,j,k)=vx2(nxb:nxe)*r1(nxb:nxe,j,k)-
     + (weights_s(nxb:nxe)*vh1(nxb:nxe)*(vt2(nxb:nxe)+ vt3(nxb:nxe))*dh1-
     + weights_p(nxb:nxe)*va1(nxb:nxe))*vx1(nxb:nxe)*2.
                r2(nxb:nxe,j,k)=vx2(nxb:nxe)*r2(nxb:nxe,j,k)-
     + (weights_s(nxb:nxe)*vh1(nxb:nxe)*(vt1(nxb:nxe)+vt3(nxb:nxe))*dh1-
     + weights_p(nxb:nxe)*va1(nxb:nxe))*vx1(nxb:nxe)*2.
                r3(nxb:nxe,j,k)=vx2(nxb:nxe)*r3(nxb:nxe,j,k)-
     + (weights_s(nxb:nxe)*vh1(nxb:nxe)*(vt1(nxb:nxe)+vt2(nxb:nxe))*dh1-
     + weights_p(nxb:nxe)*va1(nxb:nxe))*vx1(nxb:nxe)*2.
                r4(nxb:nxe,j,k)=vx2(nxb:nxe)*r4(nxb:nxe,j,k)+
     + weights_s1(nxb:nxe)*vh2_r4(nxb:nxe)*dh1*vt4(nxb:nxe)*vx1(nxb:nxe)*2.
                r5(nxb:nxe,j,k)=vx2(nxb:nxe)*r5(nxb:nxe,j,k)+
     + weights_s2(nxb:nxe)*vh2_r5(nxb:nxe)*dh1*vt5(nxb:nxe)*vx1(nxb:nxe)*2.
                r6(nxb:nxe,j,k)=vx2(nxb:nxe)*r6(nxb:nxe,j,k)+
     + weights_s3(nxb:nxe)*vh2_r6(nxb:nxe)*dh1*(vt6(nxb:nxe)+vt7(nxb:nxe))*
     + vx1(nxb:nxe)*2.

c                r1(nxb:nxe,j,k)=vx2(nxb:nxe)*r1(nxb:nxe,j,k)-
c     + (weights_s(nxb:nxe)*(vh1(nxb:nxe)*(vt2(nxb:nxe)+ vt3(nxb:nxe))*dh1-
c     + va1(nxb:nxe)))*vx1(nxb:nxe)*2.
c                r2(nxb:nxe,j,k)=vx2(nxb:nxe)*r2(nxb:nxe,j,k)-
c     + (weights_s(nxb:nxe)*(vh1(nxb:nxe)*(vt1(nxb:nxe)+vt3(nxb:nxe))*dh1-
c     + va1(nxb:nxe)))*vx1(nxb:nxe)*2.
c                r3(nxb:nxe,j,k)=vx2(nxb:nxe)*r3(nxb:nxe,j,k)-
c     + (weights_s(nxb:nxe)*(vh1(nxb:nxe)*(vt1(nxb:nxe)+vt2(nxb:nxe))*dh1-
c     + va1(nxb:nxe)))*vx1(nxb:nxe)*2.
c                r4(nxb:nxe,j,k)=vx2(nxb:nxe)*r4(nxb:nxe,j,k)+
c     + (weights_s(nxb:nxe)*vh2_r4(nxb:nxe)*dh1*vt4(nxb:nxe))*vx1(nxb:nxe)*2.
c                r5(nxb:nxe,j,k)=vx2(nxb:nxe)*r5(nxb:nxe,j,k)+
c     + (weights_s(nxb:nxe)*vh2_r5(nxb:nxe)*dh1*vt5(nxb:nxe))*vx1(nxb:nxe)*2.
c                r6(nxb:nxe,j,k)=vx2(nxb:nxe)*r6(nxb:nxe,j,k)+
c     + (weights_s(nxb:nxe)*vh2_r6(nxb:nxe)*dh1*(vt6(nxb:nxe)+vt7(nxb:nxe)))*
c     + vx1(nxb:nxe)*2.

                xx(nxb:nxe,j,k)=xx(nxb:nxe,j,k)+dt*r1(nxb:nxe,j,k)
                yy(nxb:nxe,j,k)=yy(nxb:nxe,j,k)+dt*r2(nxb:nxe,j,k)
                zz(nxb:nxe,j,k)=zz(nxb:nxe,j,k)+dt*r3(nxb:nxe,j,k)
                xy(nxb:nxe,j,k)=xy(nxb:nxe,j,k)+dt*r4(nxb:nxe,j,k)
                xz(nxb:nxe,j,k)=xz(nxb:nxe,j,k)+dt*r5(nxb:nxe,j,k)
                yz(nxb:nxe,j,k)=yz(nxb:nxe,j,k)+dt*r6(nxb:nxe,j,k)
50        continue
        enddo
      enddo

      return
      end


!
!
!                             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.
!
!
subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)

  use psb_base_mod
  use amg_s_ilu_solver, amg_protect_name => amg_s_ilu_solver_bld

  Implicit None

  ! Arguments
  type(psb_sspmat_type), intent(in), target           :: a
  Type(psb_desc_type), Intent(inout)                  :: desc_a
  class(amg_s_ilu_solver_type), intent(inout)         :: sv
  integer(psb_ipk_), intent(out)                      :: info
  type(psb_sspmat_type), intent(in), target, optional :: b
  class(psb_s_base_sparse_mat), intent(in), optional  :: amold
  class(psb_s_base_vect_type), intent(in), optional   :: vmold
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
!!$    real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
  character(len=20)   :: name='s_ilu_solver_bld', ch_err

  info=psb_success_
  call psb_erractionsave(err_act)
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()
  ctxt       = desc_a%get_context()
  call psb_info(ctxt, me, np)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' start'

  n_row  = desc_a%get_local_rows()

  nrow_a = a%get_nrows()
  nztota = a%get_nzeros()
  if (present(b)) then
    nztota = nztota + b%get_nzeros()
  end if

  call sv%l%csall(n_row,n_row,info,nztota)
  if (info == psb_success_) call sv%u%csall(n_row,n_row,info,nztota)
  if(info /= psb_success_) then
    info=psb_err_from_subroutine_
    ch_err='psb_sp_all'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  if (allocated(sv%d)) then
    if (size(sv%d) < n_row) then
      deallocate(sv%d)
    endif
  endif
  if (.not.allocated(sv%d))  allocate(sv%d(n_row),stat=info)

  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
    goto 9999
  endif

  select case(sv%fact_type)

  case (psb_ilu_t_)
    !
    ! ILU(k,t)
    !
    select case(sv%fill_in)

    case(:-1)
      ! Error: fill-in <= -1
      call psb_errpush(psb_err_input_value_invalid_i_,&
           & name,i_err=(/ithree,sv%fill_in,izero,izero,izero/))
      goto 9999

    case(0:)
      ! Fill-in >= 0
      call psb_ilut_fact(sv%fill_in,sv%thresh,&
           & a, sv%l,sv%u,sv%d,info,blck=b)
    end select
    if(info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='psb_ilut_fact'
      call psb_errpush(info,name,a_err=ch_err)
      goto 9999
    end if

  case(psb_ilu_n_,psb_milu_n_)
    !
    ! ILU(k) and MILU(k)
    !
    select case(sv%fill_in)
    case(:-1)
      ! Error: fill-in <= -1
      call psb_errpush(psb_err_input_value_invalid_i_,&
           & name,i_err=(/ithree,sv%fill_in,izero,izero,izero/))
      goto 9999
    case(0)
      ! Fill-in 0
      ! Separate implementation of ILU(0) for better performance.
      ! There seems to be a problem with the separate implementation of MILU(0),
      ! contained into psb_ilu0_fact. This must be investigated. For the time being,
      ! resort to the implementation of MILU(k) with k=0.
      if (sv%fact_type == psb_ilu_n_) then
        call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,&
             & sv%d,info,blck=b)
      else
        call psb_iluk_fact(sv%fill_in,sv%fact_type,&
             & a,sv%l,sv%u,sv%d,info,blck=b)
      endif
    case(1:)
      ! Fill-in >= 1
      ! The same routine implements both ILU(k) and MILU(k)
      call psb_iluk_fact(sv%fill_in,sv%fact_type,&
           & a,sv%l,sv%u,sv%d,info,blck=b)
    end select
    if (info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='psb_iluk_fact'
      call psb_errpush(info,name,a_err=ch_err)
      goto 9999
    end if

  case default
    ! If we end up here, something was wrong up in the call chain.
    info = psb_err_input_value_invalid_i_
    call psb_errpush(psb_err_input_value_invalid_i_,name,&
         & i_err=(/ithree,sv%fact_type,izero,izero,izero/))
    goto 9999

  end select

  call sv%l%set_asb()
  call sv%l%trim()
  call sv%u%set_asb()
  call sv%u%trim()
  call sv%dv%bld(sv%d,mold=vmold)

  if (present(amold)) then
    call sv%l%cscnv(info,mold=amold)
    call sv%u%cscnv(info,mold=amold)
  end if

  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' end'

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return
end subroutine amg_s_ilu_solver_bld


text:
      REAL FUNCTION AIS (NP,I,L,RR,ZZ)
C
C     THIS ROUTINE CALCULATES THE SINGLE PRECISION  DELTA(IJ) INTEGRALS
C     FOR AXISYMMETRIC SOLIDS IN SDR2. CALCULATIONS DONE IN DOUBLE
C     PRECISION
C
C     INPUT
C     NP  = NUMBER OF POINTS (3 OR 4.  MORE THAN 4 WILL FAIL FSN-5)
C     I,L = THE INTEGRAL DESIRED (I SERIES STARTS WITH -1)
C     R   = RADIUS ARRAY (NP LONG)
C     Z   = Z-CORD ARRAY (NP LONG)
C
C     OUTPUT
C     DKL = DESIRED INTEGRAL
C
      INTEGER          NAM(2)
      DOUBLE PRECISION A,AJ,AR,BETA,DR,DZ,DFACT,DZJ,DL1,DKL,EPS,FACT,
     1                 GKL,ZERO,ONE,TWO,THREE,PR,RA,RB,R(4),RAK,RBK,
     2                 ZA,ZB,Z(4)
      DIMENSION        RR(4),ZZ(4)
      DATA       EPS / .01D0 /
      DATA       ZERO, ONE,TWO,THREE / 0.0D0, 1.0D0, 2.0D0, 3.0D0 /
      DATA       NAM / 4HAIS , 1H    /
C
      DO 5 M = 1,NP
      R(M) = RR(M)
    5 Z(M) = ZZ(M)
      DKL  = ZERO
      L1   = L + 1
      L2   = L + 2
      DL1  = L1
      K    = I + 1
C
C     LOOP ON NUMBER OF POINTS
C
      IF (R(1) .LE. ZERO) GO TO 300
      DO 200 M = 1,NP
      J  = M + 1
      IF (M .EQ. NP) J = 1
      RA = R(M)
      RB = R(J)
      ZA = Z(M)
      ZB = Z(J)
      DR = RB - RA
      DZ = ZB - ZA
C
C     TEST IF RADIUS IS .LE. 0 (DRIVER SHOULD FIND THIS)
C
      IF (RB .LE. ZERO) GO TO 300
      GKL = ZERO
      PR  = RA + RB
      AR  = PR/TWO
C
C     CHECK FOR APPROXIMATION, DR/AVE(R)
C
      IF ( DABS ( DR/AR ) .LT. EPS ) GO TO 70
C
      A = ZA*DR - RA*DZ
      BETA = A/DR
C
C     CHECK FOR BETA .EQ. 0 CASE
C
      IF (DABS(BETA/AR) .GT. EPS) GO TO 10
C
      IF (DZ .EQ. ZERO) GO TO 200
      LK  = L + K + 1
      AR  = LK
      GKL = (DZ/DR)**L1*(RA**LK-RB**LK)/(DL1*AR)
      GO TO 200
C
C     GENERAL CASE
C
   10 RAK = RA**K
      RBK = RB**K
      IF (K) 300,20,30
C
C     GENERAL CASE, K.EQ.0, CONSTANT TERM
C
   20 GKL = DLOG(RA/RB)/DL1
      GO TO 40
C
C     GENERAL CASE, CONSTANT TERM
C
   30 AR  = K*L1
      GKL = (RAK-RBK)/AR
C
C     GENERAL CASE, SUMMATION
C
   40 IF (DZ .EQ. ZERO) GO TO 65
      LFACT = 1
C
C     CALCULATE FACTORIAL (L+1)
C
      DO 50 J = 2,L
   50 LFACT = LFACT*J
      FACTL = LFACT
      JFACT = 1
      AJ    = ONE
      DZJ   = ONE
      LMJF  = LFACT*L1
      DO 60 J = 1,L1
      JFACT = JFACT*J
C
C     CALCULATE (L+1-J) FACTORIAL IN LMJF
C
      LMJF  = LMJF /(L2-J)
      FACT  = FACTL/FLOAT (JFACT*LMJF)
      DFACT = K + J
      DFACT = FACT /DFACT
      AJ    = AJ *A
      RAK   = RAK*RA
      RBK   = RBK*RB
      DZJ   = DZJ*DZ
   60 GKL   = GKL + (DFACT*DZJ*(RAK-RBK))/AJ
C
   65 GKL   = GKL*BETA**L1
      GO TO 200
C
C     APPROXIMATE CODE
C
   70 CONTINUE
      IF (DR .EQ. ZERO) GO TO 200
      DZJ = L1*L2
      RBK = ZB**L1
      J   = K - 1
      GKL =-DR*AR**J*RBK/DL1
C
      IF (DZ .EQ. ZERO) GO TO 200
      GKL = GKL + (((TWO*RA+RB)/THREE)**J*DR*DABS(ZA**L2-RBK*ZB))
     1            /(DZJ*DZ)
C
  200 DKL = DKL + GKL
C
C     ALL DONE
C
  210 AIS = DKL
      RETURN
C
C     ERROR
C
  300 CALL MESAGE (-7,K,NAM)
      GO TO 210
      END

ext: f


text:
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPHERP     SPHERICAL TRANSFORM TO A POINT
C   PRGMMR: KANAMITSU     ORG: NP23          DATE: 1992-04-16
C
C ABSTRACT: TRANSFORMS A SPECTRAL FIELD TO A POINT ON THE GLOBE.
C
C PROGRAM HISTORY LOG:
C 1992-04-16  IREDELL
C
C USAGE:    CALL SPHERP(WAVE,MAXWV,IROMB,IDY,MLT,FAC,RLAT,RLON,VALUE)
C
C   INPUT ARGUMENT LIST:
C     WAVE     - IF IDIR < 0, COMPLEX (KMAX) FIELD TO TRANSFORM,
C                WHERE KMAX=(MAXWV+1)*(IROMB+1)*MAXWV+2)/2.
C                WAVE STARTS AT THE GLOBAL MEAN COMPONENT AND THEN
C                CONTAINS THE ZONALLY SYMMETRIC COMPONENTS.
C     MAXWV    - SPECTRAL TRUNCATION OF THE WAVE
C     IROMB    - IROMB = 0 FOR TRIANGULAR TRUNCATION
C                IROMB = 1 FOR RHOMBOIDAL TRUNCATION
C     IDIR     - 0 FOR NORMAL, 1 FOR LEGENDRE DERIVATIVE
C     MLT      - MULTIPLICATION OPTION ON WAVE VALUES.
C                MLT = 0 FOR NO MULTIPLICATION,
C                MLT = 1 TO MULTIPLY WAVE BY FAC
C                MLT = -1 TO MULTIPLY WAVE BY SQRT(-1)*FAC
C     FAC      - IF MLT.NE.0, REAL (KMAX) FIELD TO MULTIPLY WAVE.
C     RLAT     - LATITUDE IN DEGREES NORTH
C     RLON     - LONGITUDE IN DEGREES EAST
C
C   OUTPUT ARGUMENT LIST:
C     VALUE    - TRANSFORMED VALUE
C
C   SUBPROGRAMS CALLED:
C     (LEGENDR2) - COMPUTE LEGENDRE POLYNOMIALS
C
C   REMARKS: FORTRAN 9X EXTENSIONS ARE USED.
C     FPP CAN BE USED TO INLINE SUBPROGRAMS.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM-SP
C
C$$$
CFPP$ EXPAND(LEGENDR2)
      SUBROUTINE SPHERP(WAVE,MAXWV,IROMB,IDY,MLT,FAC,RLAT,RLON,VALUE)
      REAL WAVE((MAXWV+1)*((IROMB+1)*MAXWV+2))
      REAL FAC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL PNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL DNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL FR(0:MAXWV),FI(0:MAXWV)
      PARAMETER(PI=3.1415926535898,PI180=PI/180.)
      SINLAT=SIN(RLAT*PI180)
      CALL LEGENDR2(SINLAT,MAXWV,IROMB,PNM,DNM)
      IF(IDY.NE.0) THEN
        PNM=DNM
      ENDIF
      IF(MLT.NE.0) THEN
        PNM=PNM*FAC
      ENDIF
      K=0
      DO M=0,MAXWV
        FR(M)=0.
        FI(M)=0.
      ENDDO
      DO M=0,MAXWV
        DO N=M,IROMB*M+MAXWV
          FR(M)=FR(M)+WAVE(2*K+1)*PNM(K+1)
          FI(M)=FI(M)+WAVE(2*K+2)*PNM(K+1)
          K=K+1
        ENDDO
      ENDDO
      VALUE=0.
      RLAM=RLON*PI180
      IF(MLT.GE.0) THEN
        DO M=1,MAXWV
          VALUE=VALUE+FR(M)*COS(RLAM*M)-FI(M)*SIN(RLAM*M)
        ENDDO
      ELSE
        DO M=1,MAXWV
          VALUE=VALUE-FI(M)*COS(RLAM*M)-FR(M)*SIN(RLAM*M)
        ENDDO
      ENDIF
      VALUE=FR(0)+2.*VALUE
      RETURN
      END
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  LEGENDR2   COMPUTE LEGENDRE POLYNOMIALS.
C   PRGMMR: IREDELL       ORG: NP23          DATE: 1992-04-16
C
C ABSTRACT: EVALUATES THE NORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
C   AT A GIVEN COSINE OF COLATITUDE.
C
C PROGRAM HISTORY LOG:
C 1992-04-16  IREDELL
C
C USAGE:    CALL LEGENDR2(COSCLT,MAXWV,IROMB,PNM,DNM)
C
C   INPUT ARGUMENT LIST:
C     COSCLT   - COSINE OF COLATITUDE AT WHICH TO EVALUTATE.
C     MAXWV    - SPECTRAL TRUNCATION
C     IROMB    - IROMB = 0 FOR TRIANGULAR TRUNCATION
C                IROMB = 1 FOR RHOMBOIDAL TRUNCATION
C
C   OUTPUT ARGUMENT LIST:
C     PNM      - REAL (KMAX) POLYNOMIAL VALUES,
C     DNM      - REAL (KMAX) POLYNOMIAL DERIVATIVES,
C                WHERE KMAX=(MAXWV+1)*(IROMB+1)*MAXWV+2)/2.
C
C   REMARKS: FORTRAN 9X EXTENSIONS ARE USED.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM-SP
C
C$$$
      SUBROUTINE LEGENDR2(COSCLT,MAXWV,IROMB,PNM,DNM)
      DIMENSION PNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      DIMENSION DNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      DIMENSION EPS(0:100)
      SINCLT=SQRT(1.-COSCLT**2)
      NM=0
      K=1
      PNM(K)=SQRT(0.5)
CDIR$ NEXTSCALAR
      DO M=1,MAXWV-NM*(1-IROMB)
        PNMK=PNM(K)
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        PNM(K)=PNMK*SINCLT*SQRT(FLOAT(2*M+1)/FLOAT(2*M))
      ENDDO
      NM=1
CDIR$ IVDEP
      DO M=0,MAXWV-NM*(1-IROMB)
        N=M+NM
        EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        PNM(K)=COSCLT*PNM(K-1)/EPSNM
        DNM(K-1)=EPSNM*PNM(K)*(N-1)
        EPS(M)=EPSNM
      ENDDO
      DO NM=2,MAXWV
CDIR$ IVDEP
        DO M=0,MAXWV-NM*(1-IROMB)
          N=M+NM
          EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
          K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
          PNM(K)=(COSCLT*PNM(K-1)-EPS(M)*PNM(K-2))/EPSNM
          DNM(K-1)=EPSNM*PNM(K)*(N-1)-EPS(M)*PNM(K-2)*N
          EPS(M)=EPSNM
        ENDDO
      ENDDO
CDIR$ IVDEP
      DO M=0,MAXWV
        NM=MAXWV-M*(1-IROMB)
        N=M+NM
        EPS(M)=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        N=N+1
        EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        K=K+1
        PNMK=(COSCLT*PNM(K-1)-EPS(M)*PNM(K-2))/EPSNM
        DNM(K-1)=EPSNM*PNMK*(N-1)-EPS(M)*PNM(K-2)*N
      ENDDO
      RETURN
      END

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
module transport_properties

  use iso_c_binding
  use chemistry_module
  use variables_module

  implicit none

  double precision, parameter :: Pr = 0.72d0   ! Prandtl number
  double precision, parameter :: Sc = 0.72d0   ! Schmidt number
  double precision, parameter :: C_S = 1.458d-5  ! constant in Sutherland's law
  double precision, parameter :: T_S = 110.4d0   ! Sutherland temperature

  private

  public get_trans_prop_3d

contains

  subroutine get_trans_prop_3d(tlo,thi,lo,hi,q,mu,xi,lam,Ddiag,ng) bind(c,name='get_trans_prop_3d')
    integer, intent(in) :: lo(3), hi(3), ng, tlo(3), thi(3)
    double precision,intent(in   )::    q(lo(1)-ng:hi(1)+ng,lo(2)-ng:hi(2)+ng,lo(3)-ng:hi(3)+ng,nprim)
    double precision,intent(inout)::   mu(lo(1)-ng:hi(1)+ng,lo(2)-ng:hi(2)+ng,lo(3)-ng:hi(3)+ng)
    double precision,intent(inout)::   xi(lo(1)-ng:hi(1)+ng,lo(2)-ng:hi(2)+ng,lo(3)-ng:hi(3)+ng)
    double precision,intent(inout)::  lam(lo(1)-ng:hi(1)+ng,lo(2)-ng:hi(2)+ng,lo(3)-ng:hi(3)+ng)
    double precision,intent(inout)::Ddiag(lo(1)-ng:hi(1)+ng,lo(2)-ng:hi(2)+ng,lo(3)-ng:hi(3)+ng,nspecies)

    integer :: i, j, k, n, iwrk
    double precision :: rwrk, cp, Wbar, Wbarinv, Yt(nspecies)
    double precision :: Prinv, Scinv

    Prinv = 1.d0/Pr
    Scinv = 1.d0/Sc

    do k=tlo(3),thi(3)
       do j=tlo(2),thi(2)
          do i=tlo(1),thi(1)

             mu(i,j,k) = C_S * sqrt(q(i,j,k,qtemp)) * q(i,j,k,qtemp) &
                  / (q(i,j,k,qtemp) + T_S)

             xi(i,j,k) = 0.d0

             do n=1,nspecies
                Yt(n) = q(i,j,k,qy1+n-1)
             end do

             call ckmmwy(Yt, iwrk, rwrk, Wbar)
             Wbarinv = 1.d0 / Wbar

             call ckcpbs(q(i,j,k,qtemp), Yt, iwrk, rwrk, cp)

             lam(i,j,k) = mu(i,j,k) * cp * Prinv

             do n=1,nspecies
                Ddiag(i,j,k,n) = mu(i,j,k) * molecular_weight(n) * Wbarinv * Scinv
             end do

          end do
       end do
    end do

  end subroutine get_trans_prop_3d

end module transport_properties

max_stars_repo_name: ylunalin/amrex

max_issues_repo_name: ylunalin/amrex

max_forks_repo_name: ylunalin/amrex


text:
Located just south of Academic Surge on La Rue, the Watershed Sciences building houses the eponymous UC Davis Center for Watershed Sciences, which in turn researches the ecological health of rivers, streams, and deltas in California. Interestingly, only the first floor of the building was devoted to this purpose. The second floor was left as an unfinished shell, with nothing but walls and electrical sockets.

Early in 2009, several labs and offices in the Computer Science department were moved from Kemper Hall to this building, and took up residence on the second floor. Laboratories to be found here include the Security Lab and the Networking Lab.

Misc. Photos

 Construction Photos

 Cracks

This relatively new building has a pretty serious pattern of cracks in the wiki:WikiPedia:Stucco stucco.

I have seen this happen, often because of it not drying properly suggesting a kinda off mixture for the stucco. They should have redone it, but then again they are on budget, hopefully these are not too terribly deep cracks. ~Davepoole

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


SUBROUTINE stats (stat_q)

! ----------------------------------------------------------------------
! SUBROUTINE: stats.f90
! ----------------------------------------------------------------------
! Purpose:
!  Computation of statistical qunatities based on a provided data series
!  Root Mean Square (RMS), Maximum, Minimum, Mean, Standard deviation
! ----------------------------------------------------------------------
! Input arguments:
! - array1:         Input allocatable array
! - mdl_array:		Module for the input/output allocatable arrays
!
! Output arguments:
! - stat_q:			Statistics quantities: RMS, Max, Min, Mean, sdev
! ----------------------------------------------------------------------
! Dr. Thomas D. Papanikolaou, Geoscience Australia              May 2016
! ----------------------------------------------------------------------

      USE mdl_precision
      USE mdl_arr
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
! IN

! OUT
      REAL (KIND = prec_d) :: stat_q(5)
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Local variables declaration
! ----------------------------------------------------------------------
      REAL (KIND = prec_d) :: min_val, max_val, RMS, Mean, sdev
      REAL (KIND = prec_d) :: sum_dx2, sum_dx, sum_dxx
      INTEGER (KIND = prec_int4) :: sz1, sz2, i
! ----------------------------------------------------------------------

![size_dx size2] = size(dx);
      sz1 = SIZE (array1,DIM=1)
      !sz2 = SIZE (array1,DIM=2)

      min_val = array1(1)
      max_val = array1(1)
      sum_dx2 = 0.0D0
      sum_dx = 0.0D0

      DO i = 1 , sz1
	     sum_dx2 = sum_dx2 + array1(i)**2
	     sum_dx = sum_dx + array1(i)
         if ( array1(i) < min_val ) then
		     min_val = array1(i)
		 end if
		 if ( array1(i) > max_val ) then
		     max_val = array1(i)
		 end if
      END DO

! RMS
      RMS = sqrt(sum_dx2 / sz1)

! Mean
      Mean = sum_dx / sz1

! ----------------------------------------------------------------------
! Standard Deviation
      sum_dxx = 0.0D0
      DO i = 1 , sz1
	     sum_dxx = sum_dxx + (array1(i) - Mean)**2
      END DO

      sdev = sqrt(sum_dxx / sz1)

      sdev = sqrt(sum_dxx / (sz1 - 1) )
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
      stat_q(1) = RMS
      stat_q(2) = max_val
      stat_q(3) = min_val
      stat_q(4) = Mean
      stat_q(5) = sdev
! ----------------------------------------------------------------------

END


text:
      SUBROUTINE BORT2(STR1,STR2)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    BORT2
C   PRGMMR: KEYSER           ORG: NP22       DATE: 2003-11-04
C
C ABSTRACT: THIS SUBROUTINE PRINTS (TO STDOUT) TWO GIVEN ERROR STRINGS
C   AND THEN CALLS BUFR ARCHIVE LIBRARY SUBROUTINE BORT_EXIT TO ABORT
C   THE APPLICATION PROGRAM CALLING THE BUFR ARCHIVE LIBRARY SOFTWARE.
C   IT IS SIMILAR TO BUFR ARCHIVE LIBRARY SUBROUTINE BORT, EXCEPT BORT
C   PRINTS ONLY ONE ERROR STRING.
C
C PROGRAM HISTORY LOG:
C 2003-11-04  D. KEYSER  -- ORIGINAL AUTHOR
C DART $Id$
C
C USAGE:    CALL BORT2 (STR1, STR2)
C   INPUT ARGUMENT LIST:
C     STR1     - CHARACTER*(*): FIRST ERROR MESSAGE TO BE PRINTED TO
C                STANDARD OUTPUT
C     STR2     - CHARACTER*(*): SECOND ERROR MESSAGE TO BE PRINTED TO
C                STANDARD OUTPUT
C
C   OUTPUT FILES:
C     UNIT 06  - STANDARD OUTPUT PRINT
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT_EXIT
C    THIS ROUTINE IS CALLED BY: ELEMDX   PARSEQ   PARUSR   PARUTG
C                               RDUSDX   SEQSDX   STRING   UFBINT
C                               UFBOVR   UFBREP   UFBSTP   VALX
C                               Normally not called by any application
C                               programs but it could be.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CHARACTER*(*) STR1, STR2
      PRINT*
      PRINT*,'**************BUFR ARCHIVE LIBRARY ABORT*****************'
      PRINT*,STR1
      PRINT*,STR2
      PRINT*,'**************BUFR ARCHIVE LIBRARY ABORT*****************'
      PRINT*
      CALL BORT_EXIT
      END

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess