text:
! Subroutine ETRANS

! 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.

! Banks & Nagy 2-stream electron transport code
! Adapted by Stan Solomon, 1986, 1988
! Uses variable altitude and energy grids
! Updated comments and removed artifacts, scs, 2005
! Moved common blocks into use-associated variables defined in cglow.f90, btf, 2015
! Changed input to subroutine impit to argument list, scs, 2015
! Modernized to remove upper case and numbered line statements, scs, 2016
! Refactored for f90, scs, 2016

! Subroutine EXSECT called first time only to calculate electron impact cross sections.

! Definitions:
! Use-associated variables (formerly in common blocks):  See glow.f90, cglow.f90, and exsect.f
! psi    first term of parabolic d.e., = 1
! alpha  second term "; cm-1
! beta   third term  "; cm-2
! gamma  forth term  "; cm-4 s-1 ev-1
! delz   altitude increments; cm
! del2   sum of altitude increment and next higher increment; cm
! dela   average of "
! delp   product of dela and next higher delz
! delm   product of dela and delz
! dels   product of delz and next higer delz
! den    dummy array for transfer of calculated downward flux
! fac    factor for extrapolating production rate, = 0
! prod   sum of photoelectron production and secondary electrons from protons; cm-3 s-1 ev-1
! eprod  energy of "; ev cm-3
! t1     elastic collision term; cm-1
! t2     elastic + inelastic collision term; cm-1
! tsa    total energy loss cross section for each species; cm2
! produp upward cascade + secondary production; cm-3 s-1 ev-1
! prodwn downward "
! phiup  upward flux; cm-2 s-1 ev-1
! phidwn downward "
! tsigne thermal electron collision term; cm-1
! secion total ionization rate; cm-3 s-1
! secp   secondary electron production; cm-3 s-1 ev-1
! r1     ratio term for calculating upward flux; cm-2 s-1 ev-1
! expt2  exponential term for calculating upward flux
! produa collection array for calculating produp; cm-3 s-1 ev-1
! prodda  "                               prodwn
! phiinf downward flux at top of atmos., divided by avmu; cm-2 s-1 ev-1
! potion ionizaition potential for each species; ev
! avmu   cosine of the average pitch angle

! Array dimensions:
! jmax    number of altitude levels
! nbins   number of energetic electron energy bins
! nmaj    number of major species
! nei     number of states produced by electron impact

  subroutine etrans

    use cglow,only: nmaj,nbins,jmax,nei,ierr,jlocal, &
                    dip,ener,del,aglw,eheat,sion,phitop,zz,pespec, &
                    sespec,zte,ze,zmaj,uflx,dflx,tez,efrac           ! formerly /cglow/
    use cglow,only: siga,sigs,pe,sigex,sec,iimaxx,pin                ! formerly /cxsect/
    use cglow,only: ww                                               ! formerly /cxpars/

    implicit none

    integer,save :: ifirst=1
    integer :: ii,ib,ibb,i,n,jj,j,k,jjj4,iv,ll,kk,im,iq
    real :: prod(jmax), eprod(jmax), t1(jmax), t2(jmax), tsa(nmaj), &
            produp(jmax,nbins), prodwn(jmax,nbins), &
            phiup(jmax), phidwn(jmax), tsigne(jmax), taue(jmax), &
            secion(jmax), secp(nmaj,jmax), r1(jmax), expt2(jmax), &
            produa(jmax), prodda(jmax), phiinf(nbins), potion(nmaj), &
            alpha(jmax),beta(jmax),gamma(jmax),psi(jmax),del2(jmax), &
            delp(jmax),delm(jmax),dels(jmax),den(jmax), &
            delz(jmax),dela(jmax)
    real :: sindip,rmusin,phiout,dag,et,eet,fluxj,edep,epe,ephi,aprod,ein,eout,fac
    real,parameter :: avmu=0.5

    potion = (/16.,16.,18./)
    ierr = 0
    fac = 0.
    sindip = sin(dip)
    rmusin = 1. / sindip / avmu
    psi(1)   = 1.
!
! First call only:  calculate cross-sectons:
!
    if (ifirst == 1) then
      call exsect (ener, del)
      ifirst = 0
    endif
!
! Zero variables:
!
    alpha(1) = 0.
    beta(1) = 0.
    gamma(1) = 0.
    phiout = 0.0
    eheat(:) = 0.0
    eprod(:) = 0.0
    secion(:) = 0.0
    sion(:,:) = 0.0
    aglw(:,:,:) = 0.0
    produp(:,:) = 1.0e-20
    prodwn(:,:) = 1.0e-20
!
! Divide downward flux at top of atmos. by average pitch angle cosine:
!
    phiinf(:) = phitop(:) / avmu
!
! Calcualte delta z's:
!
    delz(1) = zz(2)-zz(1)
    do i=2,jmax
      delz(i) = zz(i)-zz(i-1)
    enddo
    do i=1,jmax-1
      del2(i) = delz(i)+delz(i+1)
      dela(i) = del2(i)/2.
      delp(i) = dela(i)*delz(i+1)
      delm(i) = dela(i)*delz(i)
      dels(i) = delz(i)*delz(i+1)
    enddo
    del2(jmax) = del2(jmax-1)
    dela(jmax) = dela(jmax-1)
    delp(jmax) = delp(jmax-1)
    delm(jmax) = delp(jmax-1)
    dels(jmax) = dels(jmax-1)
!
! Top of energy loop:
!
    do j=nbins,1,-1
!
! Calculate production:
!
      do i = 1, jmax
        prod(i) = (pespec(j,i)+sespec(j,i)) * rmusin / del(j)
        eprod(i) = eprod(i) + prod(i) * ener(j) * del(j) / rmusin
      enddo
!
! Total energy loss cross section for each species:
!
      tsa(:) = 0.0
      if (j > 1) then
        do k = 1, j-1
          do i = 1, nmaj
            tsa(i) = tsa(i) + siga(i,k,j) * (del(j-k)/del(j))
          enddo
        enddo
      else
        do i=1,nmaj
          tsa(i) = tsa(i) + siga(i,1,j) + 1.e-18
        enddo
      endif
!
! Thermal electron energy loss:
!
      jjj4 = j - 1
      if (j == 1) jjj4 = 1
      dag = ener(j) - ener(jjj4)
      if (dag <= 0.0) dag = del(1)
!
      do i = 1, jmax
        et = 8.618e-5 * zte(i)
        eet = ener(j) - et
        if (eet <= 0.0) then
          tsigne(i) = 0.0
        else
          tsigne(i) = ((3.37e-12*ze(i)**0.97)/(ener(j)**0.94)) &
                    * ((eet)/(ener(j) - (0.53*et))) ** 2.36
        endif
        tsigne(i) = tsigne(i) * rmusin / dag
      enddo
!
! Collision terms:
!
      do i = 1, jmax
        t1(i) = 0.0
        t2(i) = 0.0
        do iv = 1, nmaj
          t1(i) = t1(i) + zmaj(iv,i) * sigs(iv,j) * pe(iv,j)
          t2(i) = t2(i) + zmaj(iv,i) * (sigs(iv,j)*pe(iv,j) + tsa(iv))
        enddo
        t1(i) = t1(i) * rmusin
        t2(i) = t2(i) * rmusin + tsigne(i)
      enddo
!
! Bypass next section if local calculation was specified:
!
      if (jlocal /= 1) then
!
! Solve parabolic d.e. by Crank-Nicholson method to find downward flux:
!
        do i = 2, jmax-1
          psi(i) = 1.
          alpha(i) = (t1(i-1) - t1(i+1)) / (del2(i) * t1(i))
          beta(i) = t2(i) * (t1(i+1) - t1(i-1)) / (t1(i) * del2(i)) &
                  - (t2(i+1) - t2(i-1)) / del2(i) - t2(i)**2 + t1(i)**2
          if (prod(i) < 1.e-30) prod(i) = 1.e-30
          if (prodwn(i,j) < 1.e-30) prodwn(i,j) = 1.e-30
          gamma(i) = (prod(i)/2.0) * (-t1(i) - t2(i) - alpha(i) &
                   - (prod(i+1) - prod(i-1))/prod(i)/del2(i)) &
                   + prodwn(i,j) * (-alpha(i) - t2(i) &
                   - (prodwn(i+1,j)-prodwn(i-1,j))/prodwn(i,j)/del2(i)) &
                   - produp(i,j) * t1(i)
        enddo
        if (abs(beta(2)) < 1.e-20) then
          beta(2) = 1.e-20
          ierr = 2
        endif
        phidwn(2) = gamma(2) / beta(2)
        den(1) = phidwn(2)
        fluxj = phiinf(j)
        call impit(jmax,fluxj,fac,alpha,beta,gamma,psi,del2,delp,delm,dels,den)
        phidwn(:) = den(:)
!
! Apply lower boundary condition: phiup=phidwn.  Should be nearly zero.
!
        phiup(1) = phidwn(1)
!
! Integrate back upward to calculate upward flux:
!
        do i = 2, jmax
          r1(i) = (t1(i)*phidwn(i) + (prod(i)+2.*produp(i,j))/2.) / t2(i)
          taue(i) = t2(i)*delz(i)
          if (taue(i) > 60.) taue(i)=60.
          expt2(i) = exp(-taue(i))
        enddo
        do i=2,jmax
          phiup(i) = r1(i) + (phiup(i-1)-r1(i)) * expt2(i)
        enddo

      else
!
! Local calculation only:
!
        do i = 1, jmax
          if (t2(i) <= t1(i)) then
            ierr = 1
            t2(i) = t1(i) * 1.0001
          endif
          phiup(i) = (prod(i)/2.0 + produp(i,j)) / (t2(i) - t1(i))
          phidwn(i) = (prod(i)/2.0 + prodwn(i,j)) / (t2(i) - t1(i))
        enddo

      endif
!
! Multiply fluxes by average pitch angle cosine and put in arrays:
!
      do i=1,jmax
        uflx(j,i) = phiup(i) * avmu
        dflx(j,i) = phidwn(i) * avmu
      enddo
!
!  Calculate outgoing electron energy flux for conservation check:
!
      phiout = phiout + phiup(jmax) * del(j) * ener(j)
!
! Cascade production:
!
      if (j > 1) then
        do k = 1, j-1
          ll = j - k
          produa(:)=0.
          prodda(:)=0.
          do n = 1, nmaj
            do i=1,jmax
              produa(i) = produa(i) &
                         + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phidwn(i) &
                         + (1. - pin(n,j))*siga(n,k,j)*phiup(i))
              prodda(i) = prodda(i) &
                         + zmaj(n,i) * (siga(n,k,j)*pin(n,j)*phiup(i) &
                         + (1. - pin(n,j))*siga(n,k,j)*phidwn(i))
            enddo
          enddo
          do i=1,jmax
            produp(i,ll) = produp(i,ll) + produa(i) * rmusin
            prodwn(i,ll) = prodwn(i,ll) + prodda(i) * rmusin
          enddo
        enddo
      endif
      kk = j - 1
      if (kk > 0) then
        do i = 1, jmax
          produp(i,kk) = produp(i,kk)+tsigne(i)*phiup(i)*(del(j)/del(kk))
          prodwn(i,kk) = prodwn(i,kk)+tsigne(i)*phidwn(i)*(del(j)/del(kk))
        enddo
      endif
!
! Electron heating rate:
!
      dag = del(j)
      do i = 1, jmax
        eheat(i) = eheat(i) + tsigne(i) * (phiup(i)+phidwn(i)) * dag**2
      enddo
!
! Electron impact excitation rates:
!
      do ii = 1, jmax
        do i = 1, nmaj
          do ibb = 1, nei
            aglw(ibb,i,ii) = aglw(ibb,i,ii) + (phiup(ii) + phidwn(ii)) &
                             * sigex(ibb,i,j) * del(j) * zmaj(i,ii)
          enddo
        enddo
      enddo
!
! Calculate production of secondaries into k bin for energy j bin and add to production:
!
      do k = 1, iimaxx(j)
        do n = 1, nmaj
          do i = 1, jmax
            secp(n,i) = sec(n,k,j) * zmaj(n,i) * (phiup(i) + phidwn(i))
            sion(n,i) = sion(n,i) + secp(n,i) * del(k)
            secion(i) = secion(i) + secp(n,i) * del(k)
            produp(i,k) = produp(i,k) + (secp(n,i)*.5*rmusin)
            prodwn(i,k) = prodwn(i,k) + (secp(n,i)*.5*rmusin)
          enddo
        enddo
      enddo

    enddo       ! bottom of energy loop

      eheat(:) = eheat(:) / rmusin
!
! Calculate energy deposited as a function of altitude and total energy deposition:
!
    edep = 0.
    do im=1,jmax
      tez(im) = eheat(im)
      do ii=1,nmaj
        tez(im) = tez(im) + sion(ii,im)*potion(ii)
        do iq=1,nei
          tez(im) = tez(im) + aglw(iq,ii,im)*ww(iq,ii)
        enddo
      enddo
      edep = edep + tez(im) * dela(im)
    enddo
!
! Calculate energy input, output, and fractional conservation:
!
    epe = 0.0
    ephi = 0.0
    do i = 2, jmax
      aprod = sqrt(eprod(i)*eprod(i-1))
      epe = epe + aprod * delz(i)
    enddo
    do jj = 1, nbins
      ephi = ephi + phiinf(jj) * ener(jj) * del(jj) / rmusin
    enddo
    ein = ephi + epe
    phiout = phiout / rmusin
    eout = edep + phiout
    efrac = (eout - ein) / ein

    return

  end subroutine etrans

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

! Subroutine impit solves parabolic differential equation by implicit Crank-Nicholson method

  subroutine impit(jmax,fluxj,fac,alpha,beta,gamma,psi,del2,delp,delm,dels,den)

    implicit none

    integer,intent(in) :: jmax
    real,intent(in) :: fluxj,fac,alpha(jmax),beta(jmax),gamma(jmax), &
               psi(jmax),del2(jmax),delp(jmax),delm(jmax),dels(jmax)
    real,intent(out) :: den(jmax)

    integer :: i1,i,kk,jk
    real :: dem,k(jmax),l(jmax),a(jmax),b(jmax),c(jmax),d(jmax)

    i1 = jmax - 1
    do i = 1, i1
      a(i) = psi(i) / delp(i) + alpha(i) / del2(i)
      b(i) = -2. * psi(i) / dels(i) + beta(i)
      c(i) = psi(i) / delm(i) - alpha(i) / del2(i)
      d(i) = gamma(i)
    enddo
    k(2) = (d(2) - c(2)*den(1)) / b(2)
    l(2) = a(2) / b(2)
    do i = 3, i1
      dem = b(i) - c(i) * l(i-1)
      k(i) = (d(i) - c(i)*k(i-1)) / dem
      l(i) = a(i) / dem
    enddo
    den(i1) = (k(i1) - l(i1)*fluxj) / (1. + l(i1)*fac)
    den(jmax) = den(i1)
    do kk = 1, jmax-3
      jk = i1 - kk
      den(jk) = k(jk) - l(jk) * den(jk + 1)
    enddo

    return

  end subroutine impit

max_stars_repo_name: scivision/NCAR-GLOW

max_issues_repo_name: space-physics/NCAR-GLOW

max_forks_repo_name: space-physics/NCAR-GLOW


text:
!#############################################################################
	subroutine SEM
!#############################################################################
	USE IFPORT
	USE module_SEM
	USE vars
	USE multidata
	USE mpi
	IMPLICIT NONE
	DOUBLE PRECISION :: VOL,Ly,Lz,ENNE,XMIN,XMAX,YMIN,TI
	DOUBLE PRECISION :: YMAX,ZMIN,ZMAX,PI,U0,HU,RIZ,RDIVz,UAVE(3)
	DOUBLE PRECISION :: SIGMA_VALUE,NE_SEM,MAXVSEM,MINVSEM,rrandom
	INTEGER :: DIVY,DIVZ,IY,IZ,II,N,I,J,M,IT,IGLOBAL
	INTEGER,allocatable,dimension(:)::lsy,lsz,ley,lez
!THE BOX DIMENSIONS ARE DIFINED AS [XLENGHT] * [Ly] * [Lz]
![DIVX], [DIVY] & [DIVZ] ARE THE NUMBERS OF SPACIAL POINTS.
	PI = 2*ACOS(0.0D0) ; Ly  = yen-yst ; Lz= zen-zst
	U0 = ubulk 	; TI=TI_SEM

	ALLOCATE(iddom(jdom*kdom),ljdom(jdom*kdom),lkdom(jdom*kdom))
	ALLOCATE(elemyst(jdom*kdom),elemzst(kdom*jdom))
	ALLOCATE(elemyen(jdom*kdom),elemzen(kdom*jdom))
	ALLOCATE(ley(jdom),lez(kdom),lsy(jdom),lsz(kdom))

	 ljdom=0 	; lkdom=0 	;  DIVY = 0 	; DIVZ = 0
	 elemyst=0 	; elemyen=0	;  elemzst=0 	; elemzen=0
	 I=0 ;ley=0 ; lez=0	;  lsy=0 		; lsz=0

!ID OF THE BLOCK
	 DO N=1,kdom ; DO J=1,jdom
	 I=I+1
	  iddom(I)=idom*(J-1)+(jdom*idom*(N-1))
	 ENDDO ; ENDDO
!DIVISION ON THE Y AND Z DIRECTIONS
        lsy(1)=1
	  ley(1)=lsy(1)+NINT((ycor(iddom(1),2)-ycor(iddom(1),1))/g_dy)+1
	 Do J=2,jdom
	  lsy(J)=ley(J-1)+1
	  ley(J)=lsy(J)+NINT((ycor(iddom(J),2)-ycor(iddom(J),1))/g_dy)+1
	  DIVY = DIVY + NINT((ycor(iddom(J),2)-ycor(iddom(J),1))/g_dy) +1
	 Enddo
        lsz(1)=1
	  lez(1)=lsz(1)+NINT((zcor(iddom(1),2)-zcor(iddom(1),1))/g_dz)+1
	 DO N=1,kdom
	  lsz(N)=lez(N-1)+1
	  lez(N)=lsz(N)+NINT((zcor(iddom(N),2)-zcor(iddom(N),1))/g_dz)+1
	  DIVZ = DIVZ + NINT((zcor(iddom(N),2)-zcor(iddom(N),1))/g_dz) +1
	 ENDDO
!THE DIVISIONS FOR EACH OF THE DOMAINS IS DETERMINED:
	 I=0
	 DO N=1,kdom  ;  DO J=1,jdom
	 I=I+1
	  elemyst(I) = lsy(J)  ;  elemzst(I) = lsz(N)
	  elemyen(I) = ley(J)  ;  elemzen(I) = lez(N)
	  ljdom(I) = elemyen(I) -  elemyst(I) + 1
	  lkdom(I) = elemzen(I) -  elemzst(I) + 1
	 enddo ; enddo

	write(6,*)'Divisions :',DIVY,DIVZ, ' Number of BLOCKS  :',jdom,kdom
	write(6,'(10i4)')lsy(:)
	write(6,'(10i4)')ley(:)-lsy(:)
	write(6,'(10i4)')lsz(:)
	write(6,'(10i4)')lez(:)-lsz(:)
!THE AMPLITUDE OF THE VORTICES--> turbulence integral length scale
	allocate (SIGMA(3,DIVY,DIVZ))

!    ******     -> change numbers accordingly      ******  !
	SIGMA(1,:,:) = 0.252D0
	SIGMA(2,:,:) = 0.1485d0
	SIGMA(3,:,:) = 0.1125d0

!THE NUMBER OF EDDIES IS THE INLET SURFACE DIVIDED BY THE SURFACE OF EACH TURBULENT SPOT
![N](INTEGER) AND [ENNE](DOUBLE PRECISION) REPRESENT THE NUMBER OF EDDIES
	NE_SEM = Ly/Sigma(2,1,1)*Lz/Sigma(3,1,1)

	N = INT(NE_SEM)	;  ENNE = REAL(N)

 	if(myrank.eq.0) write(6,*)'The number of SEM eddies is :',N

!  [REYNOLDS(6)] IS A VECTOR WITH THE SIX ELEMENTS OF REYNOLDS STRESSES.
!  |REYNOLDS(1)  REYNOLDS(2)  REYNOLDS(4)|
!  |REYNOLDS(2)  REYNOLDS(3)  REYNOLDS(5)|
!  |REYNOLDS(4)  REYNOLDS(5)  REYNOLDS(6)|
	REYNOLDS=(/(TI*U0)**2, 0.0D0,(TI*U0)**2, 0.0D0, 0.0D0,(TI*U0)**2/)  ![M/S]

!ALLOCATION OF THE EDDIES VECTOR.
![VSEM(DIVX,DIVY,DIVZ,3)] IS THE INSTANTANEOUS VELOCITY VECTOR IN THE POINT WITH X,Y,Z COMPONENTS
![X_EDDY(3,N)] IS THE N-TH EDDY LOCATION X,Y,Z
![EPSILO(3,N)] IS THE N-TH EDDY INTENSITY IN X,Y,Z
![MOLT(3,N)]   IS THE MATRIX PRODUCT [R(3,3)] * [EPSILO(3,N)]
![Ksem(N)]     IS A VECTOR USED TO CHECK WHITCH EDDY IS OUTSIDE THE BOX AFTER THE CONVECTION

	ALLOCATE(Vsem(DIVY,DIVZ,3),Usem(DIVY,DIVZ,3),Ksem(N))
	ALLOCATE(X_EDDY(3,N),EPSILO(3,N),MOLT(3,N))
	PRINT *, "SEM total time interval = ", DT * ITMAX_SEM, " (s)"

	XMIN=1.D8 ; XMAX=1.D-8 ; YMIN=1.D8 ; YMAX=1.D-8
	ZMIN=1.D8 ; ZMAX=1.D-8

!DEFINITION OF EDDY LENGTH SCALE AND INITIAL VELOCITY FIELD
	DO IY=1,DIVY
	  DO IZ=1,DIVZ
	   Usem(IY,IZ,:)=(/ U0 ,0.D0,0.D0/)
	   UAVE(:) = UAVE(:) + Usem(IY,IZ,:)
!CALCULATION OF THE EDDY BOX PARAMETERS
	XMIN=MIN(0.0D0-SIGMA(1,IY,IZ),XMIN);XMAX=MAX(0.0D0+SIGMA(1,IY,IZ),XMAX)
	YMIN=MIN(0.0D0-SIGMA(2,IY,IZ),YMIN);YMAX=MAX(Ly   +SIGMA(2,IY,IZ),YMAX)
	ZMIN=MIN(0.0D0-SIGMA(3,IY,IZ),ZMIN);ZMAX=MAX(Lz   +SIGMA(3,IY,IZ),ZMAX)
	  END DO
	END DO

	UAVE = UAVE / (DIVY * DIVZ)
	VOL  = (XMAX - XMIN) * (YMAX - YMIN) * (ZMAX - ZMIN)

	call RANDOM_SEED()
!GENERATION OF THE EDDY LOCATION INSIDE THE BOX AND INITIALIZATION OF THE [Ksem] VECTOR
	DO II=1,N

		call RANDOM_NUMBER(rrandom)
	  X_EDDY(1,II) = (XMAX - XMIN) * rrandom + XMIN
	  X_EDDY(2,II) = (YMAX - YMIN) * rrandom + YMIN
	  X_EDDY(3,II) = (ZMAX - ZMIN) * rrandom + ZMIN
	  Ksem(II) = 0
!INITIALIZATION OF THE INTENSITIES. FOR EVERY DIRECTION THE AVERAGE INTENSITY VALUE IS CALCULATED AND
!IT IS FORCED TO BE LOWER THAN THE [VLIM] VALUE
		call RANDOM_NUMBER(rrandom)
	  EPSILO(1,II) = (rrandom*2.0D0 - 1.0D0)
	  EPSILO(2,II) = (rrandom*2.0D0 - 1.0D0)
	  EPSILO(3,II) = (rrandom*2.0D0 - 1.0D0)
	END DO
!INITIALIZATION OF THE [R(3,3)] MATRIX WITH THE CHOLENSKY DECOMPOSITION
!OF THE REYNOLDS STRESS TENSOR
	R = 0.d0
	R(1,1) = DSQRT(REYNOLDS(1))
	R(2,1) = REYNOLDS(2) / R(1,1)
	R(2,2) = DSQRT(REYNOLDS(3) - R(2,1)*R(2,1))
	R(3,1) = REYNOLDS(4) / R(1,1)
	R(3,2) = (REYNOLDS(5) - R(2,1)*R(3,1)) / R(2,2)
	R(3,3) = DSQRT(REYNOLDS(6) - R(3,1)*R(3,1) - R(3,2)*R(3,2))
!BEGINNING OF TIME ITERATIONS
	DO IT=1,ITMAX_SEM			!PARALLELISE THIS LOOP
!PRINTINGS OF GLOBAL VELOCITY AND OF CONVECTION VELOCITY
	 Do I=1,jdom*kdom
		  IGLOBAL=500+I
	   WRITE (FILEGLOBAL,'(A13,i4.4,a1,I5.5,A4)')
     &		'inflow/Inlet_',iddom(I),'_',IT,'.dat'
	   OPEN(UNIT=IGLOBAL, FILE=FILEGLOBAL,STATUS="UNKNOWN", ACTION="WRITE")
	   WRITE (IGLOBAL,*)'Variables=up,vp,wp'
	   WRITE (IGLOBAL,*)
     &  'zone  i=',ljdom(I),', j=',lkdom(I),', k=1 f=point'
	 Enddo

	  MOLT = MATMUL(R,EPSILO) ! !aij*epsij   MATRIX MULTIPLICATION
	  MAXVSEM=0.D0 ; MINVSEM=1000000.D0
!------BEGINNING OF SPATIAL ITERATION
	  DO IY = 1,DIVY
	    DO IZ = 1,DIVZ
!X_POINT = GRID POINT COORDINATES
	   X_POINT=(/0.0D0,IY*Ly/DIVY+YMIN,IZ*Lz/DIVZ+ZMIN/)
	   Vsem(IY,IZ,:)=(/ 0.d0, 0.d0,  0.d0 /)
!------------BEGINNING OF EDDIES ITERATIONS
		DO II=1,N
		TEMP(:) = DABS(X_POINT(:) - X_EDDY(:,II))
	  IF (TEMP(1).LT.SIGMA(1,IY,IZ) .AND.
     &      TEMP(2).LT.SIGMA(2,IY,IZ) .AND.
     &      TEMP(3).LT.SIGMA(3,IY,IZ)) THEN
	        Vsem(IY,IZ,:)=Vsem(IY,IZ,:)+MOLT(:,II)
     &        *(DSQRT(1.5D0)**3.0D0)*DSQRT(VOL)
     &        *(1.0D0-TEMP(1)/SIGMA(1,IY,IZ))
     &        *(1.0D0-TEMP(2)/SIGMA(2,IY,IZ))
     &        *(1.0D0-TEMP(3)/SIGMA(3,IY,IZ))
     &        /DSQRT(SIGMA(1,IY,IZ)*SIGMA(2,IY,IZ)*SIGMA(3,IY,IZ))
		 END IF
		END DO
! Instananeous velocity=mean velocity + SEM fluctuation velocity
	  Vsem(IY,IZ,:) = Vsem(IY,IZ,:) / DSQRT(ENNE)
	  if(Vsem(IY,IZ,1).gt.MAXVSEM) MAXVSEM=Vsem(IY,IZ,1)
	  if(Vsem(IY,IZ,1).le.MINVSEM) MINVSEM=Vsem(IY,IZ,1)
	    END DO
	  END DO

	 IF(MOD(IT,50).eq.0) write(6,'(a,i6,a8,3f13.6)')
     &   'ITERATION ',IT,' Time(s)',(IT-1)*dt,MAXVSEM,MINVSEM
!------END OF SPATIAL ITERATIONS
	      Do I=1,jdom*kdom
	        IGLOBAL=500+I
		   Do M=elemzst(I),elemzen(I)
		    Do J=elemyst(I),elemyen(I)
                  WRITE(500+I,'(3E15.6)')Vsem(J,M,:)
	  	  enddo ;enddo
	    	 CLOSE (UNIT=IGLOBAL)
	  	Enddo
!--------BEGINNING OF EDDIES CONVECTION ITERATIONS
	  DO II=1,N
!RE-CALCULATION OF THE EDDIES POSITION. IF ANY EDDY GOES BEYOND THE BOX LIMITS
!IT IS RESTARTED AT THE SURFACE FACING THE EXIT.
	    X_EDDY(1,II) = X_EDDY(1,II) + UAVE(1) * dt
	    X_EDDY(2,II) = X_EDDY(2,II) + UAVE(2) * dt
	    X_EDDY(3,II) = X_EDDY(3,II) + UAVE(3) * dt
!AFTER THE EDDIES CONVECTIONS ARE NECESSARY SOME TESTS TO CHECK IF ANY EDDY IS NOW OUTSIDE
!THE SEM BOX DEFINED EARLIER
!WHEN A EDDY IS RE-GENERATE THE Ksem(I) FACTOR ASSUME THE 1 VALUE. THIS VALUE IS USED LATER TO
!GENERATE A NEW INTENSITY FOR THE NEW EDDY

	    IF (X_EDDY(1,II) > XMAX) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = XMIN
		X_EDDY(2,II) = (YMAX - YMIN) * rrandom + YMIN
		X_EDDY(3,II) = (ZMAX - ZMIN) * rrandom + ZMIN
		Ksem(II) = 1
	    ELSE IF (X_EDDY(1,II) < XMIN) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = XMAX
		X_EDDY(2,II) = (YMAX - YMIN) * rrandom + ZMIN
		X_EDDY(3,II) = (ZMAX - ZMIN) * rrandom + ZMIN
		Ksem(II) = 1
	    ELSE IF (X_EDDY(2,II) > YMAX) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = (XMAX - XMIN) * rrandom + XMIN
		X_EDDY(2,II) = YMIN
		X_EDDY(3,II) = (ZMAX - ZMIN) * rrandom + ZMIN
		Ksem(II) = 1
	    ELSE IF (X_EDDY(2,II) < YMIN) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = (XMAX - XMIN) * rrandom + XMIN
		X_EDDY(2,II) = YMAX
		X_EDDY(3,II) = (ZMAX - ZMIN) * rrandom + ZMIN
		Ksem(II) = 1
	    ELSE IF  (X_EDDY(3,II) > ZMAX) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = (XMAX - XMIN) * rrandom + XMIN
		X_EDDY(2,II) = (YMAX - YMIN) * rrandom + ZMIN
		X_EDDY(3,II) = ZMIN
		Ksem(II) = 1
	    ELSE IF (X_EDDY(3,II) < ZMIN) THEN
		call RANDOM_NUMBER(rrandom)
		X_EDDY(1,II) = (XMAX - XMIN) * rrandom + XMIN
		X_EDDY(2,II) = (YMAX - YMIN) * rrandom + ZMIN
		X_EDDY(3,II) = ZMAX
		Ksem(II) = 1
	    END IF
!INTENSITY GENERATION FOR THE RE-CREATED EDDIES. WE ARE USING THE Ksem FACTOR AS EXPLAINED FATOR.
	    IF (Ksem(II)== 1) THEN
		call RANDOM_NUMBER(rrandom)
		EPSILO(3,II) = (rrandom*2.0D0 - 1.0D0)
		EPSILO(2,II) = (rrandom*2.0D0 - 1.0D0)
		EPSILO(1,II) = (rrandom*2.0D0 - 1.0D0)
	    END IF
	    Ksem(II) = 0
	  END DO
!---------END OF EDDIES CONVECTIONS ITERATIONS
	END DO
!----END OF TIME ITERATIONS

	END SUBROUTINE

ext: for


text:
      subroutine runY_0000i(k,l,i1,Xtwiddle,Gtwiddle,Shat5zz,N0)
      implicit none
C---  Expression for D0000i obtained from 5.50, following the comment after
C     5.50 on how to add adding additional "00" pairs
C---  (similar to Eq. 5.56b but with "00" added)
      include 'pvDnames.f'
      include 'pvDv.f'
      include 'Darraydef.f'
      include 'Darrays.f'
      integer ep,N0,k,l,i1,np
      parameter(np=3)
      double precision Xtwiddle(0:np,0:np),Gtwiddle(np,np)
      double complex Shat5zz(np,z2max,-2:0)

      if ((i1 .eq. l) .or. (i1 .eq. 0)) then
      return
      endif

      do ep=-2,0

      Dv(dzzzzi(i1)+N0,ep)=
     . (-2*Gtwiddle(k,i1)*Dv(dzzzzi(l)+N0,ep)
     . +Gtwiddle(k,1)*Shat5zz(1,z2(l,i1),ep)
     . +Gtwiddle(k,2)*Shat5zz(2,z2(l,i1),ep)
     . +Gtwiddle(k,3)*Shat5zz(3,z2(l,i1),ep)
     . +Xtwiddle(0,k)*Dv(dzzii(z2(l,i1))+N0,ep)
     . -Xtwiddle(0,0)*Dv(dzziii(z3(k,l,i1))+N0,ep))/(2*Gtwiddle(k,l))

      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      subroutine killop

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine performs the kill operation

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name         |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bio_ms(:)    |kg/ha         |land cover/crop biomass (dry weight)
!!    curyr        |none          |current year of simulation
!!    hrupest(:)  |none           |pesticide use flag:
!!                                | 0: no pesticides used in HRU
!!                                | 1: pesticides used in HRU
!!    icr(:)       |none          |sequence number of crop grown within the
!!                                |current year
!!    ihru         |none          |HRU number
!!    ncrops(:,:,:)|
!!    npmx        |none           |number of different pesticides used in
!!                                |the simulation
!!    nro(:)       |none          |sequence number for year in rotation
!!    nyskip       |none          |number of years to skip output printing/
!!                                |summarization
!!    plt_pst(:,:)|kg/ha          |pesticide on plant foliage
!!    sol_fon(:,:) |kg N/ha       |amount of nitrogen stored in the fresh
!!                                |organic (residue) pool
!!    sol_fop(:,:) |kg P/ha       |amount of phosphorus stored in the fresh
!!                                |organic (residue) pool
!!    sol_pst(:,:,1)|kg/ha        |pesticide in first layer of soil
!!    sol_rsd(:,:) |kg/ha         |amount of organic matter in the soil
!!                                |classified as residue
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name         |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bio_ms(:)    |kg/ha         |land cover/crop biomass (dry weight)
!!    idorm(:)     |none          |dormancy status code:
!!                                |0 land cover growing (not dormant)
!!                                |1 land cover dormant
!!    igro(:)      |none          |land cover status code:
!!                                |0 no land cover currently growing
!!                                |1 land cover growing
!!    laiday(:)    |m**2/m**2     |leaf area index
!!    ncrops(:,:,:)|
!!    phuacc(:)    |none          |fraction of plant heat units accumulated
!!    plantn(:)    |kg N/ha       |amount of nitrogen in plant biomass
!!    plantp(:)    |kg P/ha       |amount of phosphorus in plant biomass
!!    plt_pst(:,:) |kg/ha         |pesticide on plant foliage
!!    sol_fon(:,:) |kg N/ha       |amount of nitrogen stored in the fresh
!!                                |organic (residue) pool
!!    sol_fop(:,:) |kg P/ha       |amount of phosphorus stored in the fresh
!!                                |organic (residue) pool
!!    sol_pst(:,:,1)|kg/ha        |pesticide in first layer of soil
!!    sol_rsd(:,:) |kg/ha         |amount of organic matter in the soil
!!                                |classified as residue
!!    strsw(:)     |none          |fraction of potential plant growth achieved
!!                                |on the day where the reduction is caused by
!!                                |water stress
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    j           |none          |HRU number
!!    k           |none          |counter
!!    resnew      |
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Max

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer :: j, k
      real :: resnew

      !!by zhang
      !!====================
      real :: BLG1, BLG2, BLG3,  CLG, sf
      real :: sol_min_n, resnew_n, resnew_ne
      real :: LMF, LSF, LSLF, LSNF,LMNF
      orgc_f = 0.
      BLG1 = 0.
      BLG2 = 0.
      BLG3 = 0.
      CLG = 0.
      sf = 0.
      sol_min_n = 0.
      resnew = 0.
      resnew_n = 0.
      resnew_ne = 0.
      LMF = 0.
      LSF = 0.
      LSLF = 0.
      LSNF = 0.
      LMNF = 0.
      !!by zhang
      !!====================

      j = 0
      j = ihru

!      if (curyr > nyskip) then
!       ncrops(icr(j),j) = ncrops(icr(j),j) + 1
!      endif

	!! 22 January 2008
	resnew = 0.
	rtresnew = 0.
      resnew = bio_ms(j) * (1. - rwt(j))
	rtresnew = bio_ms(j) * rwt(j)
	call rootfr

	!! update residue, N, P on soil surface
      sol_rsd(1,j) = resnew + sol_rsd(1,j)
      sol_fon(1,j) = plantn(j) * (1. - rwt(j)) + sol_fon(1,j)
      sol_fop(1,j) = plantp(j) * (1. - rwt(j)) + sol_fop(1,j)
      sol_rsd(1,j) = Max(sol_rsd(1,j),0.)
	sol_fon(1,j) = Max(sol_fon(1,j),0.)
	sol_fop(1,j) = Max(sol_fop(1,j),0.)

            !!insert new biomss by zhang
            !!=================================
            if (cswat == 2) then
	          !!all the lignin from STD is assigned to LSL,
	            !!add STDL calculation
	          !!
	          !sol_LSL(k,ihru) = sol_STDL(k,ihru)
	          !CLG=BLG(3,JJK)*HUI(JJK)/(HUI(JJK)+EXP(BLG(1,JJK)-BLG(2,JJK)*&HUI(JJK))
	          ! 52  BLG1 = LIGNIN FRACTION IN PLANT AT .5 MATURITY
                ! 53  BLG2 = LIGNIN FRACTION IN PLANT AT MATURITY
                !CROPCOM.dat BLG1 = 0.01 BLG2 = 0.10
                !SUBROUTINE ASCRV(X1,X2,X3,X4)
                !EPIC0810
                !THIS SUBPROGRAM COMPUTES S CURVE PARMS GIVEN 2 (X,Y) POINTS.
                !USE PARM
                !XX=LOG(X3/X1-X3)
                !X2=(XX-LOG(X4/X2-X4))/(X4-X3)
                !X1=XX+X3*X2
                !RETURN
                !END
                !HUI(JJK)=HU(JJK)/XPHU

                BLG1 = 0.01/0.10
                BLG2 = 0.99
                BLG3 = 0.10
                XX = log(0.5/BLG1-0.5)
                BLG2 = (XX -log(1./BLG2-1.))/(1.-0.5)
                BLG1 = XX + 0.5*BLG2
                CLG=BLG3*phuacc(j)/(phuacc(j)+EXP(BLG1-BLG2*phuacc(j)))

	          !if (k == 1) then
		        sf = 0.05
	          !else
		        !sf = 0.1
	          !end if

               !kg/ha
	          sol_min_n = 0.
	          sol_min_n = (sol_no3(1,j)+sol_nh3(1,j))

	          resnew = resnew
	          resnew_n = ff1 * (plantn(j) - yieldn)
        	    resnew_ne = resnew_n + sf * sol_min_n

       	        !Not sure 1000 should be here or not!
        	    !RLN = 1000*(resnew * CLG/(resnew_n+1.E-5))
        	    RLN = (resnew * CLG/(resnew_n+1.E-5))
        	    RLR = MIN(.8, resnew * CLG/(resnew+1.E-5))

        	    LMF = 0.85 - 0.018 * RLN
        	    if (LMF <0.01) then
        	        LMF = 0.01
        	    else
        	        if (LMF >0.7) then
        	            LMF = 0.7
        	        end if
        	    end if
	          !if ((resnew * CLG/(resnew_n+1.E-5)) < 47.22) then
		        !    LMF = 0.85 - 0.018 * (resnew * CLG/(resnew_n+1.E-5))
	          !else
		        !    LMF = 0.
	          !end if

	          LSF =  1 - LMF

	          sol_LM(1,j) = sol_LM(1,j) + LMF * resnew
	          sol_LS(1,j) = sol_LS(1,j) + LSF * resnew

	          !here a simplified assumption of 0.5 LSL
	          LSLF = 0.0
	          LSLF = CLG

	          sol_LSL(1,j) = sol_LSL(1,j) + RLR* LSF * resnew
	          sol_LSC(1,j) = sol_LSC(1,j) + 0.42*LSF * resnew

	          sol_LSLC(1,j) = sol_LSLC(1,j) + RLR*0.42*LSF * resnew
	          sol_LSLNC(1,j) = sol_LSC(1,j) - sol_LSLC(1,j)

                !X3 = MIN(X6,0.42*LSF * resnew/150)

	          if (resnew_n >= (0.42 * LSF * resnew /150)) then
		         sol_LSN(1,j) = sol_LSN(1,j) + 0.42 * LSF * resnew / 150
		         sol_LMN(1,j) = sol_LMN(1,j) + resnew_n -
     &                         (0.42 * LSF * resnew / 150) + 1.E-25
	          else
		         sol_LSN(1,j) = sol_LSN(1,j) + resnew_n
		         sol_LMN(1,j) = sol_LMN(1,j) + 1.E-25
	          end if

	          !LSNF = sol_LSN(1,j)/(sol_LS(1,j)+1.E-5)

	          sol_LMC(1,j) = sol_LMC(1,j) + 0.42 * LMF * resnew
	          !LMNF = sol_LMN(1,j)/(sol_LM(1,j) + 1.E-5)

                !update no3 and nh3 in soil
                sol_no3(1,j) = sol_no3(1,j) * (1-sf)
                sol_nh3(1,j) = sol_nh3(1,j) * (1-sf)
            end if
            !!insert new biomss by zhang
            !!===============================

	!! allocate dead roots, N, P to soil layers
	do l=1, sol_nly(j)
	 sol_rsd(l,j) = sol_rsd(l,j) + rtfr(l) * rtresnew
	 sol_fon(l,j) = sol_fon(l,j) + rtfr(l) * plantn(j) * rwt(j)
	 sol_fop(l,j) = sol_fop(l,j) + rtfr(l) * plantp(j) * rwt(j)

              !!insert new biomss by zhang
              !!==============================
              if (cswat == 2) then
	          !!all the lignin from STD is assigned to LSL,
	            !!add STDL calculation
	          !!
	          !sol_LSL(k,ihru) = sol_STDL(k,ihru)
	          !CLG=BLG(3,JJK)*HUI(JJK)/(HUI(JJK)+EXP(BLG(1,JJK)-BLG(2,JJK)*&HUI(JJK))
	          ! 52  BLG1 = LIGNIN FRACTION IN PLANT AT .5 MATURITY
                ! 53  BLG2 = LIGNIN FRACTION IN PLANT AT MATURITY
                !CROPCOM.dat BLG1 = 0.01 BLG2 = 0.10
                !SUBROUTINE ASCRV(X1,X2,X3,X4)
                !EPIC0810
                !THIS SUBPROGRAM COMPUTES S CURVE PARMS GIVEN 2 (X,Y) POINTS.
                !USE PARM
                !XX=LOG(X3/X1-X3)
                !X2=(XX-LOG(X4/X2-X4))/(X4-X3)
                !X1=XX+X3*X2
                !RETURN
                !END
                !HUI(JJK)=HU(JJK)/XPHU

                BLG1 = 0.01/0.10
                BLG2 = 0.99
                BLG3 = 0.10
                XX = log(0.5/BLG1-0.5)
                BLG2 = (XX -log(1./BLG2-1.))/(1.-0.5)
                BLG1 = XX + 0.5*BLG2
                CLG=BLG3*phuacc(j)/(phuacc(j)+EXP(BLG1-BLG2*phuacc(j)))

	          if (l == 1) then
		        sf = 0.05
	          else
		        sf = 0.1
	          end if

               !kg/ha
	          sol_min_n = 0.
	          sol_min_n = (sol_no3(l,j)+sol_nh3(l,j))

	          resnew = rtfr(l) * rtresnew
	          resnew_n = rtfr(l) *ff2 * (plantn(j) - yieldn)
        	    resnew_ne = resnew_n + sf * sol_min_n
        	        !Not sure 1000 should be here or not!
        	    !RLN = 1000*(resnew * CLG/(resnew_n+1.E-5))
        	    RLN = (resnew * CLG/(resnew_n+1.E-5))
        	    RLR = MIN(.8, resnew * CLG/1000/(resnew/1000+1.E-5))

        	    LMF = 0.85 - 0.018 * RLN
        	    if (LMF <0.01) then
        	        LMF = 0.01
        	    else
        	        if (LMF >0.7) then
        	            LMF = 0.7
        	        end if
        	    end if
	          !if ((resnew * CLG/(resnew_n+1.E-5)) < 47.22) then
		        !    LMF = 0.85 - 0.018 * (resnew * CLG/(resnew_n+1.E-5))
	          !else
		        !    LMF = 0.
	          !end if

	          LSF =  1 - LMF

	          sol_LM(l,j) = sol_LM(l,j) + LMF * resnew
	          sol_LS(l,j) = sol_LS(l,j) + LSF * resnew

	          !here a simplified assumption of 0.5 LSL
	          !LSLF = 0.0
	          !LSLF = CLG

	          sol_LSL(l,j) = sol_LSL(l,j) + RLR*resnew
	          sol_LSC(l,j) = sol_LSC(l,j) + 0.42*LSF * resnew

	          sol_LSLC(l,j) = sol_LSLC(l,j) + RLR*0.42*resnew
	          sol_LSLNC(l,j) = sol_LSC(l,j) - sol_LSLC(l,j)

                !X3 = MIN(X6,0.42*LSF * resnew/150)

	          if (resnew_ne >= (0.42 * LSF * resnew /150)) then
		         sol_LSN(l,j) = sol_LSN(l,j) + 0.42 * LSF * resnew / 150
		         sol_LMN(l,j) = sol_LMN(l,j) + resnew_ne -
     &                         (0.42 * LSF * resnew / 150) + 1.E-25
	          else
		         sol_LSN(l,j) = sol_LSN(l,j) + resnew_ne
		         sol_LMN(l,j) = sol_LMN(l,j) + 1.E-25
	          end if

	          !LSNF = sol_LSN(l,j)/(sol_LS(l,j)+1.E-5)

	          sol_LMC(l,j) = sol_LMC(l,j) + 0.42 * LMF * resnew
	          !LMNF = sol_LMN(l,j)/(sol_LM(l,j) + 1.E-5)

                !update no3 and nh3 in soil
                sol_no3(l,j) = sol_no3(l,j) * (1-sf)
                sol_nh3(l,j) = sol_nh3(l,j) * (1-sf)
            end if
            !!insert new biomss by zhang
            !!===============================

	end do

      if (hrupest(j) == 1) then
        do k = 1, npmx
          sol_pst(k,j,1) = sol_pst(k,j,1) + plt_pst(k,j)
          plt_pst(k,j) = 0.
        end do
      end if

      bio_hv(icr(j),j) = bio_ms(j) + bio_hv(icr(j),j)
      !bio_yrms(j) = bio_yrms(j) + bio_ms(j) / 1000.

	!! reset variables
      igro(j) = 0
      idorm(j) = 0
      bio_ms(j) = 0.
      rwt(j) = 0.
	plantn(j) = 0.
      plantp(j) = 0.
      strsw(j) = 1.
      laiday(j) = 0.
      hvstiadj(j) = 0.
      phuacc(j) = 0.
!      phubase(j) = 0.
	rtfr = 0. ! Resetting roots fraction per layer array

      return
      end


text:
-makelib ies_lib/xilinx_vip -sv \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi4stream_vip_axi4streampc.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi_vip_axi4pc.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/xil_common_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi4stream_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi_vip_pkg.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi4stream_vip_if.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/axi_vip_if.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/clk_vip_if.sv" \
  "C:/Xilinx/Vivado/2019.1/data/xilinx_vip/hdl/rst_vip_if.sv" \
-endlib
-makelib ies_lib/xil_defaultlib -sv \
  "C:/Xilinx/Vivado/2019.1/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
  "C:/Xilinx/Vivado/2019.1/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib ies_lib/xpm \
  "C:/Xilinx/Vivado/2019.1/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies_lib/axi_infrastructure_v1_1_0 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_vip_v1_1_5 -sv \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/d4a8/hdl/axi_vip_v1_1_vl_rfs.sv" \
-endlib
-makelib ies_lib/processing_system7_vip_v1_0_7 -sv \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/8c62/hdl/processing_system7_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/design_1/ip/design_1_processing_system7_0_0/sim/design_1_processing_system7_0_0.v" \
-endlib
-makelib ies_lib/lib_cdc_v1_0_2 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/ef1e/hdl/lib_cdc_v1_0_rfs.vhd" \
-endlib
-makelib ies_lib/proc_sys_reset_v5_0_13 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/8842/hdl/proc_sys_reset_v5_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/design_1/ip/design_1_rst_ps7_0_100M_0/sim/design_1_rst_ps7_0_100M_0.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/fxd_arith_pkg.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/rtpm_pkg.vhd" \
  "../../../bd/design_1/ipshared/9bb0/hdl/PowerMonitoringIP_v1_0_S00_AXI.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/ac.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/calc.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/sp_ram.vhd" \
  "../../../bd/design_1/ipshared/9bb0/src/top.vhd" \
  "../../../bd/design_1/ipshared/9bb0/hdl/PowerMonitoringIP_v1_0.vhd" \
  "../../../bd/design_1/ip/design_1_PowerMonitoringIP_0_1/sim/design_1_PowerMonitoringIP_0_1.vhd" \
-endlib
-makelib ies_lib/generic_baseblocks_v2_1_0 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/b752/hdl/generic_baseblocks_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_4 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/1f5a/simulation/fifo_generator_vlog_beh.v" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_4 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/1f5a/hdl/fifo_generator_v13_2_rfs.vhd" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_4 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/1f5a/hdl/fifo_generator_v13_2_rfs.v" \
-endlib
-makelib ies_lib/axi_data_fifo_v2_1_18 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/5b9c/hdl/axi_data_fifo_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_register_slice_v2_1_19 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/4d88/hdl/axi_register_slice_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_protocol_converter_v2_1_19 \
  "../../../../PowerMonitoring.srcs/sources_1/bd/design_1/ipshared/c83a/hdl/axi_protocol_converter_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/design_1/ip/design_1_auto_pc_0/sim/design_1_auto_pc_0.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/design_1/sim/design_1.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  glbl.v
-endlib

max_stars_repo_name: al-ludwig/Run-Time-Power-Monitoring

max_issues_repo_name: al-ludwig/Run-Time-Power-Monitoring

max_forks_repo_name: al-ludwig/Run-Time-Power-Monitoring


    Subroutine x04ccfn(uplo, diag, n, a, title, errbuf, ifail)

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     Prints a real triangular matrix stored in packed vector form.
!     Easy to use driver for X04CDFN.

!     .. Use Statements ..
      Use lapack_example_aux, Only: p01abf, x04cdfn
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Character (6), Parameter :: srname = 'X04CCF'
!     .. Scalar Arguments ..
      Integer :: ifail, n
      Character (1) :: diag, uplo
      Character (200) :: errbuf
      Character (*) :: title
!     .. Array Arguments ..
      Real (Kind=dp) :: a(*)
!     .. Local Scalars ..
      Integer :: ierr, ifcdf, nrec
!     .. Local Arrays ..
      Character (1) :: clabs(1), rlabs(1)
      Character (80) :: rec(2)
!     .. Intrinsic Procedures ..
      Intrinsic :: char
!     .. Executable Statements ..
      Continue
      errbuf(1:1) = char(0)
!
!     Check for incorrect arguments.
!
      ierr = 0
      nrec = 0
!
      If (uplo/='L' .And. uplo/='l' .And. uplo/='U' .And. uplo/='u') Then
        ierr = 1
        nrec = 1
        Write (rec, Fmt=110) uplo
        Write (errbuf, Fmt=120) srname, 1, 99999, 1, uplo
      Else If (diag/='U' .And. diag/='u' .And. diag/='N' .And. diag/='n' .And. &
          diag/='B' .And. diag/='b') Then
        ierr = 2
        nrec = 1
        Write (rec, Fmt=130) diag
        Write (errbuf, Fmt=140) srname, 2, 99998, 1, diag
      End If
!
      If (ierr/=0) Go To 100
!
!     End of argument checking.
!
      ifcdf = 1
      Call x04cdfn(uplo, diag, n, a, ' ', title, 'I', rlabs, 'I', clabs, 80, &
        0, errbuf, ifcdf)
!
100   Continue
      errbuf(200:200) = char(0)
      ifail = p01abf(ifail, ierr, srname, nrec, rec)
!
      Return
!
110   Format (' ** On entry, UPLO is not valid: UPLO = ''', A, '''.')
120   Format (A, 3(1X,I5), 1X, A)
130   Format (' ** On entry, DIAG is not valid: DIAG = ''', A, '''.')
140   Format (A, 3(1X,I5), 1X, A)
    End Subroutine


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.

C $Id: ipopt_getdata.f,v 1.2 2002/07/24 03:00:53 andreasw Exp $

      subroutine IPOPT_GETDATA(FFUNCS, CFUNCS, NOCG, NORES, NONEGCURV)
      implicit none
      integer FFUNCS, CFUNCS, NOCG, NORES, NONEGCURV
      include 'TIMER.INC'
      include 'IPOPT.INC'
      integer fevals, cevals, nboundsfree
      common /EVALS/ fevals, cevals
      FFUNCS = fevals
      CFUNCS = cevals
      NOCG = COUNT_CG
      NORES = COUNT_RESTO_ITER
      NONEGCURV = COUNT_NEG_CURV
C
C     Also close output file if open
C
      if( QCNR.gt.0 ) close(QCNR)
      return
      end

      subroutine CHECK_FLAGFILE( EX )
      integer EX
      logical exist
      inquire(file='FLAG', exist=exist)
      if( exist ) then
         EX = 1
      else
         EX = 0
      endif
      return
      end

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE CHGINT (IOLD, INEW, LIST, LLIST)
C=======================================================================

C   --*** CHGINT *** (GJOIN) Changes all occurrences in list
C   --   Written by Amy Gilkey - revised 09/29/87
C   --
C   --CHGINT changes all occurrences of an integer in a list to a given
C   --integer.
C   --
C   --Parameters:
C   --   IOLD - IN - the integer to change
C   --   INEW - IN - the integer to replace IOLD
C   --   LIST - IN/OUT - the list of integers
C   --   LLIST - IN - the length of LIST

      INTEGER LIST(LLIST)

      DO 100 I = 1, LLIST
         IF (LIST(I) .EQ. IOLD) THEN
            LIST(I) = INEW
         END IF
  100 CONTINUE

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      SUBROUTINE ORDER (GPLST,ID,REST,GRIDS,IDTAB,LCOR,B1,B2,B3)
C
      LOGICAL         SPILL
      INTEGER         GRIDS(1),ID(1),IDTAB(2),TP,GPLST(1),IGRD(2),SCR4,
     1                ISYM(14),ITYPE(14),HOLD(3),ELID,SILS(34),REST(2),
     2                EST,SIL,SCR2,ECPT,B1,B2,B3,THREE(3),OFFSET
      COMMON /BLANK / NGP,SKP(11),EST,SKIP1(3),SIL,SKIP2(5),ECPT,OES1,
     1                SCR1,SCR2,NEWOES,SCR4
      EQUIVALENCE     (THREE(1),IFLAG),(THREE(2),NELMT),(THREE(3),IGDPT)
      EQUIVALENCE     (KQ4,ISYM(13)),(KT3,ISYM(14))
      DATA    ISYM  / 2HSH,2HT1,2HTB,2HTP,2HTM,2HQP,2HQM,2HT2,2HQ2,2HQ1,
     1                2HM1,2HM2,2HQ4,2HT3/    ,KBAR/2HBR/
      DATA    ITYPE / 4,6,7,8,9,15,16,17,18,19,62,63,64,83/
      DATA    NTYPE / 14 /
C
C     BUILD A TABLE FOR GPECT POINTERS TO ELID AND ITS ORDERED GRID PTS
C
      SPILL = .FALSE.
      J = 1
      I = 3
      IDTAB(1) = 0
      JSPILL = 1
      LCORX  = LCOR
      NEWIN  = SCR4
      NEWOUT = SCR2
 2    CALL READ (*130,*12,EST,TP,1,0,M)
      OFFSET = 0
      IF (TP .EQ. KBAR) OFFSET = 6
      IF (TP.EQ.KT3 .OR. TP.EQ.KQ4) OFFSET = 1
 3    CALL FREAD (EST,NGPPE,1,0)
      IDTAB(I-1) = NGPPE
C
C     SKIP PAST THE NON-CONTOUR ELEMENTS
C
      DO 4 K = 1,NTYPE
      IF (TP .EQ. ISYM(K)) GO TO 8
 4    CONTINUE
 6    CALL FREAD (EST,ELID,1,0)
      IF (ELID .EQ. 0) GO TO 2
      J = 1 + NGPPE + OFFSET
      CALL FREAD (EST,0,-J,0)
      GO TO 6
C
C     CONSTRUCT IDTAB  1. 0, 2.NGPPE, 3.ELID, 4.ELIDPTR, 5.REPEAT.  3,4
C     FOR ALL ELEMENTS OF THIS TYPE, 6.REPEAT 1-5 FOR ALL ELEMENTS IN
C     THE SET. CONSTRUCT GRIDS  1-NGPPE. GRIDS FOR 1ST ELEMENT, NEXT.
C     REPEAT 1ST FOR ALL ELEMENTS IN THE IDTAB
C
 8    CALL READ (*12,*12,EST,IDTAB(I),2,0,M)
      I = I + 2
      IF (IDTAB(I-2) .NE. 0) GO TO 10
C
C     END OF ELEMENTS OF THIS TYPE
C
      TP = IDTAB(I-1)
      GO TO 3
C
 10   CALL FREAD (EST,GRIDS(J),NGPPE,0)
      IF (OFFSET .NE. 0) CALL FREAD (EST,0,-OFFSET,0)
      J = J + NGPPE
      IF (I .GE. LCORX) GO TO 14
      GO TO 8
C
C     TABLE FIT INTO CORE
C
 12   CALL BCKREC (EST)
      GO TO 16
C
C     SPILL OCCURS - TABLE DID NOT FIT
C
 14   SPILL = .TRUE.
C
C     END OF TABLE
C
 16   LIDTAB = I - 1
      IF (LIDTAB .LE. 2) GO TO (130,140), JSPILL
      LGRIDS = J - 1
      LASTNG = NGPPE
      GO TO (18,140), JSPILL
 18   CALL OPEN (*130,ECPT,GPLST(B2),0)
      CALL GOPEN (SCR2,GPLST(B1),1)
      CALL FWDREC (*120,ECPT)
      IGDPT = 0
 20   IGDPT = IGDPT + 1
      IEOR  = 0
      IF (GPLST(IGDPT) .NE. 0) GO TO 25
      CALL FWDREC (*120,ECPT)
      GO TO 20
 25   NELMT = 0
      IFLAG =-1
C
C      ECPT--1. PIVOT POINT, 2. DEG.FREEDOM, 3. -LENGTH, 4. ELID POINTER
C      5. ELTYPE, 6.SILS (THERE ARE (LENGTH-2) OF THEM), 7. REPEAT ITEMS
C      (3-6) FOR ALL ELEMENTS ATTACHED TO PIVOT, 8. EOR, 9. REPEAT ITEMS
C      (1-8) FOR ALL GRIDS IN THE PROBLEM.
C
      CALL READ (*120,*120,ECPT,IGRD,2,0,M)
 30   CALL READ (*120,*75,ECPT,LENGTH,1,0,M)
      CALL FREAD (ECPT,SILS,-LENGTH,0)
      TP = SILS(2)
      DO 32 I = 1,NTYPE
      IF (TP .EQ. ITYPE(I)) GO TO 33
 32   CONTINUE
      GO TO 30
C
C     MATCH ELIDPTR WITH ITS ELID AND GRID POINTS IF POSSIBLE
C
 33   J = 1
      DO 50 I = 1,LIDTAB,2
      IF (IDTAB(I)) 40,35,40
 35   NGPPE = IDTAB(I+1)
      GO TO 50
 40   IF (IDTAB(I+1) .EQ. SILS(1)) GO TO 55
      J = J + NGPPE
 50   CONTINUE
C
C     IF NOT IN THE TABLE, IS THERE SPILL(IE IS TABLE NOT COMPLETE).
C     NO SPILL, SKIP HIM.  YES SPILL, FLAG HIM.
C
      IF (.NOT.SPILL) GO TO (30,145), JSPILL
      ELID  =-SILS(1)
      NELMT = NELMT + 1
      GO TO 70
C
C     FOUND ELEMENT IN THE TABLE
C
 55   ELID = IDTAB(I)
      DO 60 I = 1,NGPPE
      K = J + I - 1
      IF (IGDPT .EQ. GRIDS(K)) GO TO 65
 60   CONTINUE
 65   IAFTER = I - (I/NGPPE)*NGPPE + J
      IBEFOR = J + I - 2
      IF (I .EQ. 1) IBEFOR = IBEFOR + NGPPE
      NELMT = NELMT + 1
      REST(2*NELMT-1) = GRIDS(IAFTER)
      REST(2*NELMT  ) = GRIDS(IBEFOR)
 70   ID(NELMT) = ELID
      IF (NELMT .LT. LCOR/2) GO TO (30,145), JSPILL
      GO TO 80
 75   IF (NELMT .EQ. 0) GO TO (20,140), JSPILL
      IEOR = 1
C
C     ORDER ELEMENTS IF WE HAVE REACHED END OF EST FILE
C
 80   IF (SPILL) GO TO 112
      IF (NELMT .LE.2)  GO TO 110
      INDEX = 3
      IALL  = 2*NELMT
      IONE  = REST(1)
      ITWO  = REST(2)
 85   IF (IONE .EQ. ITWO) GO TO 105
      DO 90 I = INDEX,IALL,2
      IF (ITWO .EQ. REST(I)) GO TO 95
 90   CONTINUE
      GO TO 110
 95   IF (I .EQ. INDEX) GO TO 100
      J = (INDEX+1)/2
      K = (I+1)/2
      HOLD(1) = ID(J)
      ID(J)   = ID(K)
      ID(K)   = HOLD(1)
      HOLD(2) = REST(INDEX  )
      HOLD(3) = REST(INDEX+1)
      REST(INDEX  ) = REST(I  )
      REST(INDEX+1) = REST(I+1)
      REST(I  ) = HOLD(2)
      REST(I+1) = HOLD(3)
 100  INDEX = INDEX + 2
      ITWO  = REST(INDEX-1)
      IF (INDEX .LT. IALL) GO TO 85
      IF (IONE  .NE. ITWO) GO TO 110
C
C     INTERIOR ELEMENTS
C
 105  CALL WRITE (NEWOUT,THREE,3,0)
      CALL WRITE (NEWOUT,ID,NELMT,1)
      IF (IGDPT .LT. NGP) GO TO (20,140), JSPILL
      GO TO 120
C
C     BORDER ELEMENTS
C
 110  IFLAG = -2
 112  CALL WRITE (NEWOUT,THREE,3,0)
      J = -1
      DO 115 I = 1,NELMT
      J = J + 2
      CALL WRITE (NEWOUT,ID(I),1,0)
      CALL WRITE (NEWOUT,REST(J),2,0)
 115  CONTINUE
      IQ = 2*NELMT
      CALL WRITE (NEWOUT,0,0,1)
      GO TO (118,140), JSPILL
 118  IF (IEOR) 25,25,119
 119  IF (IGDPT .LT.NGP) GO TO 20
 120  CALL CLOSE (ECPT,1)
 125  CALL WRITE (NEWOUT,0,1,1)
      CALL CLOSE (NEWOUT,1)
C
C     IF NO SPILL -  RETURN
C
 130  IF (.NOT.SPILL) GO TO 170
C
C    COME HERE IF WE HAVE SPILL
C
      I = NEWOUT
      NEWOUT = NEWIN
      NEWIN  = I
      CALL GOPEN (NEWIN,GPLST(B1),0)
      CALL GOPEN (NEWOUT,GPLST(B2),1)
      JSPILL = 2
      NGPPE  = LASTNG
      IDTAB(1) = 0
      IDTAB(2) = LASTNG
      I = 3
      J = 1
      SPILL = .FALSE.
      GO TO 8
C
C     TABLE CONSTRUCTED SO RETURN HERE
C
 140  CALL READ (*160,*160,NEWIN,THREE,3,0,M)
      NELMT = 0
 145  CALL READ (*160,*75,NEWIN,SILS(1),3,0,M)
      IF (SILS(1)) 150,150,155
 150  SILS(1) = -SILS(1)
      GO TO 33
 155  ELID  = SILS(1)
      NELMT = NELMT + 1
      REST(2*NELMT-1) = SILS(2)
      REST(2*NELMT  ) = SILS(3)
      GO TO 70
C
C     END OF FILE
C
 160  CALL CLOSE (NEWIN,1)
      GO TO 125
C
C     OUTPUT FILE MUST BE SCRATCH 2
C
 170  IF (NEWOUT .EQ. SCR2) RETURN
      CALL GOPEN  (NEWOUT,GPLST(B1),0)
      CALL GOPEN  (SCR2,GPLST(B2),1)
      CALL CPYFIL (NEWOUT,SCR2,REST,LCOR,M)
      CALL CLOSE  (SCR2,1)
      CALL CLOSE  (NEWOUT,1)
      RETURN
      END


text:
subroutine GaussArrayRestart(a,array,itask)
   !------------------------------------------------------------------------
   !    This routine reads the initial values from the restart file
   !    corresponding to gauss point arrays and writes them to a new restart
   !------------------------------------------------------------------------
   use typre
   use def_parame
   use Mod_iofile
   use Mod_int2str
   use Mod_PhysicalProblem
   implicit none
   class(PhysicalProblem) :: a
   integer(ip), intent(in)    :: itask

   integer(ip) :: ndime,npoin,oldNpoin,ipoin,idofn,ndofn,ielem,nelem
   type(r3p)    :: array(:)
   real(rp)    :: rdummy
   integer(ip) :: iaux,iaux2,iaux3
   real(rp), allocatable :: auxvesgs(:)

   call a%Timer%Total%Tic
   call a%Timer%Restar%Tic

   select case (itask)

   case (1)

      call a%Mesh%GetNdime(ndime)
      call a%Mesh%GetNpoin(npoin)
      call a%Mesh%GetNelem(nelem)

      iaux = 0
      do ielem = 1,nelem
         iaux = iaux+size(array(ielem)%a,3)*ndime
      enddo
      call a%Memor%alloc(iaux,auxvesgs,'auxvesgs','nsi_restar')

      read (a%lun_rstar) (auxvesgs(ipoin),ipoin=1,OldNpoin)

      iaux = 0
      do ielem = 1,nelem
         iaux2 = size(array(ielem)%a,3)
         do iaux3 = 1,iaux2
            array(ielem)%a(1:ndime,1,iaux3)= auxvesgs(iaux+1:iaux+ndime)
            iaux = iaux+ndime
         enddo
      enddo
      call a%Memor%dealloc(iaux,auxvesgs,'auxvesgs','nsi_restar')

   case (2)

      call a%Mesh%GetNdime(ndime)
      call a%Mesh%GetNpoin(npoin)
      call a%Mesh%GetNelem(nelem)

      iaux = 0
      do ielem = 1,nelem
         iaux = iaux+size(array(ielem)%a,3)*ndime
      enddo
      call a%Memor%alloc(iaux,auxvesgs,'auxvesgs','nsi_restar')

      iaux = 0
      do ielem = 1,nelem
         iaux2 = size(array(ielem)%a,3)
         do iaux3 = 1,iaux2
           auxvesgs(iaux+1:iaux+ndime) = array(ielem)%a(1:ndime,1,iaux3)
           iaux = iaux+ndime
         enddo
      enddo

      write (a%lun_rstar) (auxvesgs(ipoin),ipoin=1,npoin)

      call a%Memor%dealloc(iaux,auxvesgs,'auxvesgs','nsi_restar')

   end select

   call a%Timer%Total%Toc
   call a%Timer%Restar%Toc

100 format(5x,'ERROR:   ',a)
101 format(5x,'WARNING: ',a)

end subroutine GaussArrayRestart

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


!
!     ------------------------------------------------------------------
!     F I N D
!     ------------------------------------------------------------------
!
      CHARACTER*(3) FUNCTION FIND (I, OF, EL)
!...Translated by Pacific-Sierra Research 77to90  4.3E  21:26:11  11/16/01
!...Switches:
!
!  ---  THIS ROUTINE FINDS ELECTRONS IN ONE OF THREE LISTS
!
      IMPLICIT NONE
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: NWD = 128
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: I
      CHARACTER , INTENT(IN) :: OF(NWD,2)*3
      CHARACTER , INTENT(IN) :: EL(NWD)*3
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!
      IF (I <= NWD) THEN
         FIND = EL(I)
      ELSE IF (I <= NWD*2) THEN
         FIND = OF(I-NWD,1)
      ELSE
         FIND = OF(I-2*NWD,2)
      ENDIF
      RETURN
      END FUNCTION FIND


text:
subroutine sld_outerr(a)
!------------------------------------------------------------------------
! NAME
!    sld_outerr
! DESCRIPTION
!    This routine checks if there are errors and warnings
!------------------------------------------------------------------------
  use typre
  use Mod_int2str
  use Mod_Solids
  implicit none
  integer(ip) :: ierro=0,iwarn=0

  class(SolidsProblem) :: a

   if(ierro==1) then
      call runend(adjustl(trim(int2str(ierro)))//' ERROR HAS BEEN FOUND')
   else if(ierro>=2) then
      call runend(adjustl(trim(int2str(ierro)))//' ERRORS HAVE BEEN FOUND')
   end if

!Formats
100 format(5x,'ERROR:   ',a)
101 format(5x,'WARNING: ',a)

end subroutine sld_outerr


text:
20     PRINT *, ' Enter x, l (l .lt. 0 to stop)'
       READ *, X,L

       IF (L .LT. 0) THEN
           STOP
       ELSE IF (L .EQ. 0) THEN
           PL=0.
       ELSE IF (L .EQ. 1) THEN
           PL=X
       ELSE
           PM=1.
           PZ=X
           DO 10 IL=1,L-1
              PP=((2*IL+1)*X*PZ-IL*PM)/(IL+1)
              PM=PZ
              PZ=PP
10         CONTINUE
           PL=PZ
       END IF

       PRINT *,X,L,PL
       GOTO 20
       END


text:
c=======================================================================
c=== RBio/RBrread_64 ===================================================
c=======================================================================

c RBio: a MATLAB toolbox for reading and writing sparse matrices in
c Rutherford/Boeing format.
c Copyright (c) 2007, Timothy A. Davis, Univ. of Florida

c-----------------------------------------------------------------------
c RB*read:  read a Rutherford/Boeing matrix
c-----------------------------------------------------------------------

	subroutine RBrread
     $	    (nrow, ncol, nnz, ptrfmt, indfmt, valfmt,
     $	    mkind, skind, Ap, Ai, Ax, nzeros, w, cp, info, nw, nnz1)
	integer*8
     $	    nrow, ncol, nnz, mkind, skind, p, j, i, alen, llen, k,
     $	    ilast, nzeros, info, nw, nnz1
	integer*8
     $	    Ap (ncol+1), Ai (nnz1), w (nw), cp (ncol+1)
	double precision Ax (nnz1), x
	character ptrfmt*16, indfmt*16, valfmt*20
	character*4 s

c	----------------------------------------------------------------
c	read the column pointers and row indices
c	----------------------------------------------------------------

	call RBpattern (ptrfmt, indfmt, nrow, ncol, nnz, skind,
     $	    Ap, Ai, w, cp, info, nw)
	if (info .ne. 0) then
c	    error: pattern is invalid
	    return
	endif

c	----------------------------------------------------------------
c	read the values
c	----------------------------------------------------------------

	if (mkind .eq. 1) then

c	    pattern-only matrix, set all values to 1
	    do 10 i = 1, nnz
		Ax (i) = 1
10	    continue

	elseif (mkind .eq. 3) then

c	    Read in nnz integer values, then convert them to real.
c	    use Ai as workspace.  If the matrix is symmetric, Ai is
c	    twice as big as nnz, so use Ai (nnz+1...2*nnz) as workspace.
c	    Otherwise, for an iua matrix, use Ai (1..nnz) and then
c	    rewind the file and read Ap and Ai back in again.

	    if (skind .le. 0) then
		k = 0
	    else
		k = nnz
	    endif
	    read (7, valfmt, err = 94, end = 94) (Ai(p), p = k+1, k+nnz)
	    do 5 p = 1, nnz
		Ax (p) = Ai (p+k)
5	    continue
	    if (skind .le. 0) then
c		now that Ai has been destroyed, rewind the file and read
c		in Ap and Ai again.  Skip the 4-line header first.
		rewind (unit = 7)
		read (7, 15) (s (k:k), k = 1,4)
15		format (a1 / a1 / a1 / a1)
		call RBpattern (ptrfmt, indfmt, nrow, ncol, nnz, skind,
     $		    Ap, Ai, w, cp, info, nw)
		if (info .ne. 0) then
c		    error: pattern is invalid.  This 'cannot' happen,
c		    because the pattern was already read in above.
		    return
		endif
	    endif

	else

c	    read nnz values
	    read (7, valfmt, err = 94, end = 94) (Ax (p), p = 1, nnz)

	endif

c	----------------------------------------------------------------
c	construct upper triangular part for symmetric matrices
c	----------------------------------------------------------------

c	If skind is zero, then the upper part is not constructed.  This
c	allows the caller to skip this part, and create the upper part
c	of a symmetric (S,H,Z) matrix.  Just pass skind = 0.

	if (skind .gt. 0) then

c	    ------------------------------------------------------------
c	    shift the matrix by adding gaps to the top of each column
c	    ------------------------------------------------------------

	    do 30 j = ncol, 1, -1

c		number of entries in lower tri. part (incl. diagonal)
		llen = Ap (j+1) - Ap (j)

c		number of entries in entire column
		alen = cp (j+1) - cp (j)

c		move the column from Ai (Ap(j) ... Ap(j+1)-1)
c		down to Ai (cp(j+1)-llen ... cp(j+1)-1), leaving a gap
c		at Ai (Ap(j) ... cp(j+1)-llen)

		do 20 k = 1, llen
		    Ai (cp (j+1) - k) = Ai (Ap (j+1) - k)
		    Ax (cp (j+1) - k) = Ax (Ap (j+1) - k)
20		continue
30	    continue

c	    ------------------------------------------------------------
c	    populate the upper triangular part
c	    ------------------------------------------------------------

c	    create temporary column pointers to point to the gaps
	    do 40 j = 1, ncol
		w (j) = cp (j)
40	    continue

	    do 60 j = 1, ncol

c		scan the entries in the lower tri. part, in
c		Ai (cp(j+1)-llen ... cp(j+1)-1)
		llen = Ap (j+1) - Ap (j)
		do 50 k = 1, llen

c		    get the A(i,j) entry in the lower triangular part
		    i = Ai (cp (j+1) - k)
		    x = Ax (cp (j+1) - k)

c		    add A(j,i) as the next entry in column i (excl diag)
		    if (i .ne. j) then
			p = w (i)
			w (i) = w (i) + 1
			Ai (p) = j

			if (skind .eq. 1) then
c			    *SA matrix
			    Ax (p) = x
			elseif (skind .eq. 2) then
c			    *HA matrix
			    Ax (p) = x
			else
c			    *ZA matrix
			    Ax (p) = -x
			endif

		    endif

50		continue
60	    continue

c	    finalize the column pointers
	    do 70 j = 1, ncol+1
		Ap (j) = cp (j)
70	    continue

	endif

c	----------------------------------------------------------------
c	count the number of explicit zeros
c	----------------------------------------------------------------

	nzeros = 0
	do 90 j = 1, ncol
	    cp (j) = nzeros + 1
	    do 80 p = Ap (j), Ap (j+1)-1
		if (Ax (p) .eq. 0) then
		    nzeros = nzeros + 1
		endif
80	    continue
90	continue
	cp (ncol+1) = nzeros + 1

c	----------------------------------------------------------------
c	matrix is valid
c	----------------------------------------------------------------

	info = 0
	return

c	----------------------------------------------------------------
c	error return
c	----------------------------------------------------------------

94	info = -94
	return
	end

c-----------------------------------------------------------------------
c RB*zeros: extract explicit zero entries
c-----------------------------------------------------------------------
c
c   nrow-by-ncol: size of A and Z
c   cp: column pointers of Z on input
c   Ap, Ai, Ax: matrix with zeros on input, pruned on output
c   Zp, Zi, Zx: empty matrix on input, pattern of zeros on output
c
c-----------------------------------------------------------------------

	subroutine RBrzeros
     $	    (nrow, ncol, cp, Ap, Ai, Ax, Zp, Zi, Zx)
	integer*8
     $	    nrow, ncol, Ap (ncol+1), Ai (*), Zp (ncol+1), Zi (*),
     $	    cp (*), i, j, p, pa, pz, p1
	double precision Ax (*), x
	double precision Zx (*)

c	----------------------------------------------------------------
c	copy the column pointers if Z is being constructed
c	----------------------------------------------------------------

	do 10 j = 1, ncol+1
	    Zp (j) = cp (j)
10	continue

c	----------------------------------------------------------------
c	split the matrix
c	----------------------------------------------------------------

	pa = 1
	pz = 1
	do 30 j = 1, ncol
c	    save the new start of column j of A
	    p1 = Ap (j)
	    Ap (j) = pa
	    pz = Zp (j)
c	    split column j of A
	    do 20 p = p1, Ap (j+1)-1
		i = Ai (p)
		x = Ax (p)
		if (x .eq. 0) then
c		    copy into Z
		    Zi (pz) = i
		    Zx (pz) = 1
		    pz = pz + 1
		else
c		    copy into A
		    Ai (pa) = i
		    Ax (pa) = x
		    pa = pa + 1
		endif
20	    continue
30	continue
	Ap (ncol+1) = pa

	return
	end

c-----------------------------------------------------------------------
c RB*prune: discard explicit zero entries
c-----------------------------------------------------------------------
c
c   nrow-by-ncol: size of A
c   Ap, Ai, Ax: matrix with zeros on input, pruned on output
c
c-----------------------------------------------------------------------

	subroutine RBrprune
     $	    (nrow, ncol, Ap, Ai, Ax)
	integer*8
     $	    nrow, ncol, Ap (ncol+1), Ai (*), i, j, p, pa, pz, p1
	double precision Ax (*), x

c	----------------------------------------------------------------
c	prune the matrix
c	----------------------------------------------------------------

	pa = 1
	do 20 j = 1, ncol
c	    save the new start of column j of A
	    p1 = Ap (j)
	    Ap (j) = pa
c	    prune column j of A
	    do 10 p = p1, Ap (j+1)-1
		i = Ai (p)
		x = Ax (p)
		if (x .ne. 0) then
c		    copy into A
		    Ai (pa) = i
		    Ax (pa) = x
		    pa = pa + 1
		endif
10	    continue
20	continue
	Ap (ncol+1) = pa

	return
	end

max_stars_repo_name: andrejadd/HBR-net-inf

max_issues_repo_name: andrejadd/HBR-net-inf

max_forks_repo_name: andrejadd/HBR-net-inf


text:
*DECK GAMI
      FUNCTION GAMI (A, X)
C***BEGIN PROLOGUE  GAMI
C***PURPOSE  Evaluate the incomplete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      SINGLE PRECISION (GAMI-S, DGAMI-D)
C***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate the incomplete gamma function defined by
C
C GAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
C
C GAMI is evaluated for positive values of A and non-negative values
C of X.  A slight deterioration of 2 or 3 digits accuracy will occur
C when GAMI is very large or very small, because logarithmic variables
C are used.  GAMI, A, and X are single precision.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  ALNGAM, GAMIT, 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***END PROLOGUE  GAMI
C***FIRST EXECUTABLE STATEMENT  GAMI
      IF (A .LE. 0.0) CALL XERMSG ('SLATEC', 'GAMI',
     +   'A MUST BE GT ZERO', 1, 2)
      IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMI',
     +   'X MUST BE GE ZERO', 2, 2)
C
      GAMI = 0.0
      IF (X.EQ.0.0) RETURN
C
C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW.
      FACTOR = EXP (ALNGAM(A) + A*LOG(X) )
C
      GAMI = FACTOR * GAMIT(A, X)
C
      RETURN
      END


text: #redirect Big Blue

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*DK ocmidpt
      subroutine ocmidpt(x,x1,x2,nn)
C
C
C#######################################################################
C
C      PURPOSE -
C
C         THIS ROUTINE
C
C      INPUT ARGUMENTS -
C
C         input_message - CHARACTER STRING CONTAINING THE INPUT
C                            MESSAGE.
C
C      OUTPUT ARGUMENTS -
C
C         ierror_return - ERROR RETURN CODE (==0 ==> OK, <>0 ==> ERROR)
C
C      CHANGE HISTORY -
C
C        $Log: ocmidpt.f,v $
C        Revision 2.00  2007/11/05 19:46:03  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.1   Mon Apr 14 16:56:12 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.0   01/17/95 16:39:08   pvcs
CPVCS     Original version
C
C#######################################################################
C
      implicit real*8 (a-h, o-z)
      character*132 logmess
C
C#######################################################################
C
c
c        This routine finds the midpoint x along the shorter great
c        circle arc between points x1 and x2 on the unit sphere.
c
      real x(nn,3), x1(nn,3), x2(nn,3)
c
      do 20 j=1,3
 20   x(1,j) = x1(1,j) + x2(1,j)
      xnorm = 1./sqrt(x(1,1)**2 + x(1,2)**2 + x(1,3)**2)
      do 40 j=1,3
 40   x(1,j) = xnorm*x(1,j)
c
      goto 9999
 9999 continue
      return
      end

max_forks_repo_name: daniellivingston/LaGriT