text:
*DECK ALBETA
      FUNCTION ALBETA (A, B)
C***BEGIN PROLOGUE  ALBETA
C***PURPOSE  Compute the natural logarithm of the complete Beta
C            function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7B
C***TYPE      SINGLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C)
C***KEYWORDS  FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ALBETA computes the natural log of the complete beta function.
C
C Input Parameters:
C       A   real and positive
C       B   real and positive
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  ALNGAM, ALNREL, GAMMA, R9LGMC, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  ALBETA
      EXTERNAL GAMMA
      SAVE SQ2PIL
      DATA SQ2PIL / 0.9189385332 0467274 E0 /
C***FIRST EXECUTABLE STATEMENT  ALBETA
      P = MIN (A, B)
      Q = MAX (A, B)
C
      IF (P .LE. 0.0) CALL XERMSG ('SLATEC', 'ALBETA',
     +   'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2)
      IF (P.GE.10.0) GO TO 30
      IF (Q.GE.10.0) GO TO 20
C
C P AND Q ARE SMALL.
C
      ALBETA = LOG(GAMMA(P) * (GAMMA(Q)/GAMMA(P+Q)) )
      RETURN
C
C P IS SMALL, BUT Q IS BIG.
C
 20   CORR = R9LGMC(Q) - R9LGMC(P+Q)
      ALBETA = ALNGAM(P) + CORR + P - P*LOG(P+Q) +
     1  (Q-0.5)*ALNREL(-P/(P+Q))
      RETURN
C
C P AND Q ARE BIG.
C
 30   CORR = R9LGMC(P) + R9LGMC(Q) - R9LGMC(P+Q)
      ALBETA = -0.5*LOG(Q) + SQ2PIL + CORR + (P-0.5)*LOG(P/(P+Q))
     1  + Q*ALNREL(-P/(P+Q))
      RETURN
C
      END


text:
!======================================================================*
SUBROUTINE RISort(X,Y,N,KFLAG)
    !----------------------------------------------------------------------!
    IMPLICIT NONE
    !----------------------------------------------------------------------!
    !-------------------------------[Locals]-------------------------------!
    INTEGER             :: n, i, ij, il, iu, j, k, kk, l, m, nn, kflag
    REAL                :: r
    !----------------------------------------------------------------------*
    ! Sorts real array X and make same interchanges in integer array Y.    *
    !----------------------------------------------------------------------*
    !-------------------------------[Locals]-------------------------------!
    DIMENSION IL(21),IU(21)
    REAL    X(N),T, TT
    INTEGER Y(N),TY,TTY
    !--------------------------------[CVS]---------------------------------*
    !      $Id: RISort.f90,v 1.3 2002/10/30 16:29:33 niceno Exp $
    !      $Source: /home/muhamed/.CVSROOT/T-Rex/Library/RISort.f90,v $
    !======================================================================*
    !***FIRST EXECUTABLE STATEMENT RISORT
    NN = N
    IF (NN.GE.1) GO TO 10
    RETURN
10  KK = IABS(KFLAG)
    IF ((KK.EQ.1).OR.(KK.EQ.2)) GO TO 15
    RETURN
    !
    ! ALTER ARRAY X TO GET DECREASING ORDER IF NEEDED
    !
15  IF (KFLAG.GE.1) GO TO 30
    DO 20 I=1,NN
        X(I) = -X(I)
20  CONTINUE
30  if(KK==1) goto 100
    if(KK==2) goto 200
!
! SORT X ONLY
!
100 CONTINUE
    M=1
    I=1
    J=NN
    R=.375
110 IF (I .EQ. J) GO TO 155
    IF (R .GT. .5898437) GO TO 120
    R=R+3.90625E-2
    GO TO 125
120 R=R-.21875
125 K=I
    !                                  SELECT A CENTRAL ELEMENT OF THE
    !                                  ARRAY AND SAVE IT IN LOCATION T
    IJ = I + IFIX (FLOAT (J-I) * R)
    T=X(IJ)
    !                                  IF FIRST ELEMENT OF ARRAY IS GREATER
    !                                  THAN T, INTERCHANGE WITH T
    IF (X(I) .LE. T) GO TO 130
    X(IJ)=X(I)
    X(I)=T
    T=X(IJ)
130 L=J
    !                                  IF LAST ELEMENT OF ARRAY IS LESS THAN
    !                                  T, INTERCHANGE WITH T
    IF (X(J) .GE. T) GO TO 140
    X(IJ)=X(J)
    X(J)=T
    T=X(IJ)
    !                                  IF FIRST ELEMENT OF ARRAY IS GREATER
    !                                  THAN T, INTERCHANGE WITH T
    IF (X(I) .LE. T) GO TO 140
    X(IJ)=X(I)
    X(I)=T
    T=X(IJ)
    GO TO 140
135 TT=X(L)
    X(L)=X(K)
    X(K)=TT
    !                                  FIND AN ELEMENT IN THE SECOND HALF OF
    !                                  THE ARRAY WHICH IS SMALLER THAN T
140 L=L-1
    IF (X(L) .GT. T) GO TO 140
    !                                  FIND AN ELEMENT IN THE FIRST HALF OF
    !                                  THE ARRAY WHICH IS GREATER THAN T
145 K=K+1
    IF (X(K) .LT. T) GO TO 145
    !                                  INTERCHANGE THESE ELEMENTS
    IF (K .LE. L) GO TO 135
    !                                  SAVE UPPER AND LOWER SUBSCRIPTS OF
    !                                  THE ARRAY YET TO BE SORTED
    IF (L-I .LE. J-K) GO TO 150
    IL(M)=I
    IU(M)=L
    I=K
    M=M+1
    GO TO 160
150 IL(M)=K
    IU(M)=J
    J=L
    M=M+1
    GO TO 160
    !                                  BEGIN AGAIN ON ANOTHER PORTION OF
    !                                  THE UNSORTED ARRAY
155 M=M-1
    IF (M .EQ. 0) GO TO 300
    I=IL(M)
    J=IU(M)
160 IF (J-I .GE. 1) GO TO 125
    IF (I .EQ. 1) GO TO 110
    I=I-1
165 I=I+1
    IF (I .EQ. J) GO TO 155
    T=X(I+1)
    IF (X(I) .LE. T) GO TO 165
    K=I
170 X(K+1)=X(K)
    K=K-1
    IF (T .LT. X(K)) GO TO 170
    X(K+1)=T
    GO TO 165
!
! SORT X AND CARRY Y ALONG
!
200 CONTINUE
    M=1
    I=1
    J=NN
    R=.375
210 IF (I .EQ. J) GO TO 255
    IF (R .GT. .5898437) GO TO 220
    R=R+3.90625E-2
    GO TO 225
220 R=R-.21875
225 K=I
    !                                  SELECT A CENTRAL ELEMENT OF THE
    !                                  ARRAY AND SAVE IT IN LOCATION T
    IJ = I + IFIX (FLOAT (J-I) *R)
    T=X(IJ)
    TY= Y(IJ)
    !                                  IF FIRST ELEMENT OF ARRAY IS GREATER
    !                                  THAN T, INTERCHANGE WITH T
    IF (X(I) .LE. T) GO TO 230
    X(IJ)=X(I)
    X(I)=T
    T=X(IJ)
    Y(IJ)= Y(I)
    Y(I)=TY
    TY= Y(IJ)
230 L=J
    !                                  IF LAST ELEMENT OF ARRAY IS LESS THAN
    !                                  T, INTERCHANGE WITH T
    IF (X(J) .GE. T) GO TO 240
    X(IJ)=X(J)
    X(J)=T
    T=X(IJ)
    Y(IJ)= Y(J)
    Y(J)=TY
    TY= Y(IJ)
    !                                  IF FIRST ELEMENT OF ARRAY IS GREATER
    !                                  THAN T, INTERCHANGE WITH T
    IF (X(I) .LE. T) GO TO 240
    X(IJ)=X(I)
    X(I)=T
    T=X(IJ)
    Y(IJ)= Y(I)
    Y(I)=TY
    TY= Y(IJ)
    GO TO 240
235 TT=X(L)
    X(L)=X(K)
    X(K)=TT
    TTY= Y(L)
    Y(L)= Y(K)
    Y(K)=TTY
    !                                  FIND AN ELEMENT IN THE SECOND HALF OF
    !                                  THE ARRAY WHICH IS SMALLER THAN T
240 L=L-1
    IF (X(L) .GT. T) GO TO 240
    !                                  FIND AN ELEMENT IN THE FIRST HALF OF
    !                                  THE ARRAY WHICH IS GREATER THAN T
245 K=K+1
    IF (X(K) .LT. T) GO TO 245
    !                                  INTERCHANGE THESE ELEMENTS
    IF (K .LE. L) GO TO 235
    !                                  SAVE UPPER AND LOWER SUBSCRIPTS OF
    !                                  THE ARRAY YET TO BE SORTED
    IF (L-I .LE. J-K) GO TO 250
    IL(M)=I
    IU(M)=L
    I=K
    M=M+1
    GO TO 260
250 IL(M)=K
    IU(M)=J
    J=L
    M=M+1
    GO TO 260
    !                                  BEGIN AGAIN ON ANOTHER PORTION OF
    !                                  THE UNSORTED ARRAY
255 M=M-1
    IF (M .EQ. 0) GO TO 300
    I=IL(M)
    J=IU(M)
260 IF (J-I .GE. 1) GO TO 225
    IF (I .EQ. 1) GO TO 210
    I=I-1
265 I=I+1
    IF (I .EQ. J) GO TO 255
    T=X(I+1)
    TY= Y(I+1)
    IF (X(I) .LE. T) GO TO 265
    K=I
270 X(K+1)=X(K)
    Y(K+1)= Y(K)
    K=K-1
    IF (T .LT. X(K)) GO TO 270
    X(K+1)=T
    Y(K+1)=TY
    GO TO 265
    !
    ! CLEAN UP
    !
300 IF (KFLAG.GE.1) RETURN
    DO 310 I=1,NN
        X(I) = -X(I)
310 CONTINUE
    RETURN
END SUBROUTINE

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
subroutine sup_endste(a,itask)
   !This routine ends a time step of the incoma%pressible NS equations.
   use typre
   use def_parame
   use Mod_ThreeField
   use Mod_Mesh
   use Mod_memor
   use Mod_SupExacso
   implicit none
   class(ThreeFieldNSProblem) :: a
   type(FemMesh)    :: Mesh
   type(MemoryMan)  :: Memor
   type(SupExacso)  :: exacso

   integer(ip) :: itask
   integer(ip) :: ifout,ndime,npoin,ipoin,idime,ntens,auxstep,icomp,j
   integer(ip) :: ielem,nelem
   !Exact Values
   real(rp), allocatable   :: exvel(:),exveg(:,:),exprg(:),exsig(:),exsigr(:,:)
   real(rp)                :: expre
   !Nodal coordinates
   real(rp), pointer       :: coord(:)
   !For log-conformation model
   real(rp), allocatable   :: Expo_Vect(:), Vect(:), Expo_Matrix(:,:), Matrix(:,:), psiNew(:), Vect2(:)
   real(rp)   ::  lambda0,acvis,auxL,beta,lamb, auxL2

   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)
   ntens=(ndime-1)*(ndime-1)+2

   auxstep=3

   !Pre-CrankNicolson
   if (itask == 1) then
      !Compute dissipation and other endste elemental computations using u^(n+theta)
      !Also compute the dynamic subgrid scales
      call a%EnditeElmope('Endste')

   !Post-CrankNicolson
   elseif (itask == 2) then

      !When log-conformation model is employed a%sigma is the new variable psi.
      !We have to calculate initial sigma  (now called sigmaold) to save this field too.

      call a%Memor%alloc(ndime,ndime,Matrix,'Matrix','output')
      call a%Memor%alloc(ntens,Vect,'Vect','output')
      call a%Memor%alloc(ntens,Vect2,'Vect2','output')
      call a%Memor%alloc(ntens,Expo_Vect,'ExpoVect','output')
      call a%Memor%alloc(ntens,PsiNew,'PsiNew','output')
      call a%Memor%alloc(ndime,ndime,Expo_Matrix,'Expo_Matrix','output')
      call a%Mesh%GetNpoin(npoin)

      call a%IncrementalLambda(1,a%MatProp(1)%LawViParam(3),lamb)
      beta=a%MatProp(1)%LawViParam(2)
      acvis=a%MatProp(1)%LawViParam(1)

      if (a%LogFormulation==1) then

         if (a%MatProp(1)%LawViParam(3)==0 .or. a%MatProp(1)%LawViParam(5)==0) then
            lambda0=0.01_rp
         else
            if(a%kfl_exacs/=0) then
               lambda0=a%MatProp(1)%LawViParam(5) !SETLAMBDA0
            else
               if(a%kfl_timei==0) then
                  lambda0=lamb*a%MatProp(1)%LawViParam(5)
               else
                  lambda0=a%MatProp(1)%LawViParam(5)*a%MatProp(1)%LawViParam(3)
               end if
            end if
         end if
         auxL=((1.0_rp-beta)*acvis)/lambda0
         auxL2=((1.0_rp-beta)*acvis)/lamb

         do ipoin=1,npoin
            Vect(:)=a%sigma(:,ipoin,1)
            call sup_ComputeExponential(ndime,ntens,Vect,Expo_Vect,Expo_Matrix)

            call deter(Expo_Matrix,a%taudet(ipoin),ndime)!DETERMINANT OF CONFORMATION TENSOR

            do j=1,ndime
               Expo_Vect(j)=Expo_Vect(j)-1.0_rp
            end do
            a%sigmaold(:,ipoin)=auxL*Expo_Vect(:)

            !If it is required, the variable psi "real" is calculated
            !- in order to compare this field with other publications

            Vect2(:)=auxL*Expo_Vect(:)
            call sup_LogarithmicTransformation(ndime,ntens,auxL2,Vect2,psiNew)
            if(a%npp_stepi(28)>0) a%psiReal(:,ipoin)=psiNew(:)
            !Expo_Matrix = tau_matrix
         end do
      else if (a%LogFormulation==0) then     !IF THE STANDARD MODEL IS SET
         auxL=lamb/((1.0_rp-beta)*acvis)
         call a%Mesh%GetNpoin(npoin)
         do ipoin=1,npoin
            Vect(:)=a%sigma(:,ipoin,1)
            Vect(:)=(1.0_rp/auxL)*Vect(:)
            do j=1,ndime
               Vect(j)=Vect(j)+1.0_rp
            end do
            !pasar vector a matrix.
            call sup_SigmaMatrix(ndime,ntens,Vect,Matrix)
            call deter(Matrix,a%taudet(ipoin),ndime)
         end do
      end if

      call a%Memor%dealloc(ndime,ndime,Matrix,'Matrix','output')
      call a%Memor%dealloc(ntens,Vect,'Vect','output')
      call a%Memor%dealloc(ntens,Vect2,'Vect2','output')
      call a%Memor%dealloc(ntens,Expo_Vect,'ExpoVect','output')
      call a%Memor%dealloc(ntens,PsiNew,'PsiNew','output')
      call a%Memor%dealloc(ndime,ndime,Expo_Matrix,'Expo_Matrix','output')

      !Update a%velocity and a%pressure
      if (a%kfl_exacs==0) then
         !Higher order components
         if (a%ncomp > 3) then
            do icomp = a%ncomp, 4,-1
               a%veloc(:,:,icomp) = a%veloc(:,:,icomp-1)
               a%press(:,icomp) = a%press(:,icomp-1)
               a%sigma(:,:,icomp) = a%sigma(:,:,icomp-1)
            enddo
         endif
         !Previous time step component
         a%veloc(:,:,3) = a%veloc(:,:,1)  ! Vn = Vn+1
         a%press(:,3) = a%press(:,1)
         a%sigma(:,:,3) = a%sigma(:,:,1)

         !Boundary operation at the end
         call a%EndBouope('Endste')

         !Update the subgrid a%velocity DynamicSS
         if(a%kfl_tacsg/=0) then
            call a%Mesh%GetNelem(nelem)
            if (a%kfl_tacsg>0) then
               !Crank Nicolson schemes
               do ielem=1,nelem
                  a%vesgs(ielem)%a(:,2,:)=a%vesgs(ielem)%a(:,1,:)
                  if (a%MatProp(1)%lawvi<0) then
                     a%sisgs(ielem)%a(:,2,:)=a%sisgs(ielem)%a(:,1,:)
                  end if
                  if (a%kfl_repro>=2) then
                     a%vesgs2(ielem)%a(:,2,:)=a%vesgs2(ielem)%a(:,1,:)
                     a%vesgs3(ielem)%a(:,2,:)=a%vesgs3(ielem)%a(:,1,:)
                  end if
                  if (a%kfl_repro == 4 .and. a%MatProp(1)%lawvi<0) then
                     a%sisgs2(ielem)%a(:,2,:)=a%sisgs2(ielem)%a(:,1,:)
                     a%sisgs3(ielem)%a(:,2,:)=a%sisgs3(ielem)%a(:,1,:)
                  end if
               end do
            endif
         end if

      ! transient exact solution
      elseif(a%kfl_exacs/=0.and.a%kfl_timei==1)then

         if (a%ncomp > 3) then
            do icomp = a%ncomp, 4,-1
               a%veloc(:,:,icomp) = a%veloc(:,:,icomp-1)
               a%press(:,icomp)   = a%press(:,icomp-1)
               a%sigma(:,:,icomp) = a%sigma(:,:,icomp-1)
            enddo
         endif

         ! Allocate exact components
         call a%Memor%alloc(ndime,exvel,'exvel','sup_endste')
         call a%Memor%alloc(ndime,ndime,exveg,'exveg','sup_endste')
         call a%Memor%alloc(ntens,exsig,'exsig','sup_endste')
         call a%Memor%alloc(ntens,ndime,exsigr,'exsigr','sup_endste')
         call a%Memor%alloc(ndime,exprg,'exprg','sup_endste')

         do ipoin = 1,npoin

            a%veloc(:,:,3) = a%veloc(:,:,1)  ! Vn = Vn+1
            a%press(:,3) = a%press(:,1)
            a%sigma(:,:,3)=a%sigma(:,:,1)

            call a%Mesh%GetPointCoord(ipoin,coord)

            ! Exact solution
            call exacso%sup_ComputeSolution(ndime,coord,a%ctime,a%LogFormulation,a)
            call exacso%sup_GetPressure(ndime,expre,exprg)
            call exacso%sup_GetVelocity(ndime,exvel,exveg)
            if (a%LogFormulation/=0) then
               call exacso%sup_GetPsi(ndime,exsig,exsigr)
            else
               call exacso%sup_GetStress(ndime,a%LogFormulation,exsig,exsigr)
            end if

            if(a%istep<=auxstep)then
               do idime=1,ndime
                  a%veloc(idime,ipoin,1) = exvel(idime)  ! Vn-2 = Vn-1
               end do
               do idime=1,ntens
                  a%sigma(idime,ipoin,1) = exsig(idime)
               end do
            end if

            !exact solution
            do idime=1,ndime
               a%exveloc(idime,ipoin,1) = exvel(idime)
            end do
            do idime=1,ntens
               a%exsigma(idime,ipoin,1) = exsig(idime)
            end do
         end do
         ! Allocate exact components
         call a%Memor%dealloc(ndime,exvel,'exvel','sup_endste')
         call a%Memor%dealloc(ndime,ndime,exveg,'exveg','sup_endste')
         call a%Memor%dealloc(ntens,exsig,'exsig','sup_endste')
         call a%Memor%dealloc(ntens,ndime,exsigr,'exsigr','sup_endste')
         call a%Memor%dealloc(ndime,exprg,'exprg','sup_endste')
      end if

   endif

end subroutine sup_endste

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
subroutine stepout(nstatsteps)

use vars
use rad, only: qrad
use sgs, only: tk
use tracers
use microphysics, only: micro_print, micro_statistics
use sgs, only: sgs_print, sgs_statistics
use stat_moments, only: compmoments
use movies, only: mvmovies
use params
use hbuffer
use instrument_diagnostics, only : isccp_write, modis_write, misr_write, zero_instr_diag
implicit none

integer i,j,k,ic,jc,nstatsteps
real div, divmax, divmin
real rdx, rdy, rdz, coef
integer im,jm,km
real wmax, qnmax(1), qnmax1(1)
real(8) buffer(5), buffer1(5)

if(mod(nstep,nstatis).eq.0) then
      call statistics()
      call micro_statistics()
      call sgs_statistics()
end if

if(mod(nstep,nmovie).eq.0.and.nstep.ge.nmoviestart &
                                   .and.nstep.le.nmovieend) then
      call mvmovies()
endif

if(mod(nstep,nstatmom).eq.0.and.nstep.ge.nstatmomstart &
                                   .and.nstep.le.nstatmomend) then
      call compmoments()
endif

if(mod(nstep,nstat).eq.0) then
  if(masterproc) print *,'Writting statistics:nstatsteps=',nstatsteps

  call t_startf ('stat_out')

  ! write out instrument simulator data and compute domain means for statistics
  !bloss(2016-02-06): This needs to be called before hbuf_avg for the simulator
  !  cloud fraction to be recorded at the first statistics output.
  call isccp_write()
  call modis_write()
  call misr_write()

  call hbuf_average(nstatsteps)
  call hbuf_write(nstatsteps)
  call hbuf_flush()
  nstatsteps = 0

  call zero_instr_diag()

  call t_stopf ('stat_out')

endif
if(mod(nstep,nstat*(1+nrestart_skip)).eq.0.or.nstep.eq.nstop.or.nelapse.eq.0) then

  call write_all() ! save restart file

end if

call t_startf ('2D_out')

if(mod(nstep,nsave2D).eq.0.and.nstep.ge.nsave2Dstart &
                                   .and.nstep.le.nsave2Dend) then
  call write_fields2D()
endif

if(.not.save2Davg.or.nstep.eq.nsave2Dstart-nsave2D) call stat_2Dinit(0) ! argument of 0 means storage terms for statistics are preserved

call t_stopf ('2D_out')

call t_startf ('3D_out')
if(mod(nstep,nsave3D).eq.0.and.nstep.ge.nsave3Dstart.and.nstep.le.nsave3Dend ) then
  ! determine if the maximum cloud water exceeds the threshold
  ! value to save 3D fields:
  qnmax(1)=0.
  do k=1,nzm
    do j=1,ny
      do i=1,nx
         qnmax(1) = max(qnmax(1),qcl(i,j,k))
         qnmax(1) = max(qnmax(1),qci(i,j,k))
      end do
    enddo
  enddo
  if(dompi) then
     call task_max_real(qnmax(1),qnmax1(1),1)
     qnmax(1) = qnmax1(1)
  end if
  if(qnmax(1).ge.qnsave3D) then
    call write_fields3D()
  end if
endif
call t_stopf ('3D_out')

!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
! Print stuff out:

call t_startf ('print_out')

if(masterproc) write(*,'(a,i10,a,i3,a,f6.3,a,f6.3)') 'NSTEP = ',nstep,'   NCYCLE=',ncycle, &
                  '  CFL_adv=',cfl_adv,'  CFL_sgs=',cfl_sgs

if(mod(nstep,nprint).eq.0) then

 divmin=1.e20
 divmax=-1.e20

 rdx = 1./dx
 rdy = 1./dy

 wmax=0.
 do k=1,nzm
  coef = rho(k)*adz(k)*dz
  rdz = 1./coef
  if(ny.ne.1) then
   do j=1,ny-1*YES3D
    jc = j+1*YES3D
    do i=1,nx-1
     ic = i+1
     div = (u(ic,j,k)-u(i,j,k))*rdx + (v(i,jc,k)-v(i,j,k))*rdy + &
		  (w(i,j,k+1)*rhow(k+1)-w(i,j,k)*rhow(k))*rdz
     divmax = max(divmax,div)
     divmin = min(divmin,div)
     if(w(i,j,k).gt.wmax) then
	wmax=w(i,j,k)
	im=i
	jm=j
	km=k
     endif
    end do
   end do
  else
    j = 1
    do i=1,nx-1
    ic = i+1
     div = (u(ic,j,k)-u(i,j,k))*rdx +(w(i,j,k+1)*rhow(k+1)-w(i,j,k)*rhow(k))*rdz
     divmax = max(divmax,div)
     divmin = min(divmin,div)
     if(w(i,j,k).gt.wmax) then
	wmax=w(i,j,k)
	im=i
	jm=j
	km=k
     endif
    end do
  endif
 end do

 if(dompi) then
   buffer(1) = total_water_before
   buffer(2) = total_water_after
   buffer(3) = total_water_evap
   buffer(4) = total_water_prec
   buffer(5) = total_water_ls
   call task_sum_real8(buffer, buffer1,5)
   total_water_before = buffer1(1)
   total_water_after = buffer1(2)
   total_water_evap = buffer1(3)
   total_water_prec = buffer1(4)
   total_water_ls = buffer1(5)
 end if

!print*,rank,minval(u(1:nx,1:ny,:)),maxval(u(1:nx,1:ny,:))
!print*,rank,'min:',minloc(u(1:nx,1:ny,:))
!print*,rank,'max:',maxloc(u(1:nx,1:ny,:))

!if(masterproc) then

!print*,'--->',tk(27,1,1)
!print*,'tk->:'
!write(6,'(16f7.2)')((tk(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'p->:'
!write(6,'(16f7.2)')((p(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'u->:'
!write(6,'(16f7.2)')((u(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'v->:'
!write(6,'(16f7.2)')((v(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'w->:'
!write(6,'(16f7.2)')((w(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'qcl:'
!write(6,'(16f7.2)')((qcl(i,1,k)*1000.,i=16,31),k=nzm,1,-1)
!print*,'qpl->:'
!write(6,'(16f7.2)')((qpl(i,1,k)*1000.,i=16,31),k=nzm,1,-1)
!print*,'qcl:->'
!write(6,'(16f7.2)')((qci(i,1,k)*1000.,i=16,31),k=nzm,1,-1)
!print*,'qpl->:'
!write(6,'(16f7.2)')((qpi(i,1,k)*1000.,i=16,31),k=nzm,1,-1)
!print*,'qrad->:'
!write(6,'(16f7.2)')((qrad(i,1,k)*3600.,i=16,31),k=nzm,1,-1)
!print*,'qv->:'
!write(6,'(16f7.2)')((qv(i,1,k)*1000.,i=16,31),k=nzm,1,-1)
!print*,'t->:'
!write(6,'(16f7.2)')((t(i,1,k),i=1,16),k=nzm,1,-1)
!print*,'tabs->:'
!write(6,'(16f7.2)')((tabs(i,1,k),i=16,31),k=nzm,1,-1)

!end if

!--------------------------------------------------------
 if(masterproc) then

    print*,'DAY = ',day
    write(6,*) 'NSTEP=',nstep
    write(6,*) 'div:',divmax,divmin
    if(.not.dodynamicocean) write(6,*) 'SST=',tabs_s
    write(6,*) 'surface pressure=',pres0

 endif

 call fminmax_print('u:',u,dimx1_u,dimx2_u,dimy1_u,dimy2_u,nzm)
 call fminmax_print('v:',v,dimx1_v,dimx2_v,dimy1_v,dimy2_v,nzm)
 call fminmax_print('w:',w,dimx1_w,dimx2_w,dimy1_w,dimy2_w,nz)
 call fminmax_print('p:',p,0,nx,1-YES3D,ny,nzm)
 call fminmax_print('t:',t,dimx1_s,dimx2_s,dimy1_s,dimy2_s,nzm)
 call fminmax_print('tabs:',tabs,1,nx,1,ny,nzm)
 call fminmax_print('qv:',qv,1,nx,1,ny,nzm)
 if(dosgs) call sgs_print()
 if(docloud) then
   call fminmax_print('qcl:',qcl,1,nx,1,ny,nzm)
   call fminmax_print('qci:',qci,1,nx,1,ny,nzm)
   call micro_print()
 end if
 if(doprecip) then
   call fminmax_print('qpl:',qpl,1,nx,1,ny,nzm)
   call fminmax_print('qpi:',qpi,1,nx,1,ny,nzm)
 end if
 if(dolongwave.or.doshortwave) call fminmax_print('qrad(K/day):',qrad*86400.,1,nx,1,ny,nzm)
 if(dotracers) then
   do k=1,ntracers
     call fminmax_print(trim(tracername(k))//':',tracer(:,:,:,k),dimx1_s,dimx2_s,dimy1_s,dimy2_s,nzm)
   end do
 end if
 call fminmax_print('shf:',fluxbt*cp*rhow(1),1,nx,1,ny,1)
 call fminmax_print('lhf:',fluxbq*lcond*rhow(1),1,nx,1,ny,1)
 call fminmax_print('uw:',fluxbu,1,nx,1,ny,1)
 call fminmax_print('vw:',fluxbv,1,nx,1,ny,1)
 call fminmax_print('sst:',sstxy+t00,0,nx,1-YES3D,ny,1)

 total_water_before = total_water_before/float(nx_gl*ny_gl)
 total_water_after = total_water_after/float(nx_gl*ny_gl)
 total_water_evap = total_water_evap/float(nx_gl*ny_gl)
 total_water_prec = total_water_prec/float(nx_gl*ny_gl)
 total_water_ls = total_water_ls/float(nx_gl*ny_gl)

 if(masterproc) then

   print*,'total water budget:'
   write(*,991) total_water_before !'before (mm):    ',total_water_before
   write(*,992) total_water_after !'after (mm) :    ',total_water_after
   write(*,993) total_water_evap !'evap (mm/day):  ',total_water_evap
   write(*,994) total_water_prec !'prec (mm/day):  ',total_water_prec
   write(*,995) total_water_ls !'ls (mm/day):    ',total_water_ls
   write(*,996) (total_water_after-(total_water_before+total_water_evap+total_water_ls-total_water_prec))
   991 format(' before (mm):       ',F16.11)
   992 format(' after (mm):        ',F16.11)
   993 format(' evaporation (mm):  ',F16.11)
   994 format(' precipitation (mm):',F16.11)
   995 format(' large-scale (mm):  ',F16.11)
   996 format(' Imbalance (mm)    ',F16.11)
   print*,' imbalance (rel error):', &
     (total_water_after-(total_water_before+total_water_evap+total_water_ls-total_water_prec))/total_water_after
   print*,'evap (mm/day):',total_water_evap/dt*86400.
   print*,'prec (mm/day):',total_water_prec/dt*86400.
   print*,'ls (mm/day):',total_water_ls/dt*86400.
   print*,'imbalance (mm/day)', &
     (total_water_after-total_water_before-total_water_evap-total_water_ls+total_water_prec)/dt*86400.

 end if

end if ! (mod(nstep,nprint).eq.0)

call t_stopf ('print_out')

end

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
      subroutine wraca(r,th0,thf)
C.....subroutine wraca 01/06/75
      common/gfgel/gf(28)
      common/gfbftr/xtr,ytr
      if(abs(r).lt.1.e-6)  goto 3
      gf(10)=100.
      r1=abs(r)
      i=abs(gf(9))
      gf(i)=20.
      gf(i+1)=xtr-r1*cos(th0*0.01745329)
      gf(i+2)=ytr-r1*sin(th0*0.01745329)
      gf(i+3)=r
      gf(i+6)=th0
      gf(i+7)=thf
      if(gf(9).gt.0.) call strmod(0)
      return
 3    call gffals(5)
      return
      end

ext: f

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
program nast
  use input_module, only : read_file
  use driver_module
  use write_nast_data
  implicit none

  double precision  :: t_1, t_2
  character(len=80) :: input_file

  open(unit=66, file = 'nast.out')
  call get_command_argument(1,  input_file)
  if (len_trim(input_file) == 0) then
     write (66,*) 'No input file'
     call exit
  end if

  write(66,'(A,A)') 'Input file: ', input_file
  open (unit=11, file=input_file )

  call cpu_time(t_1)

  call write_header()
  call read_file()
  call write_options()
  call run()

  call cpu_time(t_2)
  write(66,'(/,1X,A,F7.2)') 'Total CPU time = ',t_2-t_1
  write(66,'(1X,A)') 'NAST terminated now.'

  close (11)
  close (66)

end program nast

max_stars_repo_name: svarganov/NAST

max_issues_repo_name: svarganov/NAST

max_forks_repo_name: svarganov/NAST


text:
!
!===============================================================================
!-- Common Module Tests
!> Tests for common_mod module, primairly global variables
!===============================================================================
!/****/h* modules|common/test_common_mod
!* SYNOPSIS
MODULE test_common_mod
!* PURPOSE
!*   Tests for common_mod module, primairly global variables
!* INCLUDES
!*   common_mod
!* CONTAINS
!*   Subroutine              Purpose
!*   test_vars(cnt)          - deallocated global structures/variables
!* AUTHOR
!*   Pieter Boom
!* MODIFICATION HISTORY
!*   2019/08/21: Created (PB)
!* COPYRIGHT
!*   (c) University of Manchester
!******/
!> author: Pieter Boom
!> date: 2019/08/21
!> Tests for common_mod module, primairly global variables
!===============================================================================

  !-----------------------------------------------------------------------------
  ! use statements and implicit none
  !-----------------------------------------------------------------------------
  USE common_mod

  IMPLICIT NONE

!===============================================================================
CONTAINS
!===============================================================================
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! /****s*/ /src/common/common_test|test_vars
!* SYNOPSIS
  SUBROUTINE test_vars(cnt)
!* PURPOSE
!*   test global structures/variables
!* INPUTS
!*   test counter
!* OUTPUTS
!*   test counter
!* AUTHOR
!*   Pieter Boom
!* MODIFICATION HISTORY
!*   2019/08/21: Created (PB)
!* COPYRIGHT
!*   (c) University of Manchester
!******/
!> author: Pieter Boom
!> date: 2019/08/21
!> deallocate global structures/variables
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! TO DO:
! - global variables?
! - primary element data structure
! - clean up
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

    !---------------------------------------------------------------------
    ! arguments
    !---------------------------------------------------------------------
    INTEGER, INTENT(INOUT) :: cnt             !> test counter

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !---------------------------------------------------------------------
    ! 1. precision
    !---------------------------------------------------------------------
    cnt = cnt + 1
    IF (1.d0 - 9.999999999999999d-1 > 0.d0) THEN
      WRITE(*,'(A,I5,A)')'common_test : ',cnt,&
      ' : v/      PASSED : precision accurate to 16 decimal places'
    ELSE
      WRITE(*,'(A,I5,A)')'common_test : ',cnt,&
      ' : XXXXXX  FAILED : precision NOT accurate to 16 decimal places'
    END IF

    !---------------------------------------------------------------------
    ! 2. global variables
    ! 3. global mpi
    ! 4. global io
    ! 5. global mesh
    ! 6. primary element
    ! 7. variable clean up
    !---------------------------------------------------------------------

  END SUBROUTINE
! common_test|test_vars
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

END MODULE
! common_test
!===============================================================================
!

max_stars_repo_name: ParaFEM/ParaGEMS

max_issues_repo_name: ParaFEM/ParaGEMS

max_forks_repo_name: ParaFEM/ParaGEMS


text:
*----------------------------------------------------------------------*
      subroutine gen_contr4(strict,
     &                      form_list,proto_main,
     &                      fix_vtx_in,occ_vtx,ol_map,op_info)
*----------------------------------------------------------------------*
*     generate all contraction originating from a proto-contraction
*     on input: proto_main - proto-contraction, i.e. a contraction with
*                         - resulting operator
*                         - all vertices in desired sequence
*                           (including "open line vertices" see
*                            description of occ_vtx; they can be
*                            identified from their operator number which
*                            must be the same as the result operator)
*                         - all arcs with connections which are fixed
*                           (including fixed external lines, given as
*                            usual arcs, however. They will be converted
*                            to xarcs, as soon as the actual contraction
*                            is generated).
*                         - all arcs which must be non-zero (i.e.
*                              the corresp. vertices must be connected)
*               occ_vtx - for convenience: occupations of all vertices
*                         in the sequence conforming with proto_main
*                         the open line vertices have to be set to the
*                         adjoint of the corresponding open line occu-
*                         pation such that we can treat it as if it
*                         were an operator vertex
*               ol_map  - array which tells us whether the given
*                         vertex is no open line (0)
*                         or to which vertex of the result
*                         it contributes (1...njoined_result)
*               op_info - needed for debugging only
*     on output: all new contractions appended to form-list
*
*     extension to more general targets
*
*     andreas, march 2008
*
*----------------------------------------------------------------------*
      implicit none

      integer, parameter ::
     &     ntest = 00

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
c      include 'def_operator.h'
      include 'mdef_operator_info.h'
      include 'def_formula_item.h'
      include 'ifc_operators.h'

      type(formula_item), intent(in), target ::
     &     form_list
      type(contraction), intent(in) ::
     &     proto_main
      logical, intent(in) ::
     &     strict, fix_vtx_in(proto_main%nvtx)
      integer, intent(in) ::
     &     occ_vtx(ngastp,2,proto_main%nvtx),
     &     ol_map(proto_main%nvtx)
      type(operator_info), intent(in) ::
     &     op_info

      type(formula_item), pointer ::
     &     form_pnt
      integer ::
     &     nvtx, ivtx, ij_tgt,
c dbg
     &     occ_test(ngastp,2,proto_main%nvtx),
c dbg
     &     zero_occ(ngastp,2)
      logical ::
     &     fix_vtx(proto_main%nvtx)

      type(contraction) ::
     &     proto

      if (ntest.ge.100) then
        write(lulog,*) '===================='
        write(lulog,*) ' gen_contr speaking'
        write(lulog,*) '===================='
        write(lulog,*) ' version: 4'
c        call prt_contr2(lulog,proto_main,op_info)
        call prt_contr3(lulog,proto_main,occ_vtx)
      end if

      if (abs(proto_main%fac).lt.1d-20) then
        call quit(1,'gen_contr4',
     &       'called with zero or near-zero prefactor')
      end if

      ! constant definitions of multi-dim. array do not work on
      ! all compilers:
      zero_occ = 0

      form_pnt => form_list

      nvtx = proto_main%nvtx

      ! get a copy of the proto-type for further processing
      call init_contr(proto)
      call copy_contr(proto_main,proto)
      call strip_contr(proto,1) ! remove the (x 0 [C]) type arcs
c dbg
c      write(lulog,*) 'the stripped proto:'
c      call prt_contr3(lulog,proto,occ_vtx)
c dbg
      do ivtx = 1, nvtx
        fix_vtx(ivtx) = fix_vtx_in(ivtx).or.ol_map(ivtx).ne.0
      end do

      ! call recursive kernel
      ivtx = 1
      call gen_contr_rec(ivtx,proto)

      ! release the allocated space for proto
      call dealloc_contr(proto)

      return

      contains

      recursive subroutine gen_contr_rec(ivtx,proto_in)

      implicit none

      integer, intent(in) ::
     &     ivtx
      type(contraction), intent(in) ::
     &     proto_in

      logical ::
     &     init, ok, must_connect(nvtx), must_not_connect(nvtx), reo
      integer ::
     &     nvtx_prev, jvtx, kvtx, iarc, narc_old, ieqvfac, ierr
      integer ::
     &     ivtx_reo(nvtx), ivtx_reo2(nvtx), ol_map_reo(nvtx),
     &     occ_ol_prev(ngastp,2), occ_ol_rem(ngastp,2),
     &     occ_ol_vtx(ngastp,2,nvtx),
     &     occ_ol_scr(ngastp,2,nvtx),
     &     occ_ex(ngastp,2), occ_dx(ngastp,2),
     &     occ_cnt2prev(ngastp,2),
     &     occ_cnt2prev_min(ngastp,2),
     &     occ_cnt2prev_max(ngastp,2),
     &     occ_cnt2rem_max(ngastp,2),
     &     occ_sum(ngastp,2)
      type(contraction) ::
     &     proto, proto_new
      logical, external ::
     &     next_part_connection, next_dist2, check_contr

      integer, allocatable ::
     &     occ_conn(:,:,:)

      if (ntest.ge.100) then
        write(lulog,*) '-------------------------------------------'
        write(lulog,*) ' recursive part, ivtx = ',ivtx
        write(lulog,*) '-------------------------------------------'
      end if

      ! get working copies of proto-contraction
      call init_contr(proto)
      call copy_contr(proto_in,proto)

      ! get currently unconnected vertex lines
      ! _vtx:  per vertex
      ! _prev: summed over 1:ivtx-1
      ! _rem : summed over ivtx+1:nvtx
      call gen_contr4_unconn(occ_ol_prev,occ_ol_rem,occ_ol_vtx,
     &     ivtx,proto,occ_vtx)
      if (ntest.ge.100) then
        write(lulog,*) 'current proto-contraction:'
        call prt_contr3(lulog,proto,occ_vtx)
        write(lulog,*) 'occ_ol_vtx:'
        call wrt_occ_n(lulog,occ_ol_vtx,nvtx)
        write(lulog,*) 'occ_ol_prev,occ_ol_rem:'
        call wrt_occ2(lulog,occ_ol_prev,occ_ol_rem)
      end if
      ! get unconnected excitation and deexcitation part of operator
      occ_ex = iocc_xdn(1,occ_ol_vtx(1:ngastp,1:2,ivtx))
      occ_dx = iocc_xdn(2,occ_ol_vtx(1:ngastp,1:2,ivtx))

      ! max. possible contractions with previous operators
      occ_cnt2prev_max = iocc_overlap(occ_ex,.true.,
     &     iocc_xdn(2,occ_ol_prev),.false.)
      ! max. possible contractions with remaining operators
      occ_cnt2rem_max = iocc_overlap(occ_dx,.false.,
     &     iocc_xdn(1,occ_ol_rem),.true.)

      if (ntest.ge.100) then
        write(lulog,*) '[EX],[DX],[C2prev],[C2rem]'
        call wrt_occ4(lulog,occ_ex,occ_dx,
     &                      occ_cnt2prev_max,occ_cnt2rem_max)
      end if

      ! validity checks
      ! (a) do not overshoot the number of excitations
      !   [EX_current] + [EX_previous] - [CNT2prev_max]
      !         element-wise <= 0
      occ_sum = occ_ex + iocc_xdn(1,occ_ol_prev)
     &     - iocc_dagger(occ_cnt2prev_max)

      ok = iocc_bound('<=',occ_sum,.false.,zero_occ,.false.)
c dbg
c      call wrt_occ2(lulog,occ_sum,occ_tgt_ex)
c      write(lulog,*) 'ok 1: ',ok
c dbg
      if (ok) then
        ! (b) do not overshoot the number of deexcitations in target
        ! [DX_current] - [CNT2rem_max] element-wise <= 0
        occ_sum = occ_dx - occ_cnt2rem_max

        ok = iocc_bound('<=',occ_sum,.false.,zero_occ,.false.)
c dbg
c        call wrt_occ(lulog,occ_sum)
c        write(lulog,*) 'ok 2: ',ok
c dbg
      end if

      if (ok) then
        ! set up array: which operators must be connected with
        ! current one?
        must_connect(1:ivtx-1) = .false.
        do iarc = 1, proto%narc
          if (proto%arc(iarc)%link(2).eq.ivtx .and.
     &         proto%arc(iarc)%occ_cnt(1,1).lt.0)
     &         must_connect(proto%arc(iarc)%link(1)) = .true.
        end do

        ! which operators must not be connect with current one?
        ! (any connections between operators for which connection
        ! was fixed on input, as we here are dealing with the
        ! additional connections, only. If a 0-contraction was
        ! given on input, any connection between these operators
        ! is avoided effectively.
        must_not_connect(1:ivtx-1) = .false.
        do iarc = 1, proto%narc
          if (proto%arc(iarc)%link(2).eq.ivtx .and.
     &        (strict.and.proto%arc(iarc)%occ_cnt(1,1).ge.0 .or.
     &    .not.strict.and.iocc_zero(proto%arc(iarc)%occ_cnt)) )
     &         must_not_connect(proto%arc(iarc)%link(1)) = .true.
        end do

        ! a must_connect overrides must_not_connect
        do jvtx = 1, ivtx-1
          must_not_connect(jvtx) = must_not_connect(jvtx).and.
     &                            .not.must_connect(jvtx)
        end do

        ! also, we should avoid contractions between open lines:
        if (ol_map(ivtx).ne.0) then
          do jvtx = 1, ivtx-1
            must_not_connect(jvtx) =
     &           must_not_connect(jvtx).or.ol_map(jvtx).ne.0
          end do
        end if

        ! count number of previous operators to which connection is
        ! possible
        ok = .true.
        nvtx_prev = 0
        do jvtx = 1, ivtx-1
          if (iocc_nonzero(occ_ol_vtx(1:ngastp,1:2,jvtx))) then
            nvtx_prev = nvtx_prev+1
            ivtx_reo(nvtx_prev) = jvtx
          else if (must_connect(jvtx)) then
            ok = .false.
          end if
        end do
c dbg
c        write(lulog,*) 'ok 3: ',ok
c        write(lulog,*) 'must_connect:     ',must_connect(1:ivtx-1)
c        write(lulog,*) 'must_not_connect: ',must_not_connect(1:ivtx-1)
c dbg

        if (ok) then
          ! (we might set up occ_cnt2prev_min here)

          ! loop over all possible connections
          ! we start with the unconnected case
          occ_cnt2prev_min = 0
          allocate(occ_conn(ngastp,2,max(1,nvtx_prev)))
          occ_cnt2prev = occ_cnt2prev_min
          cnn_dist: do
            if (ntest.ge.100) then
              write(lulog,*) 'occ_cnt2prev: (nvtx_prev = ',nvtx_prev,')'
              call wrt_occ(lulog,occ_cnt2prev)
            end if

            init = .true.
            connect: do while(next_part_connection(occ_conn,
     &           init,occ_cnt2prev,max(1,nvtx_prev),
     &           occ_cnt2prev_min,occ_cnt2prev))
              init = .false.

              if (ntest.ge.100) then
                write(lulog,*) 'current partition:'
                call wrt_occ_n(lulog,occ_conn,max(1,nvtx_prev))
              end if

              ! check connections
              ! (a) do connections fit to vertices?
              ok = .true.
              do jvtx = 1, nvtx_prev
                ok = ok.and.
     &               iocc_bound('<=',
c     &               occ_conn(1:ngastp,1:2,ivtx_reo(jvtx)),.false.,
c     &               occ_ol_vtx(1:ngastp,1:2,jvtx),.false.)
     &               occ_conn(1:ngastp,1:2,jvtx),.false.,
     &               occ_ol_vtx(1:ngastp,1:2,ivtx_reo(jvtx)),.false.)
              end do
c dbg
c              write(lulog,*) 'check (a): ',ok
c dbg

              ! (b) are all operators connected that must be connected?
              if (ok) then
                do jvtx = 1, nvtx_prev
                  if (must_connect(ivtx_reo(jvtx)))
     &               ok = ok.and.
     &               iocc_nonzero(occ_conn(1:ngastp,1:2,jvtx))
                end do
              end if
c dbg
c              write(lulog,*) 'check (b): ',ok
c dbg

              ! (c) no connections between operators for which
              !  connection was fixed (exept we generated a
              ! 0-contraction, which will be ignored below anyway)
              if (ok) then
                do jvtx = 1, nvtx_prev
                  if (must_not_connect(ivtx_reo(jvtx)))
     &                 ok = ok.and..not.
     &                     iocc_nonzero(occ_conn(1:ngastp,1:2,jvtx))
                end do
              end if
c dbg
c              write(lulog,*) 'check (c): ',ok
c dbg

              ! set up prototype
              if (ok) then
                ! create a new instance, as some of the information
                ! on proto must be retained for further "connect" loops
                call init_contr(proto_new)
                call copy_contr(proto,proto_new)

                narc_old = proto_new%narc
                do jvtx = 1, nvtx_prev
                  kvtx = ivtx_reo(jvtx)
                  ! loop over existing arcs and replace proto-arc
                  ! by actual connection
                  ok = .false.
                  do iarc = 1, narc_old
                    if (proto_new%arc(iarc)%link(1).eq.kvtx .and.
     &                  proto_new%arc(iarc)%link(2).eq.ivtx) then
c dbg
c                      write(lulog,*) '(a) add arc at ',iarc
c                      call wrt_occ(lulog,occ_conn(1:ngastp,1:2,jvtx))
c dbg
                      ! only if this was a proto-connection (indicated
                      ! by a negative number on first entry)
                      ! else: do nothing, as we have above checked that
                      ! only an additional 0-connection remains
c                      if (proto_new%arc(iarc)%occ_cnt(1,1).lt.0)
c     &                     proto_new%arc(iarc)%occ_cnt =
c     &                     occ_conn(1:ngastp,1:2,jvtx)
                      ! do not completely understand that behaviour,
                      ! trying instead:
                      if (proto_new%arc(iarc)%occ_cnt(1,1).lt.0) then
                        proto_new%arc(iarc)%occ_cnt =
     &                     occ_conn(1:ngastp,1:2,jvtx)
                      else
                        proto_new%arc(iarc)%occ_cnt =
     &                       proto_new%arc(iarc)%occ_cnt +
     &                       occ_conn(1:ngastp,1:2,jvtx)
                      end if
                      ok = .true.
                    end if
                  end do

                  ! in a few cases two or more contractions between
                  ! the same vertices end up on different arcs:
                  ! merge these and adapt narc accordingly
                  if (proto_new%narc.gt.0)
     &               call contr_clean_arcs(proto_new%arc,proto_new%narc)

                  if (.not.ok.and.
     &                 iocc_nonzero(occ_conn(1:ngastp,1:2,jvtx))) then
                    call resize_contr(proto_new,nvtx,
     &                   proto_new%narc + 1,0,0)
                    proto_new%narc = proto_new%narc + 1
                    iarc = proto_new%narc
c dbg
c                      write(lulog,*) '(b) add arc at ',iarc
c                      call wrt_occ(lulog,occ_conn(1:ngastp,1:2,jvtx))
c dbg
                    proto_new%arc(iarc)%link(1) = kvtx
                    proto_new%arc(iarc)%link(2) = ivtx
                    proto_new%arc(iarc)%occ_cnt =
     &                     occ_conn(1:ngastp,1:2,jvtx)
                  end if
                end do

                if (ivtx.lt.nvtx) then
                  ! next recursion level
                  call gen_contr_rec(ivtx+1,proto_new)!,
c     &               occ_ol_prev,occ_ol_rem,occ_ol_vtx)
                  if (ntest.ge.100) then
                    write(lulog,*) '- back in level ',ivtx
                  end if
                else
                  ! final check for validity:
                  ! no (de)excitations do remain?
                  ! [EX_current] + [EX_prev] - [CNT2prev_max^+] == 0
                  occ_sum = occ_ex + iocc_xdn(1,occ_ol_prev) -
     &                 iocc_dagger(occ_cnt2prev)
                  ok = .not.iocc_nonzero(occ_sum)
c dbg
c                  call wrt_occ(lulog,occ_sum)
c                  write(lulog,*) 'EX : ',ok
c dbg
                  if (ok) then
                    ! [DX_current] + [DX_prev] - [CNT2prev_max] == 0
                    occ_sum = occ_dx + iocc_xdn(2,occ_ol_prev) -
     &                   occ_cnt2prev
                    ok = .not.iocc_nonzero(occ_sum)
c dbg
c                    call wrt_occ(lulog,occ_sum)
c                    write(lulog,*) 'DX : ',ok
c dbg
                  end if

                  if (ok) then

                    ! make topological analysis
                    call topo_contr(ieqvfac,reo,ivtx_reo2,
     &                   proto_new,occ_vtx,fix_vtx)

                    ok = ieqvfac.gt.0
c dbg
c                    write(lulog,*) 'check ieqvfac: ',ok,ieqvfac
c dbg
                  end if

                  ! very last check:
                  ! all (0 x [C]) and (x 0 [C]) prototype arcs
                  ! must be present (as such or partitioned)
                  if (ok) then
c dbg
c                    write(lulog,*) '+----+'
c                    call prt_contr2(lulog,proto_new,op_info)
c                    call prt_contr2(lulog,proto_main,op_info)
c                    write(lulog,*) '+----+'
c dbg
                    ok = check_contr(proto_new,proto_main)

                  end if
c dbg
c                  write(lulog,*) 'final: ',ok
c dbg

                  if (ok) then

                    if (ntest.ge.100) then
                      write(lulog,*) 'new contraction found!'
                    end if

                    ! ensure canonical order
cmh                    call canon_contr(proto_new,reo,ivtx_reo2)
                    call canon_contr(proto_new,.false.,ivtx_reo2)

                    ! reorder ol_map as well
cmh                    call perm_mult(ol_map_reo,ol_map,ivtx_reo2,nvtx)

                    ! apply prefactor from equivalent permutation of vertices
                    proto_new%fac = proto_new%fac/dble(ieqvfac)

                    ! store contraction in formula list
                    ! and advance formula pointer
                    call new_formula_item(
     &                   form_pnt,command_add_contribution,
     &                   proto_new%idx_res)

                    ! transform to actual contraction form
                    ! i.e. remove external line vertices (if any)
c dbg
c                    do kvtx = 1, nvtx
c                      occ_test(1:ngastp,1:2,kvtx) =
c     &                     occ_vtx(1:ngastp,1:2,ivtx_reo2(kvtx))
c                    end do
c dbg
                    call proto2contr(form_pnt%contr,proto_new,
cmh     &                   ol_map_reo,occ_test,op_info)
     &                   ol_map,occ_test,op_info)

                    form_pnt => form_pnt%next

                  end if
                end if

                call dealloc_contr(proto_new)
              end if

            end do connect
            if (.not.next_dist2(occ_cnt2prev,2*ngastp,
     &           occ_cnt2prev_min,occ_cnt2prev_max,+1)) exit cnn_dist
          end do cnn_dist
          deallocate(occ_conn)

        end if
      end if
      call dealloc_contr(proto)

      if (ntest.ge.100) then
        write(lulog,*) 'returning from level ',ivtx
      end if

      return
      end subroutine

      subroutine gen_contr4_unconn(occ_ol_prev,occ_ol_rem,occ_ol_vtx,
     &     ivtx,proto,occ_vtx)

      implicit none

      integer, intent(in) ::
     &     ivtx,
     &     occ_vtx(ngastp,2,*)
      integer, intent(out) ::
     &     occ_ol_prev(ngastp,2), occ_ol_rem(ngastp,2),
     &     occ_ol_vtx(ngastp,2,*)
      type(contraction), intent(in) ::
     &     proto

      type(cntr_arc), pointer ::
     &     arc(:)

      integer ::
     &     idx1, idx2, jvtx, iarc

      ! initialize with vertex occupations
      do jvtx = 1, nvtx
        occ_ol_vtx(1:ngastp,1:2,jvtx) = occ_vtx(1:ngastp,1:2,jvtx)
      end do

      ! loop over all arcs and remove contractions
      arc => proto%arc
      do iarc = 1, proto%narc
        ! ignore proto-type arcs
        if (arc(iarc)%occ_cnt(1,1).lt.0) cycle
        idx1 = arc(iarc)%link(1)
        idx2 = arc(iarc)%link(2)
        ! -- dto. --
        if (idx1.le.0.or.idx2.le.0) cycle
        occ_ol_vtx(1:ngastp,1:2,idx1) = occ_ol_vtx(1:ngastp,1:2,idx1)
     &       - arc(iarc)%occ_cnt
        occ_ol_vtx(1:ngastp,1:2,idx2) = occ_ol_vtx(1:ngastp,1:2,idx2)
     &       - iocc_dagger(arc(iarc)%occ_cnt)
      end do

      if (ivtx.le.0) return

      ! accumulate all vertices up to ivtx-1:
      occ_ol_prev = 0
      do jvtx = 1, ivtx-1
        occ_ol_prev = occ_ol_prev + occ_ol_vtx(1:ngastp,1:2,jvtx)
      end do

      ! accumulate all vertices from ivtx+1 to nvtx:
      occ_ol_rem = 0
      do jvtx = ivtx+1, nvtx
        occ_ol_rem = occ_ol_rem + occ_ol_vtx(1:ngastp,1:2,jvtx)
      end do

      return
      end subroutine

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
20070710 04:22:24 nbsp hey sol, nice to see you on the daviswiki and updating the Technocultural Studies page. Users/JessicaRockwell

20071025 05:01:54 nbsp Sol, hope you are BothWell and Wonder Full Users/AyseG

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine thom (u,v,w,mySeed)
c***********************************************************************
c     selects a scattering angle for a thomson collision and then
c     generates new direction cosines.
c***********************************************************************
      implicit double precision (a-h,o-z)
      integer mySeed(*)
      real*8 ranf
      external ranf
      dimension csn (128)

      data(csn(i),i=1,67)/
     |     1.       d0,.98952871d0,.9789466 d0,.96825135d0,.95744053d0,
     x .94651168d0,
     |     .93546224d0,.92428964d0,.91299119d0,.90156417d0,.89000577d0,
     x .87831314d0,
     |     .86648335d0,.8545134d0, .84240023d0,.83014072d0,.81773168d0,
     x .80516986d0,
     |     .79245197d0,.77957465d0,.76653448d0,.75332801d0,.73995175d0,
     x .72640218d0,
     |     .71267574d0,.69876886d0,.68467798d0,.67039954d0,.65592999d0,
     x .64126583d0,
     |     .62640361d0,.61133997d0,.59607164d0,.58059548d0,.5649085 d0,
     x .5490079d0,
     |     .53289112d0,.51655582d0,.5d0,       .48322199d0,.46622053d0,
     x .44899477d0,
     |     .43154441d0,.41386968d0,.39597143d0,.3778512 d0,.35951125d0,
     x .34095466d0,
     |     .32218536d0,.30320817d0,.2840289 d0,.26465437d0,.24509241d0,
     x .22535195d0,
     |     .20544297d0,.18537655d0,.1651648 d0,.14482089d0,.12435893d0,
     x .10379394d0,
     |     .083141761d0,.062418936d0,.041642596d0,.020830322d0,.0d0,
     |     -.020830322d0,-.041642596d0/

      data(csn(i),i=68,128)/
     |     -.062418936d0,-.083141761d0,-.10379394d0,-.12435893d0,
     x -.14482089d0,
     |     -.1651648d0,-.18537655d0,-.20544297d0,-.22535195d0,
     x -.24509241d0,
     |     -.26465437d0,-.2840289d0,-.30320817d0,-.32218536d0,
     x -.34095466d0,-.35951125d0,
     |     -.3778512d0,-.39597143d0,-.41386968d0,-.43154441d0,
     x -.44899477d0,-.46622053d0,
     |     -.48322199d0,-.5d0,-.51655582d0,-.53289112d0,-.5490079d0,
     x -.5649085d0,
     |     -.58059548d0,-.59607164d0,-.61133997d0,-.62640361d0,
     x -.64126583d0,-.65592999d0,
     |     -.67039954d0,-.68467798d0,-.69876886d0,-.71267574d0,
     x -.72640218d0,-.73995175d0,
     |     -.75332801d0,-.76653448d0,-.77957465d0,-.79245197d0,
     x -.80516986d0,-.81773168d0,
     |     -.83014072d0,-.84240023d0,-.8545134d0,-.86648335d0,
     x -.87831314d0,
     |     -.89000577d0,-.90156417d0,-.91299119d0,
     |     -.92428964d0,-.93546224d0,
     |     -.94651168d0,-.95744053d0,-.96825135d0,-.9789466d0,
     x -.98952871d0/

      index  = ranf(mySeed)*128 + 1
      sctang = csn (index)
      t3     = 6.2831853d0 * ranf(mySeed)
      t4     = sqrt(1.0d0 - sctang*sctang)
      t5     = sqrt(1.0d0 - w*w)

c--------------------update scattered photon direction
c--------------------(from tid-26607)
      temp1 = t4/t5
      temp2 = sin (t3)
      temp3 = cos (t3)
      temp4 = w*temp3
      t1    = u

c--------------------new u, v, w

      u = u*sctang + temp1*(u*temp4 - v*temp2)
      v = v*sctang + temp1*(v*temp4 + t1*temp2)
      w = w*sctang - (t4*t5)*temp3

      return
      end

max_stars_repo_name: Nuullll/llvm-test-suite

max_issues_repo_name: Nuullll/llvm-test-suite

max_forks_repo_name: Nuullll/llvm-test-suite


text: http://www.ice.ucdavis.edu/seed/ Students for the Earth and Environment or SEED was founded in 1994 and seeks to provide environmental education outreach in Davis elementary schools.


text:
subroutine setup_outfile (ncid, outfil, tag, dateid, datesecid, &
                          timeid, icesst, nlon, nlat, xlon,     &
                          xlat, iyr1out, mon1out, history)
   use prec
   use types

   implicit none

   include 'netcdf.inc'

   integer, intent(out) :: ncid                   ! output file netcdf id
   character(len=*), intent(in) :: outfil         ! output file name
   character(len=*), intent(in) :: tag            ! 'CLIM' or 'AMIP'
   integer, intent(out) :: dateid, datesecid      ! date variable ids
   integer, intent(out) :: timeid                 ! time variable id
   type(icesstparms), intent(inout) :: icesst(2)  ! struct containing ice or sst-specific info
   integer, intent(in) :: nlon                    ! number of longitudes
   integer, intent(in) :: nlat                    ! number of latitudes
   real(r8), intent(in) :: xlon(nlon)             ! longitude coordinate values
   real(r8), intent(in) :: xlat(nlat)             ! latitude coordinate values
   integer, intent(in) :: iyr1out                 ! start year written to output file
   integer, intent(in) :: mon1out                 ! start month written to output file
   character(len=*), intent(in) :: history        ! history attribute

   integer :: ret                                 ! return code
   integer :: dimids(3)                           ! dimension ids for multi-dimensioned arrays
   integer :: lonid, latid                        ! lon, lat variable ids

   character(len=32) :: name                      ! variable name
   character(len=32) :: str                       ! attribute string
!
! Open file and define dimensions
!
   call wrap_nf_create (outfil, NF_CLOBBER, ncid)
   call wrap_nf_def_dim (ncid, 'lon', nlon, dimids(1))
   call wrap_nf_def_dim (ncid, 'lat', nlat, dimids(2))
   call wrap_nf_def_dim (ncid, 'time', NF_UNLIMITED, dimids(3))
!
! Model-readable date info
!
   call wrap_nf_def_var (ncid, 'date', NF_INT, 1, dimids(3), dateid)
   call wrap_nf_put_att_text (ncid, dateid, 'long_name', 'current date (YYYYMMDD)')

   call wrap_nf_def_var (ncid, 'datesec', NF_INT, 1, dimids(3), datesecid)
   call wrap_nf_put_att_text (ncid, datesecid, 'long_name', 'current seconds of current date')
!
! Dimension variables
!
   call wrap_nf_def_var (ncid, 'lon', NF_DOUBLE, 1, dimids(1), lonid)
   call wrap_nf_put_att_text (ncid, lonid, 'long_name', 'longitude')
   call wrap_nf_put_att_text (ncid, lonid, 'units', 'degrees_east')

   call wrap_nf_def_var (ncid, 'lat', NF_DOUBLE, 1, dimids(2), latid)
   call wrap_nf_put_att_text (ncid, latid, 'long_name', 'latitude')
   call wrap_nf_put_att_text (ncid, latid, 'units', 'degrees_north')

   call wrap_nf_def_var (ncid, 'time', NF_DOUBLE, 1, dimids(3), timeid)
   write(str,100) iyr1out, mon1out
100 format('days since ',i4.4,'-',i2.2,'-01 00:00:00')
   call wrap_nf_put_att_text (ncid, timeid, 'units', trim(str))
   call wrap_nf_put_att_text (ncid, timeid, 'calendar', '365_day')

   ! history attribute
   call wrap_nf_put_att_text (ncid, NF_GLOBAL, 'history', trim(history))

!
! Define variables
!
   select case (tag)
   case ('CLIM')
!
! ice
!
      name = icesst(1)%climname
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(1)%varidclim)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidclim, 'long_name', icesst(1)%climlongname)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidclim, 'units', icesst(1)%units)

      name = icesst(1)%climname_pre
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(1)%varidclim_pre)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidclim_pre, 'long_name', icesst(1)%climlongname_pre)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidclim_pre, 'units', icesst(1)%units)
!
! sst
!
      name = icesst(2)%climname
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(2)%varidclim)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidclim, 'long_name', icesst(2)%climlongname)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidclim, 'units', icesst(2)%units)

      name = icesst(2)%climname_pre
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(2)%varidclim_pre)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidclim_pre, 'long_name', icesst(2)%climlongname_pre)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidclim_pre, 'units', icesst(2)%units)
   case ('AMIP')
!
! ice
!
      name = icesst(1)%amipname
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(1)%varidamip)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidamip, 'long_name', icesst(1)%amiplongname)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidamip, 'units', icesst(1)%units)

      name = icesst(1)%amipname_pre
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(1)%varidamip_pre)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidamip_pre, 'long_name', icesst(1)%amiplongname_pre)
      call wrap_nf_put_att_text (ncid, icesst(1)%varidamip_pre, 'units', icesst(1)%units)
!
! sst
!
      name = icesst(2)%amipname
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(2)%varidamip)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidamip, 'long_name', icesst(2)%amiplongname)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidamip, 'units', icesst(2)%units)

      name = icesst(2)%amipname_pre
      call wrap_nf_def_var (ncid, name, NF_FLOAT, 3, dimids, icesst(2)%varidamip_pre)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidamip_pre, 'long_name', icesst(2)%amiplongname_pre)
      call wrap_nf_put_att_text (ncid, icesst(2)%varidamip_pre, 'units', icesst(2)%units)
   case default
      write(6,*) 'tag ', tag, ' is an invalid tag: valid values are CLIM and AMIP'
      stop 999
   end select

   ret = nf_enddef (ncid)
   if (ret < 0) then
      call handle_error (ret)
   end if

   call wrap_nf_put_var_double (ncid, lonid, xlon)
   call wrap_nf_put_var_double (ncid, latid, xlat)

   return
end subroutine setup_outfile

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
C    @(#)kmp_wydt.f	20.5 2/28/00
        integer function kmp_wydt (i,j)

        include 'ipfinc/parametr.inc'
        include 'ipfinc/bus.inc'
        include 'ipfinc/arcntl.inc'
        include 'ipfinc/qsdup.inc'
        include 'ipfinc/qksrt.inc'

        integer MAXWYEDELTA
        parameter (MAXWYEDELTA = 500)
        common /wye_delta/ num_wye_delta, lwye_delta(3,MAXWYEDELTA),
     &                     wye_delta(MAXWYEDELTA)

        integer MAXPTIRECORDS
        parameter (MAXPTIRECORDS = 16000)
        common /scratch/ count, array(4,MAXBUS), htable_2(MAXBUS),
     &                   nextptr_2(MAXBUS), count_newbus,
     &                   newbusno(MAXBUS), count_newzone,
     &                   newzoneno(MAXCZN), count_newown,
     &                   newownno(MAXOWN),  tempc(MAXPTIRECORDS),
     &                   sort_tempc(MAXPTIRECORDS), txflag(MAXBUS),
     &                   txindex(2,MAXWYEDELTA), procnode(MAXBUS),
     &                   vangle(MAXBUS), itemp(MAXBUS)
        integer array, count, htable_2, count_newbus, newbusno,
     &          count_newzone, newzoneno, count_newown, newownno,
     &          sort_tempc, txflag, txindex, procnode
        character tempc*80

        character id1*1, id2*1
c
c       Sort fields:
c       1. base1 (High to Low)
c       2. base2 (High to Low)
c       3. bus1
c       4. bus2
c       5. id
c
        k1 = lwye_delta(1,i)
        m1 = lwye_delta(1,j)
        kmp_wydt = -100.0 * (base(k1) - base(m1))
        if (kmp_wydt .eq. 0) then
          k2 = lwye_delta(2,i)
          m2 = lwye_delta(2,j)
          kmp_wydt = -100.0 * (base(k2) - base(m2))
          if (kmp_wydt .eq. 0) then
            kmp_wydt = kompr (bus(k1), bus(m1), junk)
            if (kmp_wydt .eq. 0) then
              kmp_wydt = kompr (bus(k2), bus(m2), junk)
              if (kmp_wydt .eq. 0) then
                id1 = char (lwye_delta(3,i))
                id2 = char (lwye_delta(3,j))
                kmp_wydt = kompr (id1, id2, junk)
              endif
            endif
          endif
        endif
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
c
c -----------------------------------------------------
c
      subroutine valout (lst, lend, time, nvar, naux)
c
c     # Output the results for a general system of conservation laws
c     # in 3 dimensions
c
c     # Write the results to the file fort.q<iframe>
c     # Use format required by matlab script  plotclaw3.m
c
c     # set outaux = .true. to also output the aux arrays to fort.a<iframe>

      use amr_module
      implicit double precision (a-h,o-z)
      character*10  fname1, fname2, fname3, fname4, fname5

      logical outaux
      integer output_aux_num

      !for timing
      integer :: clock_start, clock_finish, clock_rate
      real(kind=8) cpu_start, cpu_finish

      iadd(ivar,i,j,k)   = loc     +    (ivar-1)
     &                             +    (i-1)*nvar
     &                             +    (j-1)*nvar*mitot
     &                             +    (k-1)*nvar*mitot*mjtot
      iaddaux(iaux,i,j,k) = locaux +    (iaux-1)
     &                             +    (i-1)*naux
     &                             +    (j-1)*naux*mitot
     &                             +    (k-1)*naux*mitot*mjtot

      !timing
      call system_clock(clock_start,clock_rate)
      call cpu_time(cpu_start)

      output_aux_num = 0
      do i = 1, naux
        output_aux_num = output_aux_num + output_aux_components(i)
      end do

c     # Currently outputs all aux components if any are requested!
      outaux = ((output_aux_num > 0) .and.
     .         ((.not. output_aux_onlyonce) .or. (time==t0)))

c     open(unit=77,file='fort.b',status='unknown',access='stream')

c     ### Python graphics output
c
c        ###  make the file names and open output files
      fname1 = 'fort.qxxxx'
      fname2 = 'fort.txxxx'
      fname3 = 'fort.axxxx'
      fname4 = 'fort.bxxxx'
      matunit1 = 50
      matunit2 = 60
      matunit3 = 70
      matunit4 = 71
      nstp     = matlabu
      do 55 ipos = 10, 7, -1
         idigit = mod(nstp,10)
         fname1(ipos:ipos) = char(ichar('0') + idigit)
         fname2(ipos:ipos) = char(ichar('0') + idigit)
         fname3(ipos:ipos) = char(ichar('0') + idigit)
         fname4(ipos:ipos) = char(ichar('0') + idigit)
         nstp = nstp / 10
 55   continue

      open(unit=matunit1,file=fname1,status='unknown',
     .       form='formatted')

      if (output_format == 3) then
c         # binary output
          open(unit=matunit4,file=fname4,status='unknown',
     &            access='stream')
          endif

      level = lst
      ngrids = 0
c65   if (level .gt. lfine) go to 90
 65   if (level .gt. lend) go to 90
         mptr = lstart(level)
 70      if (mptr .eq. 0) go to 80
           ngrids  = ngrids + 1
           nx      = node(ndihi,mptr) - node(ndilo,mptr) + 1
           ny      = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
           nz      = node(ndkhi,mptr) - node(ndklo,mptr) + 1
           loc     = node(store1, mptr)
           locaux  = node(storeaux,mptr)
           mitot   = nx + 2*nghost
           mjtot   = ny + 2*nghost
           mktot   = nz + 2*nghost
           write(matunit1,1001) mptr, level, nx, ny, nz
 1001 format(i5,'                 grid_number',/,
     &       i5,'                 AMR_level',/,
     &       i5,'                 mx',/,
     &       i5,'                 my',/,
     &       i5,'                 mz')

           xlow = rnode(cornxlo,mptr)
           ylow = rnode(cornylo,mptr)
           zlow = rnode(cornzlo,mptr)
           write(matunit1,1002) xlow,ylow,zlow,hxposs(level),
     &           hyposs(level),hzposs(level)
 1002 format(e18.8,'    xlow', /,
     &       e18.8,'    ylow', /,
     &       e18.8,'    zlow', /,
     &       e18.8,'    dx', /,
     &       e18.8,'    dy', /,
     &       e18.8,'    dz',/)

         if (output_format == 1) then
            do 75 k = nghost+1, mktot-nghost
            do 76 j = nghost+1, mjtot-nghost
            do 77 i = nghost+1, mitot-nghost
            do ivar=1,nvar
               if (dabs(alloc(iadd(ivar,i,j,k))) < 1d-90) then
                  alloc(iadd(ivar,i,j,k)) = 0.d0
               endif
            enddo
            write(matunit1,109)
     &            (alloc(iadd(ivar,i,j,k)), ivar=1,nvar)
  109       format(50e26.16)
   77       continue
            write(matunit1,*) ' '
   76       continue
            write(matunit1,*) ' '
   75       continue
            write(matunit1,*) ' '
         endif

         if (output_format == 3) then
c            # binary output
             i1 = iadd(1,1,1,1)
             i2 = iadd(nvar,mitot,mjtot,mktot)
c            # NOTE: we are writing out ghost cell data also, unlike ascii
             write(matunit4) alloc(i1:i2)
         endif

         mptr = node(levelptr, mptr)
         go to 70
   80 level = level + 1
      go to 65

   90 continue

c       -------------------
c       # output aux arrays
c       -------------------

      if (outaux) then
c        #  output aux array to fort.aXXXX

         level = lst
         ngrids = 0
 165     if (level .gt. lend) go to 190
            mptr = lstart(level)
 170        if (mptr .eq. 0) go to 180
              ngrids  = ngrids + 1
              nx      = node(ndihi,mptr) - node(ndilo,mptr) + 1
              ny      = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
              nz      = node(ndkhi,mptr) - node(ndklo,mptr) + 1
              loc     = node(store1, mptr)
              locaux  = node(storeaux,mptr)
              mitot   = nx + 2*nghost
              mjtot   = ny + 2*nghost
              mktot   = nz + 2*nghost

             if (output_format == 1) then
                open(unit=matunit3,file=fname3,status='unknown',
     .               form='formatted')
                write(matunit3,1001) mptr,level,nx,ny,nz
                xlow = rnode(cornxlo,mptr)
                ylow = rnode(cornylo,mptr)
                zlow = rnode(cornzlo,mptr)
                write(matunit3,1002) xlow,ylow,zlow,hxposs(level),
     &                               hyposs(level),hzposs(level)

                do 175 k = nghost+1, mktot-nghost
                do 176 j = nghost+1, mjtot-nghost
                do 177 i = nghost+1, mitot-nghost
                do ivar=1,naux
                   if (dabs(alloc(iaddaux(ivar,i,j,k))) .lt. 1d-90) then
                            alloc(iaddaux(ivar,i,j,k)) = 0.d0
                   endif
                enddo
                write(matunit3,109)
     &               (alloc(iaddaux(ivar,i,j,k)), ivar=1,naux)
 177            continue
                write(matunit3,*) ' '
 176            continue
                write(matunit3,*) ' '
 175            continue
                write(matunit3,*) ' '
             endif

             if (output_format == 3) then
c            # binary output
                open(unit=matunit3,file=fname3,status='unknown',
     &               access='stream')
                i1 = iaddaux(1,1,1,1)
                i2 = iaddaux(naux,mitot,mjtot,mktot)
c               # NOTE: we are writing out ghost cell data also, unlike ascii
                write(matunit3) alloc(i1:i2)
             endif

            mptr = node(levelptr, mptr)
            go to 170
  180 level = level + 1
      go to 165

  190 continue

      close(unit=matunit3)
      endif  !# end outputting aux array

c     --------------
c     # fort.t file:
c     --------------

      open(unit=matunit2,file=fname2,status='unknown',
     .       form='formatted')

      ndim = 3
c     # NOTE: we need to print out nghost too in order to strip
c     #       ghost cells from q when reading in pyclaw.io.binary
      write(matunit2,1000) time,nvar,ngrids,naux,ndim,nghost
 1000 format(e18.8,'    time', /,
     &       i5,'                 meqn'/,
     &       i5,'                 ngrids'/,
     &       i5,'                 naux'/,
     &       i5,'                 ndim'/,
     &       i5,'                 nghost'/,/)
c

      write(6,601) matlabu,time
  601 format('AMRCLAW: Frame ',i4,
     &       ' output files done at time t = ', d13.6,/)

      matlabu = matlabu + 1

      close(unit=matunit1)
      close(unit=matunit2)
      if (output_format == 3) then
          close(unit=matunit4)
          endif

      call system_clock(clock_finish,clock_rate)
      call cpu_time(cpu_finish)
      timeValout = timeValout + clock_finish - clock_start
      timeValoutCPU = timeValoutCPU + cpu_finish - cpu_start

      return
      end

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


!==========================================================================================
! SUBRUTINA DE CALCULO DE PROPIEDADES DE AUTO DIFUSION. DESPLAZAMIENTO CUADRATICO MEDIO Y
! COEFICIENTE DE DIFUSION DEPENDIENTE DEL TIEMPO
!
! AUTOR: MARTIN ALEJANDRO PAREDES SOSA
!==========================================================================================

Subroutine WDT(k)

  Use cte
  Implicit None

  Integer :: k, i, j, p, nmax                              !CONTADORES
  Real :: Ti, Time                                         !CONTROL DE TIEMPO
  Real :: Wtx, Wty, Wtz, Wt                                !DESPLAZAMIENTO CUADRATICO MEDIO
  Real :: Dif                                              !DIFUSION

  Open(96, File="wdt.dat")

  !TIEMPO ENTRE CONFIGURACIONES
  Ti = Real(iSave2) * dt

  !BARRIDO TEMPORAL
  TEMPO: Do i = 1 , k-1

     nmax = k - i

     Wtx = 0.0
     Wty = 0.0
     Wtz = 0.0
     Wt = 0.0

     !BARRIMIENTO DE PARTICULAS
     BPart: Do p = 1, N

        !BARRIDO EN TIEMPO
        Tiempo: Do j = 1, nmax

           Wtx = Wtx + ( CXD(p,i+j) - CXD(p,j)  ) * ( CXD(p,i+j) - CXD(p,j)  )
           wty = Wty + ( CYD(p,i+j) - CYD(p,j)  ) * ( CYD(p,i+j) - CYD(p,j)  )
           Wtz = Wtz + ( CZD(p,i+j) - CZD(p,j)  ) * ( CZD(p,i+j) - CZD(p,j)  )

        End Do Tiempo

     End Do BPart

     Time = Ti * Real(i)
     Wt = (Wtx + Wty + Wtz ) / ( Real(nmax) * Real(N) * 6.0 )
     Dif = Wt / Time

     Write(96,*) Time, Wt, Dif

  End Do TEMPO
  Close(96)

End Subroutine WDT