text:
! Subroutine QBACK

! This software is part of the GLOW model.  Use is governed by the Open Source
! Academic Research License Agreement contained in the file glowlicense.txt.
! For more information see the file glow.txt.

! Stan Solomon, 11/1988, 11/1992, 3/2005
! New version uses updated TIE-GCM and TIME-GCM qinite.F formulation, SCS, 1/2013
! Refactored to f90, SCS, 6/2016

! Estimates background ("nighttime") ionization rates.
! Four components are used:
!   Geocoronal Lyman-beta 102.6 nm (ionizes O2 only)
!   Geocoronal He I 58.4 nm
!   Geocoronal He II 30.4 nm
!   Geocoronal Lyman-alpha 121.6 nm (ionizes NO only)

! Inputs:
!   zmaj    major species O, O2, N2 at each altitude
!   zno     nitric oxide at each altitude
!   zvcd    vertical column density for each major species above each altitude
!   photoi  array of photoionization rates for each state, species, altitude
!   f107    solar 10.7 cm radio flux activity index
!   jm      number of altitude levels
!   nmaj    number of major species (3)
!   nst     number of states
! Output:
!   Array photoi is incremented by the background rates

!Other definitions:
!   al      photon flux at 102.6 nm, 58.4 nm, 30.4 nm
!   flyan   photon flux at 121.6 nm
!   sa      absorption cross sections for O, O2, N2 at each wavelength
!   si      ionization cross sections for O, O2, N2 at each wavelength
!   salyao2 absorption cross section for O2 at 121.6 nm
!   silyano ionization cross section for NO at 121.6 nm
!   bn2p    branching ratio for N2+ from ionization of N2
!   bn1p    branching ratio for N+ from ionization of N2
!   tau     optical depth
!   qbo1    production rate of O+
!   qbo2    production rate of O2+
!   qbn2    production rate of N2+
!   qbn1    production rate of N+
!   qbno    production rate of NO+

! All units cgs.

  subroutine qback (zmaj,zno,zvcd,photoi,phono,f107,jm,nmaj,nst)

    implicit none

    integer,intent(in) :: jm, nmaj, nst
    real,intent(in) :: zmaj(nmaj,jm), zno(jm), zvcd(nmaj,jm), f107
    real,intent(inout) :: photoi(nst,nmaj,jm), phono(nst,jm)

    integer :: j, l
    real :: al(3), sa(3,3), si(3,3), salyao2, silyano, bn2p, bn1p
    real :: flyan, tau, qbo1, qbo2, qbn2, qbn1, qbno

    data al /1.5e7, 1.5e6, 1.5e6/
    data sa /      0.,  1.6e-18,       0., &
             10.2e-18, 22.0e-18, 23.1e-18, &
              8.4e-18, 16.0e-18, 11.6e-18/
    data si /      0.,  1.0e-18,       0., &
             10.2e-18, 22.0e-18, 23.1e-18, &
              8.4e-18, 16.0e-18, 11.6e-18/
    data salyao2/8.0e-21/
    data silyano/2.0e-18/
    data bn2p/0.86/
    data bn1p/0.14/

! Calculate Lyman-alpha 121.6 nm geocoronal flux as a function of F10.7:

    flyan = 5.E9*(1.+0.002*(f107-65.))

! Loop over altitudes:

    do j=1,jm
      qbo1 = 0.
      qbo2 = 0.
      qbn2 = 0.
      qbn1 = 0.

! Calculate optical depth and ionization rates for major species:

      do l=1,3
        tau=(sa(1,l)*zvcd(1,j)+sa(2,l)*zvcd(2,j)+sa(3,l)*zvcd(3,j))
        qbo1 = qbo1+al(l)*si(1,l)*zmaj(1,j)*exp(-tau)
        qbo2 = qbo2+al(l)*si(2,l)*zmaj(2,j)*exp(-tau)
        qbn2 = qbn2+bn2p*(al(l)*si(3,l)*zmaj(3,j)*exp(-tau))
        qbn1 = qbn1+bn1p*(al(l)*si(3,l)*zmaj(3,j)*exp(-tau))
      enddo

! Calculate optical depth of Ly-alpha, and ionization rate of NO:

      tau = salyao2*zvcd(2,j)
      qbno = flyan*silyano*zno(j)*exp(-tau)

! Increment ionization rate array with background ionization:

      photoi(1,1,j) = photoi(1,1,j) + qbo1
      photoi(1,2,j) = photoi(1,2,j) + qbo2
      photoi(1,3,j) = photoi(1,3,j) + qbn2
      photoi(6,3,j) = photoi(6,3,j) + qbn1
      phono(1,j) = phono(1,j) + qbno

    enddo

    return
  end subroutine qback

max_stars_repo_name: scivision/NCAR-GLOW

max_issues_repo_name: space-physics/NCAR-GLOW

max_forks_repo_name: space-physics/NCAR-GLOW


!----------------------------------------------------------------------
!                           BHMIE code
!----------------------------------------------------------------------
! This is a modified version of the BHMIE code written by Bruce Draine.
!  https://www.astro.princeton.edu/~draine/scattering.html
!
! R.T. modified following points.
!- A big static array is now defined as a dynamic array.
!- Implemented the Lentz's continued fraction method to estimate
!  the logarithmic derivative of Riccati Bessel function,
!  which is safe for very large grains.
!- Implemented adjustable real and complex precision
!- Slightly modified xstop to match the Wiscombe (1980)'s criterion.
!
!                                              Ryo Tazaki (2020.12.16)
!----------------------------------------------------------------------
module types
implicit none
integer, parameter      :: dp = selected_real_kind(P=15)
end module types

subroutine bhmie_rt(x,refrel,nang,S1,S2,Qext,Qsca,Qback,gsca)
use types
implicit none
integer::nang,j,jj,n,nn,nstop,nmx
complex(kind=dp)::refrel,y,dini
complex(kind=dp),dimension(1:2*nang-1)::S1,S2
complex(kind=dp),allocatable,dimension(:)::d,ajy
complex(kind=dp)::an,bn,an1,bn1
complex(kind=dp)::xi1,xi,dpcx
real(kind=dp)::x,qext,qsca,gsca,qback
real(kind=dp)::en,fn,p,xstop,ymod,nu
real(kind=dp)::psi0,psi1,psi,chi0,chi1,chi
real(kind=dp)::dang,pii,theta
real(kind=dp),dimension(1:nang)::amu,pi0,pi1,pi,tau
real(kind=dp) realpart
real(kind=dp) imagpart
realpart(dpcx)=(dble(dpcx))
imagpart(dpcx)=(dimag(dpcx))

! some checks...
y     = refrel*x
ymod  = abs(y)
xstop = x + 4.05_dp * x ** (1.0/3.0) + 2.0_dp
nstop = nint(xstop)

!*** require nang.ge.1 in order to calculate scattering intensities
pii  = 4.0_dp*atan(1.0_dp)
dang = 0.0_dp
if(nang.gt.1)dang=0.5_dp*pii/dble(nang-1)
do j=1,nang
        theta =real(j-1,kind=dp)*dang
        amu(j)=cos(theta)
enddo
do j=1,nang
        pi0(j)=0.0_dp
        pi1(j)=1.0_dp
enddo
nn=2*nang-1
do j=1,nn
        S1(j)=cmplx(0.0_dp,0.0_dp)
        S2(j)=cmplx(0.0_dp,0.0_dp)
enddo

!----------------------------------------------------------------------
!       logarithmic derivative
!----------------------------------------------------------------------
! 1. Standard BHMIE algorithm ( Draine's code as well )
! Warming-up terms of 15 may not be sufficient for very large spheres.
!----------------------------------------------------------------------
!nmx = nint(max(xstop,ymod)) + 15
!nmx = nint(max(xstop,ymod)) + 100
!allocate(d(1:nmx))
!if(nmx .gt. 150000) then !1.5e5 is from BTD's code.
!        write(*,*) "Error: nmx > nmxx=",150000,"for |m|x=",ymod
!        stop
!endif
!calculate D_n(mx) by downward recurrence.
!Initial value is set as D(mx)=0+0i at n=NMX.
!d(nmx)=cmplx(0.0_dp,0.0_dp)
!do n=1,nmx-1
!        en=real(nmx-n+1,kind=dp)
!        d(nmx-n) = en/y - (1.0_dp/(d(nmx-n+1) + en/y))
!enddo

!----------------------------------------------------------------------
! 2. The Lentz's continued fraction method
!    Reference: W. J. Lentz, Appl. Opt. 15. 668. (1976)
!----------------------------------------------------------------------
! Finding initial value of the logarithmic derivative "d" at
! the starting order of downward recurrence (n=nmx).
!----------------------------------------------------------------------
nmx = nint(max(xstop,ymod)) + 15
allocate(d(1:nmx),ajy(1:nmx))
! Equation (5) in Lentz (1976):
nu = real(nmx,kind=dp)+0.5_dp
ajy(1) = 2.0_dp * nu / y
do j=2,nmx
       ajy(j) = (-1.0_dp)**real(j-1,kind=dp)*2.0_dp*(nu+real(j-1,kind=dp))/y
enddo
dini = ajy(nmx)
do j=1,nmx-1
        dini = ajy(nmx-j) + 1.0_dp / dini
enddo
deallocate(ajy)
! Equation (3) in Lentz (1976):
d(nmx)=-real(nmx,kind=dp)/y+dini
! usual downward recurrence:
do n=1,nmx-1
        en=real(nmx-n+1,kind=dp)
        d(nmx-n) = (en/y) - (1.0_dp/(d(nmx-n+1) + en/y))
enddo
!----------------------------------------------------------------------

psi0 = cos(x)
psi1 = sin(x)
chi0 =-sin(x)
chi1 = cos(x)
xi1  = cmplx(psi1,-chi1)
!
qsca = 0.0_dp
gsca = 0.0_dp
p    =-1.0_dp
do n=1,nstop
        en = real(n,kind=dp)
        fn = (2.0_dp*en+1.0_dp)/(en*(en+1.0_dp))

        !Calculate psi and chi via upward recurrence:
        psi = (2.0_dp*en-1.0_dp)*psi1/x - psi0
        chi = (2.0_dp*en-1.0_dp)*chi1/x - chi0
        xi = cmplx(psi,-chi)

        !*** store previous values of an and bn for use
        !    in computation of g=<cos(theta)>
        if(n.gt.1)then
                an1=an
                bn1=bn
        endif

        !*** compute an and bn:
        an = (d(n)/refrel+en/x)*psi-psi1
        an = an/((d(n)/refrel+en/x)*xi-xi1)
        bn = (refrel*d(n)+en/x)*psi-psi1
        bn = bn/((refrel*d(n)+en/x)*xi-xi1)

        !*** augment sums for qsca and g=<cos(theta)>
        qsca=qsca+(2.0_dp*en+1.0_dp)*(abs(an)**2.0+abs(bn)**2.0)
        gsca=gsca+((2.0_dp*en+1.0_dp)/(en*(en+1.0_dp)))*&
        &        (realpart(an)*realpart(bn)+imagpart(an)*imagpart(bn))
        if(n.gt.1)then
            gsca=gsca+((en-1.0_dp)*(en+1.0_dp)/en)*&
        &      (realpart(an1)*realpart(an)+imagpart(an1)*imagpart(an)+&
        &      realpart(bn1)*realpart(bn)+imagpart(bn1)*imagpart(bn))
        endif

        !c*** now calculate scattering intensity pattern
        !    first do angles from 0 to 90
        do j=1,nang
                jj=2*nang-j
                pi(j)=pi1(j)
                tau(j)=en*amu(j)*pi(j)-(en+1.0_dp)*pi0(j)
                S1(j)=S1(j)+fn*(an*pi(j)+bn*tau(j))
                S2(j)=S2(j)+fn*(an*tau(j)+bn*pi(j))
        enddo

        !*** now do angles greater than 90 using pi and tau from
        !    angles less than 90.
        !    p=1 for n=1,3,...; p=-1 for n=2,4,...
        p=-p
        do j=1,nang-1
                jj=2*nang-j
                S1(jj)=S1(jj)+fn*p*(an*pi(j)-bn*tau(j))
                S2(jj)=S2(jj)+fn*p*(bn*pi(j)-an*tau(j))
        enddo

        !Preparation for next step
        psi0 = psi1
        psi1 = psi
        chi0 = chi1
        chi1 = chi
        xi1 = cmplx(psi1,-chi1)

        !*** compute pi_n for next value of n
        !    for each angle j, compute pi_n+1
        !    from pi = pi_n , pi0 = pi_n-1
        do j=1,nang
                pi1(j)=((2.0_dp*en+1.0_dp)*amu(j)*pi(j)-(en+1.0_dp)*pi0(j))/en
                pi0(j)=pi(j)
        enddo
enddo

!*** have summed sufficient terms.
!    now compute qsca,qext,qback,and gsca
gsca=2.0_dp*gsca/qsca
qsca=(2.0_dp/(x*x))*qsca
qext=(4.0_dp/(x*x))*realpart(S1(1))
qback=(4.0_dp*(ABS(S1(2*NANG-1))/x)**2._dp)

deallocate(d)

return
end subroutine bhmie_rt


text:
	SUBROUTINE PC_INTP  ( vlev, adata, bdata, nparms, intflg,
     +			      angflg, outdat, iret )
C************************************************************************
C* PC_INTP								*
C*									*
C* This subroutine interpolates between two levels of data.  The	*
C* data are interpolated with respect to the log of the pressure.	*
C* Pressure MUST be the first variable in the input data arrays.	*
C* If errors are encountered no output data are changed.  Therefore,	*
C* data should be set to RMISSD before calling this subroutine.		*
C*									*
C* PC_INTP  ( VLEV, ADATA, BDATA, NPARMS, INTFLG, ANGFLG,		*
C*            OUTDAT, IRET )						*
C*									*
C* Input parameters:							*
C*	VLEV		REAL		Pressure to be interpolated	*
C*	ADATA (NPARMS)	REAL		Data at first level		*
C*	BDATA (NPARMS)	REAL		Data at second level		*
C*	NPARMS		INTEGER		Size of data arrays		*
C*	INTFLG (NPARMS)	LOGICAL		Interpolation flags		*
C*	ANGFLG (NPARMS)	LOGICAL		Angle interpolation flags	*
C*									*
C* Output parameters:							*
C*	OUTDAT (NPARMS)	REAL		Data interpolated to VLEV	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-24 = vlev not between levels	*
C*					-28 = invalid pressure 		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 9/84						*
C* M. desJardins/GSFC	 9/88	GEMPAK4					*
C* G. Huffman/USRA	10/89	Check both pressures for non-positive	*
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
	REAL		adata (*), bdata (*), outdat (*)
	LOGICAL		intflg (*), angflg (*)
C*
	LOGICAL		betwen
C*
        INCLUDE         'ERMISS.FNC'
C-------------------------------------------------------------------------
C*	Check that VLEV is between two pressures input.
C
	pres1 = adata (1)
	pres2 = bdata (1)
	betwen = ( ( (pres1 .lt. pres2) .and. (pres1 .lt. vlev) .and.
     +		      (vlev .lt. pres2) )  .or.
     +		    ( (pres2 .lt. pres1) .and. (pres2 .lt. vlev) .and.
     +		      (vlev .lt. pres1) ) )
C
C*	Test for errors - can't interpolate or a non-positive pressure.
C
	IF  ( .not. betwen )  THEN
	    iret = -24
	  ELSE IF  ( ( pres1 .le. 0. ) .or. ( pres2 .le. 0. ) )  THEN
	    iret = -28
	  ELSE
	    iret = 0
C
C*	    Set up interpolation.  Check for interpolation of angles
C*	    through 360 degrees.
C
	    rmult = ALOG ( vlev / adata(1) ) /
     +		    ALOG ( bdata(1) / adata(1) )
	    outdat (1) = vlev
	    DO  i = 2, nparms
		IF  ( intflg (i) )  THEN
		    IF  ( angflg (i) )  THEN
			ang1 = AMOD ( adata(i), 360. )
			ang2 = AMOD ( bdata(i), 360. )
			IF  ( ABS (ang1 - ang2) .gt. 180.)  THEN
			    IF  ( ang1 .lt. ang2 )  THEN
				ang1 = ang1 + 360.
			      ELSE
				ang2 = ang2 + 360.
			    END IF
			END IF
			ang = ang1 + (ang2 - ang1) * rmult
			outdat (i) = amod ( ang, 360. )
		      ELSE
		        outdat (i) = adata (i) +
     +				      ( bdata(i) - adata (i) ) * rmult
		    END IF
		    IF  ( ERMISS ( adata(i) ) .or.
     +				ERMISS ( bdata(i) ) )
     +					outdat (i) = RMISSD
		END IF
	    END DO
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE add_tension(amat, wten, hx, tens, tensv, fpoly,
     1   n, nb, ioff, nmat)
      USE vspline
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER n, nb, ioff, nmat
      REAL(rprec) tens, tensv, fpoly
      REAL(rprec), DIMENSION(nmat,*) :: amat
      REAL(rprec), DIMENSION(nmat) :: wten
      REAL(rprec), DIMENSION(n) :: hx
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: k, i, koff
      REAL(rprec), DIMENSION(n) :: wten0, tenshx, work
      REAL(rprec) :: delta_x, tension
C-----------------------------------------------

!
!       THIS SUBROUTINE ADDS THE TENSION TERM WTEN*(DEL(K) - DEL(K-1))
!       TO THE REST OF THE CHI-SQ AMAT, WHERE DEL(K) = (G(K+1)-G(K))/H(K)
!       AND G(K) IS THE SECOND DERIVATIVE AT THE K-TH KNOT
!
!       IOFF ALLOWS FOR ADDING TENSION SEPARATELY TO IOTA (IOFF=0)
!       AND PRESSURE (IOFF=NIOTA) SPLINES
!
!       tens:   spline tension (optionally, at the 1st pt only;see note)
!       tensv:  vbl spline tension for n-1th point (see note)
!       fpoly:  vbl spline tension form factor (note: IF tens1<>tens2
!               THEN tension(i-th point) = tens+(tensv-tens)*(i/n-1))**fpoly)

!
!       BOUNDS CHECKING
!
      IF (n + ioff > nmat) STOP '(n+ioff>nmat)'
      IF (fpoly < 0.) STOP '(fpoly<0)'
      IF (n < 1) STOP '(n < 1)'

!
!       COMPUTE TENSION COEFFICIENTS
!

      delta_x = SUM(hx(:n-1))
      tension = 0.5*(delta_x/(n))**3
      IF (fpoly.eq.0. .or. tens.eq.tensv) THEN
         tenshx(:n-1) = tens
      ELSE
         DO i = 1, n - 1
            tenshx(i) = tens + (tensv - tens)*(REAL(i - 1,rprec)/
     1         (n - 1))**fpoly
         END DO
      ENDIF

      DO i = 1,n-1
        tenshx(i) = tension * tenshx(i) / hx(i)
        work(i) = hx(i)*(wten(i+ioff) + wten(i+ioff+1))
      ENDDO
      DO i = 2,n-1
        wten0(i) = 0.5 * ( work(i) + work(i-1) )/(hx(i) + hx(i-1))
      ENDDO
      wten0(1) = wten0(2)
      wten0(n) = wten(n+ioff)
!
!       COMPUTE, FOR K = 1,N, B(K,L)*JACOBIAN(L,I) = W(K,I),
!       WHERE JACOBIAN = D[G]/D[F] and B is TRIDIAGONAL
!       SEE EQN(27) IN PHYS.PLASMAS 1, p 2277.
!
      DO k = 1, n
         koff = k + ioff
         work(:n-1) = 0
!       SET UP COEFFICIENTS IN [G(K+1)-G(K)]/h(k) - [G(K)-G(K-1)]/h(k-1)
         IF (k .eq. 1) THEN
            work(2) = tenshx(1)*wten0(2)
            work(1) = -work(2)
         ELSE IF (k .eq. n) THEN
            work(n-1) = tenshx(n-1)*wten0(n-1)
         ELSE
            work(k-1) = tenshx(k-1)*wten0(k-1)
            work(k) = -(tenshx(k)+tenshx(k-1))*wten0(k)
            work(k+1) = tenshx(k)*wten0(k+1)
         ENDIF
         IF (nb .eq. natur) work(1) = 0
         work(n) = 0
!
!       COMPUTE work(j) = work(i)*Jacobian(i,j) and add to amat(k,j)
!
         CALL jacprod (work, hx, n, nb)
         amat(koff,1+ioff:n+ioff) = amat(koff,1+ioff:n+ioff) + work(:n)
      END DO

      END SUBROUTINE add_tension


text:
module calculate_module
    use utils_module, only: dp
    implicit none
    contains

    subroutine calculate_point(loglikelihood,prior,point,settings,nlike)
        use settings_module, only: program_settings
        implicit none
        interface
            function loglikelihood(theta,phi)
                import :: dp
                real(dp), intent(in),   dimension(:) :: theta
                real(dp), intent(out),  dimension(:) :: phi
                real(dp) :: loglikelihood
            end function
        end interface
        interface
            function prior(cube) result(theta)
                import :: dp
                real(dp), intent(in), dimension(:) :: cube
                real(dp), dimension(size(cube))    :: theta
            end function
        end interface

        type(program_settings), intent(in) :: settings
        real(dp), intent(inout) , dimension(:) :: point
        integer, intent(inout) :: nlike

        real(dp),dimension(settings%nDims)    :: cube   ! Hypercube coordinate
        real(dp),dimension(settings%nDims)    :: theta  ! Physical parameters
        real(dp),dimension(settings%nDerived) :: phi    ! derived parameters
        real(dp)                              :: logL

        cube = point(settings%h0:settings%h1)

        if ( any(cube<0) .or. any(cube>1) )  then
            theta = 0
            logL  = settings%logzero
        else
            theta = prior(cube)
            logL  = loglikelihood(theta,phi)
        end if

        if(logL>settings%logzero) nlike = nlike+1

        point(settings%p0:settings%p1) = theta
        point(settings%d0:settings%d1) = phi
        point(settings%l0) = logL

    end subroutine calculate_point

    !> Calculate a posterior point from a live/phantom point
    function calculate_posterior_point(settings,point,logweight,evidence,volume) result(posterior_point)
        use settings_module,   only: program_settings
        use utils_module,      only: logincexp
        implicit none

        type(program_settings), intent(in) :: settings
        real(dp), dimension(settings%nTotal),intent(in) :: point
        real(dp),intent(in) :: logweight
        real(dp),intent(in) :: evidence
        real(dp),intent(in) :: volume
        real(dp), dimension(settings%nposterior) :: posterior_point

        ! Volume
        posterior_point(settings%pos_X)  = volume
        ! Likelihood
        posterior_point(settings%pos_l)  = point(settings%l0)
        ! Un-normalised weighting
        posterior_point(settings%pos_w)  = logweight
        ! un-normalise cumulative weighting
        posterior_point(settings%pos_Z)  = evidence
        ! Physical parameters
        posterior_point(settings%pos_p0:settings%pos_p1) = point(settings%p0:settings%p1)
        ! Derived parameters
        posterior_point(settings%pos_d0:settings%pos_d1) = point(settings%d0:settings%d1)

    end function calculate_posterior_point

    !> This function computes the similarity matrix of an array of data.
    !!
    !! Assume that the data_array can be considered an indexed array of vectors
    !! V = ( v_i : i=1,n )
    !!
    !! The similarity matrix can be expressed very neatly as
    !! d_ij = (v_i-v_j) . (v_i-v_j)
    !!      = v_i.v_i + v_j.v_j - 2 v_i.v_j
    !!
    !! The final term can be written as a data_array^T data_array, and the first
    !! two are easy to write. We can therefore calculate this in two lines with
    !! instrisic functions
    function calculate_similarity_matrix(data_array) result(similarity_matrix)

        real(dp), intent(in), dimension(:,:) :: data_array

        real(dp), dimension(size(data_array,2),size(data_array,2)) :: similarity_matrix

        integer :: i

        similarity_matrix = spread( &
            [ ( dot_product(data_array(:,i),data_array(:,i)), i=1,size(data_array,2) ) ], &
            dim=2,ncopies=size(data_array,2) )

        similarity_matrix = similarity_matrix + transpose(similarity_matrix) - 2d0 * matmul( transpose(data_array),data_array )

    end function calculate_similarity_matrix

end module calculate_module

max_stars_repo_name: annis/cosmosis

max_issues_repo_name: annis/cosmosis

max_forks_repo_name: annis/cosmosis


text:
C> @file face.f low-level initialization drivers. Eventually to be
C> superceded by nek5000 core DG handles and operators.
c-----------------------------------------------------------------------
      subroutine iface_vert_int8cmt(nx,ny,nz,fa,va,jz0,jz1,nel)
      include 'SIZE'
      integer*8 fa(nx*nz,2*ldim,nel),va(0:nx+1,0:ny+1,jz0:jz1,nel)
      integer e,f

      n = nx*nz*2*ldim*nel
      call izero8(fa,n)

      mx1 = nx+2
      my1 = ny+2
      mz1 = nz+2
      if (ldim.eq.2) mz1=1

      nface = 2*ldim
      do e=1,nel
      do f=1,nface
         call facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,f)

         if     (f.eq.1) then ! EB notation
            ky1=ky1-1
            ky2=ky1
         elseif (f.eq.2) then
            kx1=kx1+1
            kx2=kx1
         elseif (f.eq.3) then
            ky1=ky1+1
            ky2=ky1
         elseif (f.eq.4) then
            kx1=kx1-1
            kx2=kx1
         elseif (f.eq.5) then
            kz1=kz1-1
            kz2=kz1
         elseif (f.eq.6) then
            kz1=kz1+1
            kz2=kz1
         endif

         i = 0
         do iz=kz1,kz2
         do iy=ky1,ky2
         do ix=kx1,kx2
            i = i+1
            fa(i,f,e)=va(ix,iy,iz,e)
         enddo
         enddo
         enddo
      enddo
      enddo

      return
      end

!-----------------------------------------------------------------------

      subroutine setup_cmt_gs(dg_hndl,nx,ny,nz,nel,melg,vertex,gnv,gnf)

!     Global-to-local mapping for gs

      include 'SIZE'
      include 'TOTAL'

      integer   dg_hndl
      integer   vertex(*)

      integer*8 gnv(*),gnf(*),ngv

      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

      mx = nx+2
      call set_vert(gnv,ngv,mx,nel,vertex,.false.) ! lives in navier8.f

      mz0 = 1
      mz1 = 1
      if (if3d) mz0 = 0
      if (if3d) mz1 = nz+1
      call iface_vert_int8cmt(nx,ny,nz,gnf,gnv,mz0,mz1,nelt)

      nf = nx*nz*2*ldim*nelt !total number of points on faces BETTA BE 4-byte!
      call fgslib_gs_setup(dg_hndl,gnf,nf,nekcomm,np)

      return
      end

!-----------------------------------------------------------------------

      subroutine setup_cmt_commo
      include 'SIZE'
      include 'TOTAL'
      include 'DG'

      parameter(lf=lx1*lz1*2*ldim*lelt)
      common /c_is1/ glo_num_face(lf)
     $             , glo_num_vol((lx1+2)*(ly1+2)*(lz1+2)*lelt)
      integer*8 glo_num_face,glo_num_vol,ngv

      call setup_cmt_gs(dg_hndl,lx1,ly1,lz1,nelt,nelgt,vertex,
     >                  glo_num_vol,glo_num_face)
      call cmt_set_fc_ptr(nelt,lx1,ly1,lz1,ndg_face,iface_flux)

      return
      end

!-----------------------------------------------------------------------

      subroutine cmt_set_fc_ptr(nel,nx,ny,nz,nface,iface)

!     Set up pointer to restrict u to faces ! NOTE: compact
! JH062314 Now 2D so we can strip faces by element and not necessarily
!          from the whole field

      include 'SIZE'
      include 'TOTAL'

      integer nx, ny, nz, nel
      integer nface,iface(nx*nz*2*ldim,*)
      integer e,f,ef

      call dsset(nx,ny,nz) ! set skpdat. lives in connect1.f

      nxyz = nx*ny*nz
      nxz  = nx*nz
      nfpe = 2*ldim
      nxzf = nx*nz*nfpe ! red'd mod to area, unx, etc.

      do e=1,nel
      do f=1,nfpe

         ef     = eface(f)
         js1    = skpdat(1,f)
         jf1    = skpdat(2,f)
         jskip1 = skpdat(3,f)
         js2    = skpdat(4,f)
         jf2    = skpdat(5,f)
         jskip2 = skpdat(6,f)

         i = 0
         do j2=js2,jf2,jskip2
         do j1=js1,jf1,jskip1

            i = i+1
            k = i+nxz*(ef-1)           ! face numbering
            iface(k,e) = j1+nx*(j2-1)  ! cell numbering

         enddo
         enddo

      enddo
      enddo
      nface = nxzf*nel

      return
      end

!-----------------------------------------------------------------------

      subroutine full2face_cmt(nel,nx,ny,nz,iface,faces,vols)

! JH062314 Store face data from nel full elements (volume data). Merely
!          selection for the time being (GLL grid), but if we need to
!          extrapolate to faces (say, from Gauss points), this is where
!          we'd do it.

      include 'SIZE'
      include 'TOTAL'

      integer  nel,nx,ny,nz
      integer  iface(nx*nz*2*ldim,1)
      real     faces(nx*nz   ,2*ldim,nel)
      real     vols (nx,ny,nz       ,1  )
      integer  e,i,j

      n= nx*nz*2*ldim
      do e=1,nel
      do j=1,n
         i=iface(j,e)
         faces(j,1,e) = vols(i,1,1,e)
      enddo
      enddo

      return
      end

!-----------------------------------------------------------------------

      subroutine add_face2full_cmt(nel,nx,ny,nz,iface,vols,faces)

      include 'SIZE'
      include 'TOTAL'

      integer   nel,nx,ny,nz
      integer   iface(nx*nz*2*ldim,*)
      real     faces(nx*nz   ,2*ldim,*  )
      real     vols (nx,ny,nz       ,nel)
      integer  ie,i,j

      n= nx*nz*2*ldim
      do ie=1,nel
      do j=1,n
         i=iface(j,ie)
         vols(i,1,1,ie) = vols(i,1,1,ie)+faces(j,1,ie)
      enddo
      enddo

      return
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
This is the most splendid little Gift Shops gift store! If you love handmade, oneofakind clothing, you will love this place. They have cotton and wool clothing and accessories (even handmade felted slippers!) that are colorful, authentic, affordable, and made in Nepal. Prices for tops are $35$45 (These are really high quality! Cotton, handsewn, richlydyed, embroidered, each one truly unique). The owners are really nice (family owned and operated) and the atmosphere is peaceful, uncluttered, and welcoming. A taste of beautifully creative culture.

20110708 17:25:41 nbsp I bought a present here for my dad. The prices were very, very reasonable considering that the items are handmade. Beautiful, colorful, and unique gifties! ) Users/EliseKane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Tristan Harward

About: A UC Berkeley UC Undergrad MisMarketing undergrad in Computer Science, with several friends in Davis, mainly Users/CindySperry, who counts as several friends. Once considered going to UC Davis because of all the trees and other green things. (Like Bike Signal bike signals and Ducks ducks and stagnant Putah Creek.)

20051118 22:40:32 nbsp I hereby compliment him on his lovely green background. Users/CindySperry


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!     Module DELWAQ2:
!     - Encapsulate the interface of DELWQ2 and DLWQI0:
!       A, J and C are now pointers to arrays
!
      MODULE DELWAQ2

      CONTAINS

      SUBROUTINE DELWQ2 ( A, J, C, IMAXA, IMAXI, IMAXC, INIT,
     &                    ACTION, DLWQD                     )
!
!     Deltares     SECTOR WATERRESOURCES AND ENVIRONMENT
!
!            DELWAQ - Deltares WAter Quality programme
!
!                     Version 4.30           june 1997
!                     Replaces:
!                     Version 4.22           feb 1997
!                     Version 4.21           jan 1997
!                     Version 4.20           may 1996
!                     Version 4.01           may 1994
!                     Version 4.00           september 1993
!                     Release 3.05           november 1991.
!                     Release 3.0 - RWS-DIV, june 1988.
!                     Release 2.0 of july     1984.
!                     Release 1.0 of december 1981.
!
!     INFORMATION   : Deltares
!                     L. Postma,
!                     Rotterdamse weg 185,
!                     P.O. Box 177,
!                     2600 MH Delft,
!                     Netherlands.
!                     telephone (31) 15-569353
!                     telefax   (31) 15-619674
!
!     FUNCTION      : Performs the waterquality simulations on a
!                     consistent set of binary intermediate files.
!
!     LOGICAL UNITS : LUNIN  , input , binary common-block file
!                       *    , output, user console file
!
!     SUBROUTINES CALLED : DLWQI0, initialises the system
!                          DLWQN1, first   integration procedure
!                          DLWQN2, second  integration procedure
!                          DLWQN3, third   integration procedure
!                          DLWQN4, fourth  integration procedure
!                          DLWQN5, fifth   integration procedure
!                          DLWQN6, sixth   integration procedure
!                          DLWQN7, seventh integration procedure
!                          DLWQN8, eighth  integration procedure
!                          DLWQN9, nineth  integration procedure
!                          DLWQNB, tenth   integration procedure
!                          DLWQNC, 11th    integration procedure
!                          DLWQND, 12th    integration procedure
!                          DLWQNE, 13th    integration procedure
!                          DLWQNA, 14th    integration procedure
!                          DLWQNF, 15th    integration procedure
!                          DLWQNG, 16th    integration procedure
!                          DLWQNH, 17th    integration procedure
!                          DLWQNI, 18th    integration procedure
!                          DLWQNJ, 19+20   integration procedure
!                          SRSTOP, stops execution
!                          DHOPNF, opens files
!
!     PARAMETERS    :
!
!     NAME    KIND     LENGTH  DESCRIPTION
!     ---------------------------------------------------------
!     A       REAL     IMAXA   real      workspace array
!     J       INTEGER  IMAXJ   integer   workspace array
!     C       CHAR*20  IMAXC   character workspace array
!     DLWQD   TYPE(..) 1       derived type for persistent storage
!     IMAXA   INTEGER  1       maximum real      workspace array
!     IMAXI   INTEGER  1       maximum integer   workspace array
!     IMAXC   INTEGER  1       maximum character workspace array
!     INIT    LOGICAL  1       if T boot the system if F no initialisation
!     ACTION  INTEGER  1       indication of the action to be performed
!
      USE DIO_PLT_RW
      use grids
      USE DLWQI0_MOD
      USE Timers
      use m_timers_waq
      use m_couplib
      use delwaq2_data

      implicit none

      include 'actions.inc'
!
!     Declaration of arguments
!
      INTEGER       IMAXA , IMAXI , IMAXC
      INTEGER, DIMENSION(:), POINTER          :: J
      REAL, DIMENSION(:), POINTER             :: A
      CHARACTER(LEN=*), DIMENSION(:), POINTER :: C
      LOGICAL                                 :: INIT
      INTEGER                                 :: ACTION
      TYPE(DELWAQ_DATA), TARGET               :: DLWQD
      type(GridPointerColl), pointer          :: GridPs               ! collection of all grid definitions

!
!     COMMON  /  SYSN   /   System characteristics
!
      INCLUDE 'sysn.inc'
!
!     COMMON  /  SYSI  /    Timer characteristics
!
      INCLUDE 'sysi.inc'
!
!     COMMON  /  SYSA   /   Pointers in real array workspace
!
      INCLUDE 'sysa.inc'
!
!     COMMON  /  SYSJ   /   Pointers in integer array workspace
!
      INCLUDE 'sysj.inc'
!
!     COMMON  /  SYSC   /   Pointers in character array workspace
!
      INCLUDE 'sysc.inc'
!
!     PARAMETERS    :
!
!     NAME    KIND     LENGTH  DESCRIPTION
!     ---------------------------------------------------------
!     LUNIN   INTEGER  1       unit nummer of the common BOOT-file
!     IPAGE   INTEGER  1       pagelength for output in lines
!     NLUN    INTEGER  1       number of unit numbers
!     LCHMAX  INTEGER  1       length file names
!
!
!     Local declarations
!
      INTEGER, PARAMETER ::   LUNIN  =   914
      INTEGER, PARAMETER ::   IPAGE  =    64
      INTEGER, PARAMETER ::   NLUN   =    50
      INTEGER, PARAMETER ::   LCHMAX =   255
!
!           input structure for boot-file
!
      INTEGER, DIMENSION(INSIZE)  :: IN
      INTEGER, DIMENSION(IISIZE)  :: II
      EQUIVALENCE ( IN( 1) , NOSEG ) , ( II( 1) , ITSTRT )
!
      INTEGER, SAVE            :: LUN(NLUN)
      CHARACTER*(LCHMAX), SAVE :: LCHAR(NLUN)
      integer, save            :: filtype(nlun)
      CHARACTER*(LCHMAX), SAVE :: RUNID
      LOGICAL, SAVE            :: INIT2        = .TRUE. ! To suppress the start-up screen
      LOGICAL, SAVE            :: INIT_COUPLIB = .TRUE. ! Initialise the couplib library only once

!     Common to define external communications in SOBEK
!     OLCFWQ             Flag indicating ONLINE running of CF and WQ
!     SRWACT             Flag indicating active data exchange with SRW
!     RTCACT             Flag indicating output for RTC

      character*(lchmax)       :: inifil, dioconfig
      logical                  :: lfound
      integer                  :: idummy, ierr2
      real                     :: rdummy
      CHARACTER                :: cdummy
      CHARACTER*2              :: C2
      LOGICAL                  :: OLCFWQ, SRWACT, RTCACT, DDWAQ
      COMMON /COMMUN/             OLCFWQ, SRWACT, RTCACT, DDWAQ
!
      integer(4), save         :: ithndl = 0
!
!     Local variables
!
      INTEGER, SAVE            :: INDX
      INTEGER                  :: IERR
      INTEGER                  :: IMR
      INTEGER                  :: IMI
      INTEGER                  :: IMC
      INTEGER                  :: ILUN
      INTEGER                  :: IERRD
      INTEGER                  :: K
      LOGICAL                  :: NOLIC
!
      IF ( INIT ) THEN
         call timini ( )
         ! for openda-usage, where multiple instances are launched,
         ! the time module does not work correctly.
         if ( dlwqd%set_timer ) timon = .true.
         timon = .true.
         if (timon) call timstrt( "delwaq2", ithndl )
!
!        boot the system; read dimensions of sysn from delwaq03.wrk-file
!
         CALL DHGNAM(RUNID,'.mon')
         LCHAR(1) = TRIM(RUNID)//'-delwaq03.wrk'
         CALL DHOPNF ( LUNIN , LCHAR(1), 1     , 2     , IERR  )
         IF ( IERR .GT. 0 ) GOTO 999
         READ  ( LUNIN )   IN
         READ  ( LUNIN )   II
         READ  ( LUNIN )   IMR     , IMI  , IMC
         READ  ( LUNIN ) ( LUN    (K), K = 1, NOLUN )
         READ  ( LUNIN ) ( LCHAR  (K), K = 1, NOLUN )
         READ  ( LUNIN ) ( filtype(K), K = 1, NOLUN )
         DO 5 ILUN = 1, NOLUN
            CLOSE ( LUN(ILUN) )
    5    CONTINUE
         close(lunin)
!
!        start couplib
!
         if ( init_couplib ) then
             init_couplib = .false.
             call couplib_init(lunout=lun(19), idebug=0)
         endif
!
!        store number of processes and own process number in sysn
!
         mypart = myprc
         npartp = numprc
!
         if (mypart .gt. 1) then
! open lokale monitor-files voor uitvoer tijdens initialisatie
            write(lchar(19),'(a,i3.3,a)') 'part-', mypart,'.mon'
         end if

         CALL DHOPNF ( LUN(19) , LCHAR(19) , 19    , 1    , IERRD  )
         CALL SETMLU ( LUN(19) )

         IF (MYPART .EQ. 1) THEN

!      Initialise communication options SOBEK

            OLCFWQ = .FALSE.
            SRWACT = .FALSE.
            RTCACT = .FALSE.
            DDWAQ  = .FALSE.
            LCHAR(44) = ' '
            LUN(44)   = LUN(43) + 1

            call getcom ( '-i'  , 3    , lfound, idummy, rdummy,
     +                   inifil, ierr2)
            if ( lfound ) then
               if ( ierr2.ne. 0 ) then
                  inifil = ' '
               endif
            else
               inifil = 'delwaq.ini'
            endif
            open(lunin,file=inifil,status='old',err=123)
            write(lun(19),*) ' Using options from ini file : ',trim(inifil)
            call gkwini(lunin,'SimulationOptions','OnLineWQ',c2)
            if ( c2 .eq. '-1' ) then
               olcfwq = .true.
               write(lun(19),*) ' online coupling with FLOW module activated'
            endif
            call gkwini(lunin,'SimulationOptions','OutputRTC',c2)
            if ( c2 .eq. '-1' ) then
               rtcact = .true.
               write(lun(19),*) ' RTC coupling activated'
            endif
            call gkwini(lunin,'SimulationOptions','SRW',c2)
            if ( c2 .eq. '-1' ) then
               srwact = .true.
               write(lun(19),*) ' SRW coupling activated'
            endif
            call gkwini(lunin,'General','DIOConfigFile',dioconfig)
            call gkwini(lunin,'General','ProgressFile',lchar(44))
            close (lunin)
 123        continue

         ! check ddwaq from commandline

            call getcom ( '-d'  , 0    , lfound, idummy, rdummy,
     +           cdummy, ierr2)
            if ( lfound ) then
               DDWAQ = .TRUE.
               write(lun(19),*) ' DDWAQ coupling activated'
            endif

       ! initialise DIO

            IF ( OLCFWQ .OR. SRWACT .OR. RTCACT .OR. DDWAQ ) THEN
               if ( dioconfig .ne. ' ' ) then
                  write(lun(19),*) ' Using DelftIO ini file : ',trim(dioconfig)
                  CALL DIOINIT(dioconfig)
               else
                  write(lun(19),*) ' Using default DelftIO ini file'
                  CALL DIOINIT()
               endif
            ENDIF
!
!     The unlocking !
!
            IF ( INIT2 ) THEN
               INIT2 = .FALSE.
               IF ( NOQ3 .GT. 0 ) THEN
                  CALL UNLOCK (LUN(19),.TRUE.,NOLIC)
               ELSE
                  CALL UNLOCK (LUN(19),.FALSE.,NOLIC)
               ENDIF
            ENDIF

            IF (ACTION .EQ. ACTION_FULLCOMPUTATION) THEN
               WRITE(*,*)
               WRITE(*,'(A9,A)') '  runid: ',TRIM(RUNID)
               WRITE(*,*)
            ENDIF

            if ( nolic .and. noseg > 150 ) then
               write(*,'(//a)') 'Error: Authorisation problem'
               write(*,'(a)')   '       No valid license, so the number
     & of segments is limited to 150'
               call srstop(1)
            endif
         endif
!
!        end of reading master proces
!
!        initialize timers
!
         call timers_waq_init(lun(19))
         call couplib_timers_init(timer_start_couplib,
     +       couplib_max_timers, measr_idle=.false.)
!
         call timer_start(timer_total)
!
         call timer_start(timer_init)

! collaborative call to i0
         call sync_processes()
!
         IERR = 0
         gridps => dlwqd%gridps
         call dlwqi0 ( nlun   , a      , j      , c      , imaxa  ,
     &                 imaxi  , imaxc  , ipage  , lun    , lchar  ,
     &                 filtype, gridps , dlwqd  , ierr   )
!
         if (mypart .eq. 1) then
            CLOSE ( LUNIN )
            IF ( IERR .GT. 0 ) GOTO 992
!
!        end of initialisation
!
            WRITE ( * , *)
            WRITE ( * , * ) ' SIMULATION STARTED '
            WRITE ( * , * )
            WRITE ( * , * ) ' INTEGRATION ROUTINE =', intsrt
         endif
         call timer_stop(timer_init)
      ENDIF
!     SOBEK external communications ONLY implemented in scheme 10!
      IF ( OLCFWQ .OR. SRWACT .OR. RTCACT ) THEN
          IF ( INTSRT .NE. 10 .AND. (INTSRT .NE. 15 .AND. INTSRT .NE. 19) ) GOTO 991
      ENDIF

!
!     Store the local persistent variables
!
      DLWQD%II = II
      DLWQD%IN = IN

!
!     Check restrictions of Parallel computing
!
      IF ( NPARTp .GT. 1 ) THEN
!
!        Parallel computing only allowed for schemes 1, 5, 12
!
         IF ( INTSRT.NE.1 .AND. INTSRT.NE.5 .AND. INTSRT.NE.12 ) GOTO 993
!
!        Bottom layers not allowed in parallel runs
!
         IF ( NOQ4.GT.0 ) GOTO 994
!
!        Multiple grids not allowed in parallel runs
!
         IF ( NOGRID.GT.1 ) GOTO 995
!
!        Use of DelftIO coupling not allowed in parallel runs
!
         IF ( OLCFWQ .OR. SRWACT .OR. RTCACT .OR. DDWAQ ) GOTO 996
      ENDIF

!         branch to the appropriate integration option

      select case ( intsrt )

         case (  0 )     !      not transport, just processes
            call timer_start(timer_offs_intsrt+0)
            call dlwqn0 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+0)

         case (  1 )     !      backward in space and time
            call timer_start(timer_offs_intsrt+1)
            call dlwqn1 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+1)

         case (  2 )     !      modified 2nd order Runge Kutta
            call timer_start(timer_offs_intsrt+2)
            call dlwqn2 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+2)

         case (  3 )     !      2nd order Lax Wendroff
            call timer_start(timer_offs_intsrt+3)
            call dlwqn3 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+3)

         case (  4 )     !      Aternating direction implicit
            call timer_start(timer_offs_intsrt+4)
            call dlwqn4 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+4)

         case (  5 )     !      Flux corrected transport
            call timer_start(timer_offs_intsrt+5)
            call dlwqn5 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+5)

         case (  6 )     !      Direct steady state, backward differences in space
            call timer_start(timer_offs_intsrt+6)
            call dlwqn6 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+6)

         case (  7 )     !      Direct steady state, central differences in space
            call timer_start(timer_offs_intsrt+7)
            call dlwqn7 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+7)

         case (  8 )     !      Iteratively steady state, backward differences in space
            call timer_start(timer_offs_intsrt+8)
            call dlwqn8 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+8)

         case (  9 )     !      Iteratively steady state, central differences in space
            call timer_start(timer_offs_intsrt+9)
            call dlwqn9 ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+9)

         case ( 10 )     !      Fully implicit, direct method, upwind
            call timer_start(timer_offs_intsrt+10)
            call dlwqnb ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+10)

         case ( 11 )     !      Horizontal explicit upwind, vertical implicit central
            call timer_start(timer_offs_intsrt+11)
            call dlwqnc ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+11)

         case ( 12 )     !      Horizontal explicit FCT   , vertical implicit central
            call timer_start(timer_offs_intsrt+12)
            call dlwqnd ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+12)

         case ( 13 )     !      Horizontal explicit upwind, vertical implicit upwind
            call timer_start(timer_offs_intsrt+13)
            call dlwqne ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+13)

         case ( 14 )     !      Horizontal explicit FCT   , vertical implicit upwind
            call timer_start(timer_offs_intsrt+14)
            call dlwqna ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+14)

         case ( 15 )     !      GMRES, horizontal upwind, vertical upwind
            call timer_start(timer_offs_intsrt+15)
            call dlwqnf ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+15)

         case ( 16 )     !      GMRES, horizontal upwind, vertical central
            call timer_start(timer_offs_intsrt+16)
            call dlwqng ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+16)

         case ( 17 )     !      stationary GMRES, horizontal upwind, vertical upwind
            call timer_start(timer_offs_intsrt+17)
            call dlwqnh ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+17)

         case ( 18 )     !      stationary GMRES, horizontal upwind, vertical central
            call timer_start(timer_offs_intsrt+18)
            call dlwqni ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+18)

         case ( 19 )     !      TRISULA-ADI 1 (vertically upwind)
            call timer_start(timer_offs_intsrt+19)
            call dlwqnj ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+19)

         case ( 20 )     !      TRISULA-ADI 2 (vertically central)
            call timer_start(timer_offs_intsrt+20)
            call dlwqnj ( a , j , c , lun , lchar, action, dlwqd, gridps )
            call timer_stop(timer_offs_intsrt+20)

         case ( 21 )     !      Self adjusting teta method (limiter Salezac)
            call dlwqnm ( a , j , c , lun , lchar, action, dlwqd, gridps )

         case ( 22 )     !      Self adjusting teta method (limiter Boris and Book)
            call dlwqnm ( a , j , c , lun , lchar, action, dlwqd, gridps )

         case ( 23 )     !      Leonards QUICKEST
            call dlwqno ( a , j , c , lun , lchar, action, dlwqd, gridps )

         case ( 24 )     !      Local flexible time step method by Leonard Postma
            call dlwqnp ( a , j , c , lun , lchar, action, dlwqd, gridps )

         case default
            goto 990

      end select

      IF ( ACTION == ACTION_FINALISATION    .OR.
     &     ACTION == ACTION_FULLCOMPUTATION        ) THEN

!     print timer-results
!     Note: removed printing of timers to monitoring file

          call timer_stop(timer_total)
!         call timers_waq_print()

          call sync_processes()
          call couplib_stop()

          if ( timon ) then
             call timstop ( ithndl )
             call timdump ( TRIM(RUNID)//'-timers.out' )
             call timfinalize()
          endif

      endif

      return

  990 WRITE ( * , * ) ' ERROR: INTEGRATION OPTION NOT IMPLEMENTED'
      CALL SRSTOP(1)
  991 WRITE ( * , * ) ' ERROR: INTEGRATION OPTION NOT IMPLEMENTED in online mode'
      CALL SRSTOP(1)
  992 WRITE ( * , * ) ' ERROR : INITIALISATION FAILED'
      CALL SRSTOP(1)
  993 WRITE ( * , '(/,1x,a,i4,/)' )
     + ' INTEGRATION OPTION NOT IMPLEMENTED in parallel mode, npart=',npartp
      CALL SRSTOP(1)
  994 WRITE ( * , '(/,1x,a,i4,/)' )
     + ' BOTTOM LAYERS NOT SUPPORTED in parallel mode, npart=',npartp
      CALL SRSTOP(1)
  995 WRITE ( * , '(/,1x,a,i4,/)' )
     + ' MULTIPLE GRIDS NOT SUPPORTED in parallel mode, npart=',npartp
      CALL SRSTOP(1)
  996 WRITE ( * , '(/,1x,a,i4,/)' )
     + ' ON-LINE MODE NOT SUPPORTED in parallel mode, npart=',npartp
      CALL SRSTOP(1)
  999 WRITE ( * , * ) ' ERROR: NO VALID SET OF MODEL-INTERMEDIATE-FILES'
      CALL SRSTOP(1)
      END SUBROUTINE DELWQ2

      END MODULE DELWAQ2

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
The Tiki Hut is a combination Tiki and Pirate store that offers a selection of gift items.

Founded by a 87 graduate of UC Davis, they are (in their words), on a mission from Tiki... to sail our Pirate Ship into any tropical oasis possible. Escape from the boring, hot, typical day in the Valley and begin a quest into the exotic! We carry goods from around the world and loot from the finest cargo holds ever plundered. Pearls, silver and gold, finely crafted tiki mugs for boutique Rums and other fine elixirs grace the stock of the Tiki Hut in Davis. Music, art, smells, and tastes abound in the most fun new store this side of Tortuga/Yolo bypass.


text:
! ===============================================================================
Subroutine torsionalSpring(AE, BE, MaxLRows, XLL, NDim, iDofT, NodElt, Sol0, Sol1, CommonPar, &
									Param, JParam, DelT, DTm, Time)
! ===============================================================================

	use funcAux , only : getSlice, setBlockToVector, setBlockDiagonalConstantToMatrix
	use globalVariables, only : getF, NdimE
	use fibresLib

	implicit none

	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT, NodElt! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	Real*8 :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), Param(*), JParam(*), CommonPar(*)
!~ 	!   =====   END ARGUMENTS  =======

	! ======== COMMONPAR ARGUMENTS =====================
	integer iShiftUf !> @param  = nint(CommonPar(1))
	integer iShiftDeltaUf !> @var  = nint(CommonPar(2))
	integer iFType !> @var  = nint(commonPar(4))
	real*8 :: eta
	integer :: flagIncrementInside

	!====== OTHER VARIABLES ==============================
	integer , parameter :: NodG = 3
	Real*8 :: SolU(NdimE*NodG), SolUf(NdimE*NodG), Xel(NdimE*NodG)
 	Real*8 :: Kstiff(NdimE*NodG,NdimE*NodG), Res(NdimE*NodG), F(NdimE,NdimE)
 	integer :: i, ip

	iShiftUf = nint(CommonPar(1))
	iShiftDeltaUf = nint(CommonPar(2)) !!! now it's not being used
	iFType = nint(commonPar(3))
	eta = commonPar(4)
	flagIncrementInside = nint(commonPar(5))

	call getSlice(Xel,XLL,1,NodG,1 ,NdimE, Ndim)
	call getSlice(SolUf,Sol1,1,NodG,iShiftUf + 1 ,iShiftUf + NdimE, iDofT)

	call getF(F,iFtype) !!! with identity summed

	F = F - deltaKron2(1:NdimE,1:NdimE)

	do i = 1, NodG
		ip = (i-1)*NdimE
		SolU(ip+1:ip+NdimE) = matmul(F , Xel(ip+1:ip+NdimE) )
	end do

	SolU = SolU + SolUf

	call calc_Kstiff_and_Res_TorsionalSpring(Kstiff, Res, Xel, SolU, eta)

	if(iShiftUf == iShiftDeltaUf .and. flagIncrementInside == -1) then ! solver python updates outside
		Res = Res + matmul(Kstiff,SolUf)
	end if

	call setBlockToVector(BE,Res,NodG,iDofT, NdimE, iShiftDeltaUf ) !(A,B, NodG, iDofA, iDofB, iShiftA) , B in A
	call setBlockToMatrixSquare(AE,Kstiff,NodG,iDofT, NdimE, iShiftDeltaUf) !(A,B, NodG, iDofA, iDofB, iShiftA) , B in A

 	if(iShiftUf /= iShiftDeltaUf .and. flagIncrementInside > -1) then
		call setBlockDiagonalConstantToMatrix(AE,1.0d0,NodG, iDofT, NdimE, iShiftUf, iShiftUf)
		call setBlockDiagonalConstantToMatrix(AE,-1.0d0,NodG, iDofT, NdimE, iShiftUf, iShiftDeltaUf)
		call setBlockToVector(BE,SolUf,NodG,iDofT, NdimE, iShiftUf ) !(A,B, NodG, iDofA, iDofB, iShiftA) , B in A
	end if

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine torsionalSpringS(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    use funcAux, only : setCoupling_pureDisplacement, setSimbolicBlockDiagonalToMatrix, numprint

    IMPLICIT NONE
    Integer :: MaxLRows,Ndim,iDofT, iAdd
    Integer   Coupling(MaxLRows,MaxLRows)
    Real*8 CommonPar(*)

    integer , parameter :: NdimE = 2, NodG = 3
    integer :: iShiftUf, iShiftDeltaUf, flagIncrementInside

	iShiftUf = nint(CommonPar(1))
	iShiftDeltaUf = nint(CommonPar(2)) !! not being used at the moment
	flagIncrementInside = nint(commonPar(5))

	call setCoupling_pureDisplacement(Coupling,NodG,idofT,iShiftDeltaUf,NdimE)

	if(iShiftUf /= iShiftDeltaUf .and. flagIncrementInside > -1) then
		call setSimbolicBlockDiagonalToMatrix(Coupling,NodG, iDofT, NdimE, iShiftUf, iShiftUf)
		call setSimbolicBlockDiagonalToMatrix(Coupling,NodG, iDofT, NdimE, iShiftUf, iShiftDeltaUf)
	end if

end Subroutine

max_stars_repo_name: felipefr/gpmaterials

max_issues_repo_name: felipefr/gpmaterials

max_forks_repo_name: felipefr/gpmaterials


program poisson_mpi1D
	use poisson_module
	use mpi
	implicit none

	integer :: i,j,k,ret
	real, allocatable, dimension(:) :: u,u0,f,u_global
	integer, allocatable, dimension(:) :: buf
	real :: h, l2
	integer :: prank, psize, ierror, tag, status(MPI_STATUS_SIZE)
	double precision :: t_start, t_end
	integer :: n_loc, offset, n_min;
	integer :: j_start, j_end, i_start, i_end
	real :: x,y

! Initialization of the MPI library
	call MPI_INIT(ierror)
	call MPI_COMM_RANK(MPI_COMM_WORLD, prank, ierror)
	call MPI_COMM_SIZE(MPI_COMM_WORLD, psize, ierror)

	n_min = (N+1)
	h = 1. / n_min
! Divide n in psize and distribute the excess to the (n % psize) proc
	if (prank < mod(n_min,psize)) then
		n_loc = (n_min / psize) + 1
	else
		n_loc = (n_min / psize)
	endif
! Computing the offset of where in the global array the local array is
! located. (this is needed to initialize f) It could be computed locally
! without communication but this way you see a Allgather
	allocate(buf(psize))
	call MPI_ALLGATHER(n_loc, 1, MPI_INTEGER, buf, 1, MPI_INTEGER, MPI_COMM_WORLD,ierror)
	offset = 0
	do i = 1, prank
		offset = offset + buf(i)
	enddo
! add 2 for east and west ghost
	n_loc = n_loc + 2

	allocate(u(n_min*n_loc),u0(n_min*n_loc),f(n_min*n_loc))
	if (prank.eq.0) allocate(u_global(n_min*n_min))

	u = 0.
	u0 = 0.
	f = 0.

	do i = 2, n_min-1
		do j = 2, n_loc-1
!			x = real(i) * h
!			y = real(j - j_start + offset) * h
!			f(n_loc*i+j) = -2.*100. * M_PI * M_PI * sin(10.*M_PI*x) * sin(10.*M_PI*y)
			f(n_loc*i+j) = -2.*100. * M_PI * M_PI * sin(10.*M_PI*(i-1)*h) * sin(10.*M_PI*(j-1)*h)

		enddo
	enddo

!if (prank.eq.1) write(*,*) prank, f

	k = 0
	l2 = 1.
	t_start = MPI_WTIME()

!	do while (l2 > eps)
	do while (k < 400)
		l2 = 0.
		j_start = 2
		j_end   = n_loc-1
		i_start = 2
		i_end   = n_min-1

! send recv with left proc
		if (prank > 0)  then
			call MPI_SENDRECV(u0(n_min+2), n_min, MPI_REAL, prank - 1, NORTH, u0(1), n_min, MPI_REAL, prank - 1, SOUTH, MPI_COMM_WORLD, status, ierror)
!			call MPI_SENDRECV(u0(n_min:2*n_min), n_min, MPI_REAL, prank - 1, NORTH, u0(1:n_min), n_min, MPI_REAL, prank - 1, SOUTH, MPI_COMM_WORLD, status, ierror)
!				call MPI_SENDRECV(u0(:,2), n_min, MPI_REAL, prank - 1, NORTH, u0(:,1), n_min, MPI_REAL, prank - 1, SOUTH, MPI_COMM_WORLD, status, ierror)
		endif

! send recv with right proc
		if(prank < psize - 1) then
			call MPI_SENDRECV(u0(n_loc-1), n_min, MPI_REAL, prank + 1, SOUTH, u0(n_loc), n_min, MPI_REAL, prank + 1, NORTH, MPI_COMM_WORLD, status, ierror)
!			call MPI_SENDRECV(u0((n_loc-1)*n_min:(n_loc-1)*2*n_min), n_min, MPI_REAL, prank + 1, SOUTH, u0(n_loc:n_loc+n_min), n_min, MPI_REAL, prank + 1, NORTH, MPI_COMM_WORLD, status, ierror)
!			call MPI_SENDRECV(u0(:,n_loc-1), n_min, MPI_REAL, prank + 1, SOUTH, u0(:,n_loc), n_min, MPI_REAL, prank + 1, NORTH, MPI_COMM_WORLD, status, ierror)
		endif

		do i = i_start, i_end
			do j = j_start,j_end
				u(n_loc*i+j) = 0.25 * ( u0(n_loc*(i-1)+j) + u0(n_loc*(i+1)+j) + u0(n_loc*i+(j-1)) + u0(n_loc*i+(j+1)) - f(n_loc*i+j)*h*h);
				! L2 norm
				l2 = l2 + (u0(n_loc*i+j) - u(n_loc*i+j))*(u0(n_loc*i+j) - u(n_loc*i+j));
			enddo
		enddo

! reduce the 12 norm to every proc
		call MPI_ALLREDUCE(MPI_IN_PLACE, l2, 1, MPI_REAL, MPI_SUM, MPI_COMM_WORLD, ierror);

		call MPI_Gather(u(1),(n_min*(n_loc-2)),MPI_REAL, u_global,(n_min*(n_loc-2)) ,  MPI_REAL, 0, MPI_COMM_WORLD,ierror)
!		if (prank.eq.0) call write_to_file_binary1D(n_min, u_global, k, -1., 1.)

		! output
		if (prank .eq. 0) then
			write(*,*) 'iteration',k,l2
		endif
		k = k + 1
		u0 = u
	enddo
	t_end = MPI_WTIME()
	if (prank .eq. 0) then
		write(*,*) 'T',(t_end - t_start),'[s]'
		write(*,*) 'Nb steps',k,'l2',l2
		write(*,*) 'Size',n_min
	endif

	deallocate(u,u0,f)
	if (prank.eq.0) deallocate(u_global)

	call MPI_FINALIZE(ierror)

end program poisson_mpi1D


MODULE Subrotinas_MG_GS

   CONTAINS

   SUBROUTINE Atualiza_Variaveis_MG_GS_P( Lx_o, Ly_o, &
                                          Nx_o, Ny_o, &
                                          tol_o, &
                                          ite_down_o, ite_up_o, &
                                          max_ciclos_V_o, ht, Re, Tn_o )

      USE Variaveis_Solvers_P

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: Nx_o, Ny_o, ite_down_o, ite_up_o, max_ciclos_V_o, Tn_o
      DOUBLE PRECISION, INTENT(IN) :: Lx_o, Ly_o, tol_o, ht, Re

      INTEGER :: level
      DOUBLE PRECISION :: hx, hy

      Lx = Lx_o
      Ly = Ly_o
      Nx = Nx_o
      Ny = Ny_o
      tol = tol_o
      ite_down = ite_down_o
      ite_up   = ite_up_o
      max_ciclos_V = max_ciclos_V_o
      Tn = Tn_o

      ! Calculando o numero de niveis
      Levels = LOG( DBLE(Nx))/LOG( DBLE(2) )

      ! lv_if => duas dimensoes 1: inicio da malha, 2: fim da malha
      ALLOCATE( lv_if(Levels, 2), lv_nx(Levels), lv_ny(Levels), normas_L2(max_ciclos_V, Tn) )

      ! Inicio da primeira malha comeca em 1
      lv_if(1,1) = 1

      ! Fim da primeira malha termina em Nx*Ny
      lv_if(1,2) = Nx*Ny

      ! Para cada nivel de malha
      DO level = 1, Levels

         ! Calcula o Nx e Ny de cada nivel de malha
         lv_nx(level) = Nx/(2**(level-1))
         lv_ny(level) = Ny/(2**(level-1))

         ! O inicio e fim da primeira malha ja foram calculados
         IF ( level > 1 ) THEN

            ! Calculando o indice que define o fim do vetor
            lv_if(level,2) = lv_if(level-1,2) + lv_nx(level)*lv_ny(level)

            ! Calculando o indice que define o inicio do vetor
            lv_if(level,1) = lv_if(level-1,2) + 1

         END IF

         !write(*,*) "Para P", level, lv_if(level, 1), lv_if(level,2)

      END DO

      ! Alocando as variaveis
      ! lv_if( Levels, 2 ) => comprimento total do array com todas as malhas
      ALLOCATE( u( lv_if( Levels, 2 ), Tn ) , &
                f( lv_if( Levels, 2 ), Tn ) , &
                bp( lv_if( Levels, 2 ), Tn ) &
              )

     ! Para evitar lixo numerico na memoria
     u(:, :) = 0.0d0
     f(:, :) = 0.0d0
     bp(:, :) = 0.0d0

     !-------------------------- matriz A de P --------------------------------
     ALLOCATE( A(Levels) )

     ! Inicializando a matriz A de U
     DO level = 1, Levels

         hx  = Lx/DBLE( lv_nx(level) )
         hy  = Ly/DBLE( lv_ny(level) )

         CALL A(level) % init( lv_nx(level), lv_ny(level), &
                               hx, hy, ht, &
                               Re)

         ! Atribuindo um stencil a cada volume de A
         CALL A(level) % atribui_stencils()

         ! Calculandos os coeficientes e posicoes dos stencils
         CALL A(level) % stencils_p()
         !-------------------------- matriz A de P --------------------------------
     END DO

   END SUBROUTINE Atualiza_Variaveis_MG_GS_P

   SUBROUTINE RESTRICAO_P (fina, grossa, Nxf, Nxg, Nyg)
      IMPLICIT NONE

      INTEGER, INTENT(IN)                           :: Nxf, Nxg, Nyg
      DOUBLE PRECISION, DIMENSION(:), INTENT(IN)    :: fina
      DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT) :: grossa

      ! Declaracao de variaveis interna a esta rotina
      INTEGER :: pg, pf

      DO pg = 1, Nxg * Nyg

         ! Deduzido de if = 2ig - 1, jg = 2jg - 1 e Nf = 2Nxg
         pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1

         ! Fazendo a restricao dos nos internos de acordo com o stencil escolhido
         grossa(pg) = &
        ( fina(pf) + fina(pf + 1 ) + fina(pf + Nxf) + fina(pf + Nxf + 1) )/4.0d0

      END DO

   END SUBROUTINE RESTRICAO_P

   SUBROUTINE Prolongacao_P (fina, grossa, Nxf, Nxg, Nyg )
      IMPLICIT NONE

      INTEGER, INTENT(IN)                            :: Nxf, Nxg, Nyg
      DOUBLE PRECISION, DIMENSION(:), INTENT(IN)     :: grossa
      DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT)  :: fina

      ! Declaracao de variaveis interna a esta rotina
      INTEGER :: p, pf, pg

      ! Volumes no CENTRO recebem stencil 5, mas ja foram considerados
      DO p = 1, (Nxg-2)*(Nyg-2)

         pg = p + 2*FLOOR( DBLE(p-1)/DBLE(Nxg-2) ) + Nxg + 1
         pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1

         fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg - 1 ) )/16.0d0
         fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg + 1 ) )/16.0d0
         fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg - 1 ) )/16.0d0
         fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg + 1 ) )/16.0d0

      END DO

      ! Canto SOUTH WEST
      pg = 1
      pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
      ! Usando Neumann dp/dt = 0 entao p-1 = p e p-nx = p
      fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg       ) + grossa(pg           ) )/16.0d0
      fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg       ) + grossa(pg       + 1 ) )/16.0d0
      fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg     ) )/16.0d0
      fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg + 1 ) )/16.0d0

      ! Contorno SOUTH recebe stencil 2
      DO pg = 2, Nxg - 1
         pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
         ! Usando Neumann dp/dt = 0 entao p - Nx = p
         fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg       ) + grossa(pg       - 1 ) )/16.0d0
         fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg       ) + grossa(pg       + 1 ) )/16.0d0
         fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg - 1 ) )/16.0d0
         fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg + 1 ) )/16.0d0

      END DO

      ! Canto SOUTH EAST  recebe stencil 3
      pg = Nxg
      pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
      ! Usando Neumann dp/dt = 0 entao p+1 = p e p-nx = p
      fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg       ) + grossa(pg       - 1 ) )/16.0d0
      fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg       ) + grossa(pg           ) )/16.0d0
      fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg - 1 ) )/16.0d0
      fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg + 1 ) )/16.0d0

      ! Contorno WEST recebe stencil 4
      DO pg = 1 + Nxg, Nxg *( Nyg - 2 ) + 1, Nxg
          pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
          ! Usando Neumann dp/dt = 0 entao p-1 = p
          fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg     ) )/16.0d0
          fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg + 1 ) )/16.0d0
          fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg     ) )/16.0d0
          fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg + 1 ) )/16.0d0

      END DO

      ! Contorno EAST recebe stencil 6
      DO pg = 2 * Nxg, Nxg * (Nyg - 1), Nxg
            pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
            ! Usando Neumann dp/dt = 0 entao p+1 = p
            fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg - 1 ) )/16.0d0
            fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg     ) )/16.0d0
            fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg - 1 ) )/16.0d0
            fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg + Nxg ) + grossa(pg + Nxg     ) )/16.0d0

      END DO

     ! Canto NORTH WEST recebe stencil 7
     pg = Nxg *( Nyg - 1 ) + 1
     pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
     ! Usando Neumann dp/dt = 0 entao p-1 = p, p+nx = p
     fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg     ) )/16.0d0
     fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg + 1 ) )/16.0d0
     fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg       ) + grossa(pg           ) )/16.0d0
     fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg       ) + grossa(pg       + 1 ) )/16.0d0

     ! Contorno NORTH recebe stencil 8
     DO pg = Nxg *( Nyg - 1 ) + 2, Nxg * Nyg - 1
         pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
         ! Usando Neumann dp/dt = 0 entao p+nx = p
         fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg - 1 ) )/16.0d0
         fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg + 1 ) )/16.0d0
         fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg       ) + grossa(pg       - 1 ) )/16.0d0
         fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg + 1) + 3.0d0 * grossa(pg       ) + grossa(pg       + 1 ) )/16.0d0

     END DO

     ! Canto NORTH EAST recebe stencil 9
     pg = Nxg * Nyg
     pf = 2 * pg + 2 * FLOOR( DBLE(pg-1)/Nxg ) * Nxg - 1
     ! Usando Neumann dp/dt = 0 entao p+1 =p e p+nx = p
     fina(pf)           = fina(pf)           + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg - 1 ) )/16.0d0
     fina(pf+1)         = fina(pf+1)         + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg - Nxg ) + grossa(pg - Nxg     ) )/16.0d0
     fina(pf + Nxf)     = fina(pf + Nxf)     + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg - 1) + 3.0d0 * grossa(pg       ) + grossa(pg       - 1 ) )/16.0d0
     fina(pf + Nxf + 1) = fina(pf + Nxf + 1) + ( 9.0d0 * grossa(pg) + 3.0d0 * grossa(pg    ) + 3.0d0 * grossa(pg       ) + grossa(pg           ) )/16.0d0

   END SUBROUTINE Prolongacao_P

   SUBROUTINE calcula_residuo_p(A, u, b, r)

        USE ClassMatrizA

       IMPLICIT NONE

       TYPE ( MatrizA ), INTENT(IN) :: A
       DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT) :: r
       DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: u, b

       INTEGER :: p, stencil, ite, se ! se = elemento do stencil
       DOUBLE PRECISION :: ax, ux, sm, correcao ! Soma dos componentes do stencil

       DO p = 1, A % Nx * A % Ny

           ! Obtendo o stencil do volume
           stencil = A % v(p)

           ! Para cada coeficiente do stencil da matriz
           ! Lembrando que o primeiro e sempre ap e Gauss-Seidel
           ! e definido como
           ! u(p) = (aw * u(p-1) + ae * u(p+1) + as * u(p-Nx) ... + b(p)) / ap
           ! O codigo abaixo soma aw * u(p-1) + ae * u(p+1) + as * u(p-Nx) ...
           sm = 0.0d0
           DO se = 2, SIZE(A % s(stencil) % p) !(comeca em 2 porque nao considera o ap)

               ! ax =     aw,     ae,       as,      an,     aww,        ann, etc...
               ! ux = u(p-1), u(p+1),  u(p-Nx), u(p+Nx),  u(p-2),  u(p+2*Nx), etc...
               ! Logo
               ! ax * ux = aw * u(p-1), ae * u(p+1), as * u(p-Nx), etc...

               ax = A % s(stencil) % c(se)
               ux = u(p + A % s(stencil) % p(se))

               sm = sm +  ax * ux

           END DO

           ! Obtendo o ap
           ax = A % s(stencil) % c(1)

           ! Somando b(p) e subtraindo ap * up para compor o residuo
           r(p) = sm + b(p) - ax * u(p)

       END DO

   END SUBROUTINE calcula_residuo_p

   SUBROUTINE gs_p( A, u, b, iteracoes )

       USE ClassMatrizA

       IMPLICIT NONE

       TYPE ( MatrizA ), INTENT(IN) :: A
       DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT) :: u
       DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: b
       INTEGER, INTENT(IN) :: iteracoes

       INTEGER :: p, stencil, ite, se ! se = elemento do stencil
       DOUBLE PRECISION :: ax, ux, sm ! Soma dos componentes do stencil

       DO ite = 1, iteracoes

           DO p = 1, A % Nx * A % Ny

               ! Obtendo o stencil do volume
               stencil = A % v(p)

               ! Para cada coeficiente do stencil da matriz
               ! Lembrando que o primeiro e sempre ap e Gauss-Seidel
               ! e definido como
               ! u(p) = (aw * u(p-1) + ae * u(p+1) + as * u(p-Nx) ... + b(p)) / ap
               ! O codigo abaixo soma aw * u(p-1) + ae * u(p+1) + as * u(p-Nx) ...
               sm = 0.0d0
               DO se = 2, SIZE(A % s(stencil) % p) !(comeca em 2 porque nao considera o ap)

                   ! ax =     aw,     ae,       as,      an,     aww,        ann, etc...
                   ! ux = u(p-1), u(p+1),  u(p-Nx), u(p+Nx),  u(p-2),  u(p+2*Nx), etc...
                   ! Logo
                   ! ax * ux = aw * u(p-1), ae * u(p+1), as * u(p-Nx), etc...

                   ax = A % s(stencil) % c(se)
                   ux = u(p + A % s(stencil) % p(se))

                   sm = sm +  ax * ux

                   !if (A % s(stencil) % p(se) .eq. -1) then
                   !write(*,*) p, ax, A % s(stencil) % p(se)
                   !end if
               !    write(*,*) p, ax

               END DO

               ! Obtendo o ap
               ax = A % s(stencil) % c(1)

               ! Somando b(p) e divindo por ap para completar o Gauss-Seidel
               u(p) = (sm + b(p))/ax

              ! write(*,*) p, ax, b(p), u(p)

           END DO
           !write(*,*)
           !call exit()

       END DO

   END SUBROUTINE gs_p

   SUBROUTINE MG_GS_P( u_o, f_o, ciclos_V, t)

      USE Variaveis_Solvers_P

      INTEGER, INTENT(IN) :: ciclos_V, t

      ! u_o = estimativa inicial (pressao)
      DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT) :: u_o

      ! termo fonte ( div (ut, vt) )
      DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: f_o

      INTEGER :: level, ciclo, iteracao

      11 FORMAT( 1X, "Multigrid P => Ciclo:", 1X, I3, 3X, "Parada:", 1X, 1PE16.10, 3X, "Norma L2:", 1X,  1PE16.10, 3X, "(Norma L2)**(1/ciclo)", 1X, 1PE16.10)

      ! Atualizando a primeira malha do vetor u, ou seja,
      !  u (lv_if(1,1):lv_if(1,2)) = u_o(:) => estimativa inicial e a pressao
      !  calculada em um passo de tempo anterior
      u(lv_if(1,1):lv_if(1,2), t) = u_o(:)

      DO ciclo = 1, 1 !max_ciclos_V

         bp(lv_if(1,1):lv_if(1,2), t) = f_o(:)

         CALL calcula_residuo_p ( A(1), &
                                  u(lv_if(1,1):lv_if(1,2), t), &
                                  bp(lv_if(1,1):lv_if(1,2), t), &
                                  f(lv_if(1,1):lv_if(1,2), t) )
         norma_L2 = DSQRT(SUM(f(lv_if(1,1):lv_if(1,2), t) * f(lv_if(1,1):lv_if(1,2), t)))
         normas_L2(ciclo, t) = norma_L2
         !WRITE(*, 11) ciclo, normas_L2(ciclo)/normas_L2(1), normas_L2(ciclo), (normas_L2(ciclo)/normas_L2(1))**(1.0d0/ciclo)

         IF ( normas_L2(ciclo, t)/normas_L2(1, t) < tol ) THEN

            EXIT

         END IF

         DO level = 1, Levels - 1

            ! So atualiza bp = residuo, ou seja, apenas para level > 1 (o primeiro e bp = fonte)
            IF ( level > 1 ) THEN

               bp(lv_if(level,1):lv_if(level,2), t) = f(lv_if(level,1):lv_if(level,2), t)

            END IF

            CALL gs_p ( A(level), u(lv_if(level,1):lv_if(level,2), t), bp(lv_if(level,1):lv_if(level,2), t), ite_down )

            CALL calcula_residuo_p ( A(level), &
                                     u(lv_if(level,1):lv_if(level,2), t), &
                                     bp(lv_if(level,1):lv_if(level,2), t), &
                                     f(lv_if(level,1):lv_if(level,2), t) )

            CALL Restricao_P ( f(lv_if(level,1):lv_if(level,2), t),               &
                               f(lv_if(level+1, 1):lv_if(level+1, 2), t),         &
                               lv_nx(level),                                   &
                               lv_nx(level+1), lv_ny(level+1) )

         END DO

         level = Levels
         bp(lv_if(level,1):lv_if(level,2), t) = f(lv_if(level,1):lv_if(level,2), t)
         CALL gs_p ( A(level), u(lv_if(level,1):lv_if(level,2), t), bp(lv_if(level,1):lv_if(level,2), t), ite_down )

         DO level = Levels, 2, -1

             CALL Prolongacao_P ( u(lv_if(level-1,1):lv_if(level-1,2), t), &
                                  u(lv_if(level, 1):lv_if(level, 2), t),   &
                                  lv_nx(level-1),                       &
                                  lv_nx(level), lv_ny(level) )

             ! Limpando o vetor para evitar erros no proximo ciclo
             u(lv_if(level, 1):lv_if(level, 2), t) = 0.0d0
             CALL gs_p ( A(level-1), u(lv_if(level-1,1):lv_if(level-1,2), t), &
                                    bp(lv_if(level-1,1):lv_if(level-1,2), t), ite_up )

         END DO
      END DO

      u_o(:) = u(lv_if(1,1):lv_if(1,2), t)

   END SUBROUTINE MG_GS_P

   SUBROUTINE Atualiza_Variaveis_MG_GS_U( Lx_o, Ly_o, &
                                          Nx_o, Ny_o, &
                                          tol_o, &
                                          ite_down_o, ite_up_o, &
                                          max_ciclos_V_o, ht, Re, Tn_o )

      USE Variaveis_Solvers_U

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: Nx_o, Ny_o, ite_down_o, ite_up_o, max_ciclos_V_o, Tn_o
      DOUBLE PRECISION, INTENT(IN) :: Lx_o, Ly_o, tol_o, ht, Re

      INTEGER :: level
      DOUBLE PRECISION :: hx, hy

      Lx = Lx_o
      Ly = Ly_o
      Nx = Nx_o
      Ny = Ny_o
      tol = tol_o
      ite_down = ite_down_o
      ite_up   = ite_up_o
      max_ciclos_V = max_ciclos_V_o
      Tn = Tn_o

      ! Calculando o numero de niveis
      Levels = LOG( DBLE(Nx+1))/LOG( DBLE(2) ) - 1

      ! lv_if => duas dimensoes 1: inicio da malha, 2: fim da malha
      ALLOCATE( lv_if(Levels, 2), lv_nx(Levels), lv_ny(Levels), normas_L2(max_ciclos_V, Tn) )

      ! Inicio da primeira malha comeca em 1
      lv_if(1,1) = 1

      ! Fim da primeira malha termina em Nx*Ny
      lv_if(1,2) = Nx*Ny

      ! Para cada nivel de malha
      DO level = 1, Levels

         ! Calcula o Nx e Ny de cada nivel de malha
         lv_nx(level) = (Nx+1)/(2**(level-1)) - 1
         lv_ny(level) = Ny/(2**(level-1))

         ! O inicio e fim da primeira malha ja foram calculados
         IF ( level > 1 ) THEN

            ! Calculando o indice que define o fim do vetor
            lv_if(level,2) = lv_if(level-1,2) + lv_nx(level)*lv_ny(level)

            ! Calculando o indice que define o inicio do vetor
            lv_if(level,1) = lv_if(level-1,2) + 1

         END IF

         !write(*,*) "Para P", level, lv_if(level, 1), lv_if(level,2)

      END DO

      ! Alocando as variaveis
      ! lv_if( Levels, 2 ) => comprimento total do array com todas as malhas
      ALLOCATE( u( lv_if( Levels, 2 ), Tn ) , &
                f( lv_if( Levels, 2 ), Tn ) , &
                bp( lv_if( Levels, 2 ), Tn ) &
              )

     ! Para evitar lixo numerico na memoria
     u(:, :) = 0.0d0
     f(:, :) = 0.0d0
     bp(:, :) = 0.0d0

     !-------------------------- matriz A de P --------------------------------
     ALLOCATE( A(Levels) )

     ! Inicializando a matriz A de U
     DO level = 1, Levels

         hx  = Lx/DBLE( lv_nx(level) + 1 )
         hy  = Ly/DBLE( lv_ny(level) )

         CALL A(level) % init( lv_nx(level), lv_ny(level), &
                               hx, hy, ht, &
                               Re)

         ! Atribuindo um stencil a cada volume de A
         CALL A(level) % atribui_stencils()

         ! Calculandos os coeficientes e posicoes dos stencils
         CALL A(level) % stencils_u(1.0d0)
         !-------------------------- matriz A de P --------------------------------
     END DO

   END SUBROUTINE Atualiza_Variaveis_MG_GS_U

   SUBROUTINE Atualiza_Variaveis_MG_GS_V( Lx_o, Ly_o, &
                                          Nx_o, Ny_o, &
                                          tol_o, &
                                          ite_down_o, ite_up_o, &
                                          max_ciclos_V_o, ht, Re, Tn_o )

      USE Variaveis_Solvers_V

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: Nx_o, Ny_o, ite_down_o, ite_up_o, max_ciclos_V_o, Tn_o
      DOUBLE PRECISION, INTENT(IN) :: Lx_o, Ly_o, tol_o, ht, Re

      INTEGER :: level
      DOUBLE PRECISION :: hx, hy

      Lx = Lx_o
      Ly = Ly_o
      Nx = Nx_o
      Ny = Ny_o
      tol = tol_o
      ite_down = ite_down_o
      ite_up   = ite_up_o
      max_ciclos_V = max_ciclos_V_o
      Tn = Tn_o

      ! Calculando o numero de niveis
      Levels = LOG( DBLE(Nx))/LOG( DBLE(2) ) - 1

      ! lv_if => duas dimensoes 1: inicio da malha, 2: fim da malha
      ALLOCATE( lv_if(Levels, 2), lv_nx(Levels), lv_ny(Levels), normas_L2(max_ciclos_V, Tn) )

      ! Inicio da primeira malha comeca em 1
      lv_if(1,1) = 1

      ! Fim da primeira malha termina em Nx*Ny
      lv_if(1,2) = Nx*Ny

      ! Para cada nivel de malha
      DO level = 1, Levels

         ! Calcula o Nx e Ny de cada nivel de malha
         lv_nx(level) = Nx/(2**(level-1))
         lv_ny(level) = (Ny+1)/(2**(level-1)) - 1

         ! O inicio e fim da primeira malha ja foram calculados
         IF ( level > 1 ) THEN

            ! Calculando o indice que define o fim do vetor
            lv_if(level,2) = lv_if(level-1,2) + lv_nx(level)*lv_ny(level)

            ! Calculando o indice que define o inicio do vetor
            lv_if(level,1) = lv_if(level-1,2) + 1

         END IF

         !write(*,*) "Para P", level, lv_if(level, 1), lv_if(level,2)

      END DO

      ! Alocando as variaveis
      ! lv_if( Levels, 2 ) => comprimento total do array com todas as malhas
      ALLOCATE( u( lv_if( Levels, 2 ), Tn ) , &
                f( lv_if( Levels, 2 ), Tn ) , &
                bp( lv_if( Levels, 2 ), Tn ) &
              )

     ! Para evitar lixo numerico na memoria
     u(:, :) = 0.0d0
     f(:, :) = 0.0d0
     bp(:, :) = 0.0d0

     !-------------------------- matriz A de P --------------------------------
     ALLOCATE( A(Levels) )

     ! Inicializando a matriz A de U
     DO level = 1, Levels

         hx  = Lx/DBLE( lv_nx(level) )
         hy  = Ly/DBLE( lv_ny(level) + 1 )

         CALL A(level) % init( lv_nx(level), lv_ny(level), &
                               hx, hy, ht, &
                               Re)

         ! Atribuindo um stencil a cada volume de A
         CALL A(level) % atribui_stencils()

         ! Calculandos os coeficientes e posicoes dos stencils
         CALL A(level) % stencils_v(1.0d0)
         !-------------------------- matriz A de P --------------------------------
     END DO

   END SUBROUTINE Atualiza_Variaveis_MG_GS_V

END MODULE Subrotinas_MG_GS


text:
C -----------------------------------------------------------------------
C  Program to apply transformation coefficients to a table of
C  instrumental photometry
C
C  ISYA 2004                               (c) Ignasi Ribas (IEEC, Spain)
C -----------------------------------------------------------------------
C   TODO: Possibility of including positional data to keep track
C
      PROGRAM TRANSSTD
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      Open(1,file='coefficients.dat',status='old')
      Read(1,*)
      Read(1,*)
      Read(1,*)
      Read(1,100) coey1,coey2
      Read(1,100) ecoy1,ecoy2
      Read(1,*)
      Read(1,*)
      Read(1,*)
      Read(1,*)
      Read(1,100) coeb1,coeb2
      Read(1,100) ecob1,ecob2
      Close(1)
 100  Format(18x,f7.3,4x,f7.3)

      Open(1,file='stars.lst',status='old')
      Open(2,file='stars_std_pos.lst')
      Write(2,150)
 150  Format('# V_std ','  s_V  ',' (b-y)_s','  s_by ')
  1   Read(1,*,end=2) xceny,yceny,xcenb,ycenb,yinst,ery,binst,erb
      ystd=yinst+coey1+coey2*(binst-yinst)
      bystd=coeb1+coeb2*(binst-yinst)
      erby=dsqrt(erb*erb+ery*ery)
      erys=dsqrt(ery*ery+ecoy1*ecoy1+coey2*coey2*erby*erby+
     $ (binst-yinst)*(binst-yinst)*ecoy2*ecoy2)
      erbys=dsqrt(ecob1*ecob1+coeb2*coeb2*erby*erby+
     $ (binst-yinst)*(binst-yinst)*ecob2*ecob2)
      Write(2,200) ystd,erys,bystd,erbys,xceny,yceny
 200  Format(2(f8.3,f7.3),2f9.2)
      Goto 1
  2   Close(1)

      Stop
      End

max_stars_repo_name: mbadenas/Optical-Astronomy

max_issues_repo_name: mbadenas/Optical-Astronomy

max_forks_repo_name: mbadenas/Optical-Astronomy