text:
      program phase
      implicit none

c     # Grid parameters
      integer maxmx, maxmy
      parameter(maxmx = 2**12, maxmy = 2**12)

c     # Model parameters
      double precision xi, Tm, Lv, C, kappa, gma, mu, k,
     &      sigma, tau, beta, S, alpha, m, gma2, beta_p,
     &      tau_th, tanh

c     # Temporary variables
      double precision uij_n, src, d, dpdt, gphi,
     &  pij_n, pij_np1

c     # Input parameters
      double precision ax_in, ay_in, dx_in, dy_in, t_in, tstart
      double precision domain_length
      integer mx_in, my_in, nstp, izero, nchar
      character*100 fname
      logical restart

      double precision phi(maxmx, maxmy)
      double precision u(maxmx, maxmy)
      double precision theta, r, rp, s1, s2, amp

c     # More grid parameters
      integer mx, my, i, j, n, nout, nstep
      double precision ax, ay, bx, by, dx, dy, dt, tfinal, t
      double precision x(maxmx), y(maxmy)

c     # Misc
      integer iframe, iframe_start
      double precision pi, dsgn

      integer mbdcnd, nbdcnd, ierror, idimf
      double precision lambda_u, lambda_phi, pertrb
      double precision bda(maxmy), bdb(maxmy)
      double precision bdc(maxmx), bdd(maxmx)
      double precision f(maxmx, maxmy)

      double precision r0, x0, y0, kanio, pdel
      common /parm_comm/ r0, x0, y0, kanio, pdel

      integer wsize, wsize_t
c         13M + 4N + M*INT(LOG2(N))
c             log2(N) <= log2(maxmy) ~ 10
      parameter(wsize = 4*maxmy + (13 + 12)*maxmx)
      double precision work(wsize)

      double precision xlow, ylow, w

      wsize_t = 4*(maxmy) +
     &      (13 + int(log(real(my))/log(2.d0))*mx)
      if (wsize < wsize_t) then
         write(6,*) 'Increase wsize from ', wsize, ' to ',wsize_t
         stop
      endif
      pi = 4.d0*datan(1.d0)

c     # --------------------------------
c     # Set up grid
c     # --------------------------------

      open(10,file='phase.dat')
      read(10,*) mx, my
c      read(10,*) ax, bx, ay, by
      read(10,*) domain_length
      read(10,*) nstep, nout
      read(10,*) dt

      if (mx > maxmx .or. my > maxmy) then
         write(6,*) 'mx or my too big'
         stop
      endif

      ax = -domain_length/2
      bx = domain_length/2
      ay = -domain_length/2
      by = domain_length/2

      dx = (bx - ax)/mx
      dy = (by - ay)/my

      do i = 1,mx
         x(i) = ax + (i-0.5)*dx
      enddo

      do j = 1,my
         y(j) = ay + (j-0.5)*dy
      enddo

c     # ----------------------------
c     # physical parameters
c     # ----------------------------
      read(10,*) tm
      read(10,*) lv
      read(10,*) C
      read(10,*) kappa
      read(10,*) gma
      read(10,*) mu
      read(10,*) k
      read(10,*) sigma
      read(10,*) tau
      read(10,*) beta

c     # ----------------------------
c     # model parameters
c     # ----------------------------
      read(10,*) s
      read(10,*) alpha
      read(10,*) m
      read(10,*) xi
      read(10,*) kanio
      read(10,*) r0
      read(10,*) amp
      read(10,*) pdel
      read(10,*) gma2
      read(10,*) beta_p
      read(10,*) tau_th

c     # ----------------------------
c     # Get restart information
c     # ----------------------------
      read(10,*) restart
      if (restart) then
         read(10,*) fname
      endif
      close(10)

      write(6,*) 'Parameters used for this run : '
      write(6,*) '-------------------------------'
      write(6,*) 'Physical parameters '
      write(6,*) '-------------------------------'
      write(6,200) 'Tm    = ', tm
      write(6,200) 'Lv    = ', Lv
      write(6,200) 'C     = ', C
      write(6,200) 'kappa = ', kappa
      write(6,210) 'Gamma = ', gma
      write(6,200) 'mu    = ', mu
      write(6,200) 'k     = ', k
      write(6,210) 'sigma = ', sigma
      write(6,210) 'tau   = ', tau
      write(6,200) 'beta  = ', beta
      write(6,*) ' '
      write(6,*) '---------------------------'
      write(6,*) 'Non-dimensional parameters'
      write(6,*) '---------------------------'
      write(6,200) 'S       = ', S
      write(6,200) 'alpha   = ', alpha
      write(6,200) 'm       = ', m
      write(6,200) 'xi      = ', xi
      write(6,200) 'k       = ', kanio
      write(6,200) 'r0      = ', r0
      write(6,200) 'amp     = ', amp
      write(6,200) 'pdel    = ', pdel
      write(6,200) 'Gamma   = ', gma2
      write(6,200) 'beta_p  = ', beta_p
      write(6,200) 'tau_th  = ', tau_th
      write(6,*) ' '
      write(6,*) '---------------------------'
      write(6,*) 'Grid parameters'
      write(6,*) '---------------------------'
      write(6,220) 'mx   = ', mx
      write(6,220) 'my   = ', my
      write(6,200) 'ax   = ', ax
      write(6,200) 'bx   = ', bx
      write(6,200) 'ay   = ', ay
      write(6,200) 'by   = ', by
      write(6,200) 'dx   = ', dx
      write(6,200) 'dy   = ', dy

      if (restart) then
         write(6,*) 'Restart from file ', fname
      else
         write(6,*) 'Initializing run from scratch...'
      endif

 200  format(a,F10.5)
 210  format(a,E12.6)
 220  format(a,I5)

c     # --------------------------------
c     # Initialize phi
c     # --------------------------------
      write(6,*) 'Initializing phi and u...'

      if (restart) then
         open(20,file=fname)
         read(20,*) mx_in
         read(20,*) my_in
         read(20,*) ax_in
         read(20,*) ay_in
         read(20,*) dx_in
         read(20,*) dy_in
         read(20,*) t_in

         if (mx_in .ne. mx .or. my_in .ne. my) then
            write(6,*) 'Problem with restart: mx or my not right.'
            stop
         elseif (ax_in .ne. ax .or. ay_in .ne. ay) then
            write(6,*) 'Problem with restart: ax or ay not right'
            stop
         endif
         izero = iachar('0')
         tstart = t_in
         iframe_start = 0
         nstp = 1
         do i = 10,7,-1
            nchar = iachar(fname(i:i)) - izero
            iframe_start = iframe_start + nchar*nstp
            nstp = nstp*10
         enddo

         do j = 1,my
            do i = 1,mx
               read(20,*) phi(i,j), u(i,j)
            enddo
         enddo
         close(20)
      else

         x0 = (ax + bx)/2.d0
         y0 = (ay + by)/2.d0
         tstart = 0.0
         iframe = 0
         do j = 1,my
            do i = 1,mx
                xlow = ax + (i-1)*dx
                ylow = ay + (j-1)*dy
c                rp = sqrt((xlow)**2 + (ylow)**2)

                call cellave(xlow,ylow,dx,dy,w)
                u(i,j) = w-1.d0
                phi(i,j) = 1-w

c                r = sqrt((amp*(x(i) - x0))**2 + (y(j) - y0)**2)
c               theta = atan2(y(j) - y0,x(i) - x0)
c               rp = r0*(1.d0 + pdel*cos(kanio*theta))
c               d = tanh((r-rp)/xi)
c               u(i,j) = -0.5d0*(d + 1)  !! In [-1,0]
c               phi(i,j) = 0 0.5d0*(d + 1)  !! In [0,1]
            enddo
         enddo
         write(6,*) 'Done with initialization'

c     # Output frame only if we have initialized data from scratch.
         write(6,600) 0,tstart
         call out2(maxmx,maxmy,mx,my,ax,ay,dx,dy,phi,u,tstart,
     &        iframe_start)
      endif

c     % --------------------------------------
c     % Time step through results
c     % --------------------------------------

      do j = 1,my
         bda(j) = 0.d0
         bdb(j) = 0.d0
      enddo

      do i = 1,mx
         bdc(i) = 0.d0
         bdd(i) = 0.d0
      enddo

      lambda_phi = -1.d0/(m*dt)
      lambda_u = -1.d0/dt
      idimf = maxmx
      mbdcnd = 3
      nbdcnd = 3

      iframe = iframe_start
      do n = 1,nout
         t = tstart + n*dt

c         write(6,100) n,t
  100    format('Step',I5,' at time t = ',1PE12.5)

c        # Compute right hand side for phi :
         do j = 1,my
            do i = 1,mx
               pij_n = phi(i,j)
               uij_n = u(i,j)
               gphi = pij_n**2*(1.d0-pij_n)**2
               src = 30.d0*xi*S*alpha*gphi*uij_n +
     &              pij_n*(1.d0 - pij_n)*(pij_n - 0.5d0)
               f(i,j) = lambda_phi*(dt*(m/(xi**2))*src + pij_n)
            enddo
         enddo

c        # Solve equation for phi; solution will be in f
         call hstcrt(ax,bx,mx,mbdcnd,bda,bdb,ay,by,my,nbdcnd,
     &         bdc,bdd,lambda_phi,f,idimf,pertrb,ierror,work)

         if (ierror .ne. 0) then
            write(6,*) 'hstcrt : ierror /= 0; ', ierror
            stop
         endif

         if (work(1) > wsize) then
            write(6,*) 'hstcrt : work(1) > wsize; ', work(1), wsize
            stop
         endif

         if (abs(pertrb) .gt. 0) then
            write(6,*) 'pertrb > 0;  ', pertrb
            stop
         endif

c        # Set up right hand side for u
         do j = 1,my
            do i = 1,mx
               pij_n = phi(i,j) !! old phi(i,j)
               pij_np1 = f(i,j) !! new phi(i,j)
               phi(i,j) = f(i,j)  !! for next time step

               uij_n = u(i,j)
               dpdt = (pij_np1 - pij_n)/dt !! d(phi)/dt
               gphi = pij_np1**2*(1.0 - pij_np1)**2
               src = 30.d0*gphi*(-1.d0/S)*dpdt
               f(i,j) = lambda_u*(dt*src + uij_n)
            enddo
         enddo

c        # Solve equation for u; solution will be in f
         call hstcrt(ax,bx,mx,mbdcnd,bda,bdb,ay,by,my,nbdcnd,
     &        bdc,bdd,lambda_u,f,idimf,pertrb,ierror,work)

         if (ierror .ne. 0) then
            write(6,*) 'hstcrt : ierror /= 0', ierror
            stop
         endif

         if (abs(pertrb) .gt. 0) then
            write(6,*) 'pertrb > 0;  ', pertrb
            stop
         endif

         do j = 1,my
            do i = 1,mx
               u(i,j) = f(i,j)
            enddo
         enddo

c        # ---------------------------------
c        # Output results
c        # ---------------------------------

          write(6,601) n, t

         if (nstep*(n/nstep) == n) then
            iframe = iframe + 1
            write(6,600) iframe, t
            write(6,*) ' '
            call out2(maxmx,maxmy,mx,my,ax,ay,dx,dy,phi,u,t,iframe)
         endif
  600    format('Writing frame ',I5,' at time t = ',1PE12.5)
  601    format('Step ',I5,' at time t = ',1PE12.5)

      enddo

      end

      double precision function fdisc(x,y)
      implicit none

      double precision x,y, r, rp, theta

      double precision r0, x0, y0, kanio, pdel
      common /parm_comm/ r0, x0, y0, kanio, pdel

      r = sqrt((x-x0)**2 + (y-y0)**2)
      theta = atan2(y - y0,x - x0)
      rp = r0*(1.d0 + pdel*cos(kanio*theta))
      fdisc = r - rp

      end

c      double precision function tanh(x)
c      implicit none
c      double precision x
c      tanh = (exp(x) - exp(-x))/(exp(x) + exp(-x))
c      return
c      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
-timescale=1ns/1ps
//-------------------------------------
// DUT Define
//-------------------------------------

//-------------------------------------
// DUT Include
//-------------------------------------

//-------------------------------------
// DUT Filelist
//-------------------------------------

//-------------------------------------
// C Define
//-------------------------------------

//-------------------------------------
// C Include
//-------------------------------------

//-------------------------------------
// C Filelist
//-------------------------------------

//-------------------------------------
// SV Define
//-------------------------------------
+define+YUU_AHB_MAX_MASTER_NUM=1
+define+YUU_AHB_MAX_SLAVE_NUM=1
+define+YUU_AHB_MAX_ADDR_WIDTH=32
+define+YUU_AHB_MAX_DATA_WIDTH=32

//-------------------------------------
// SV Include
//-------------------------------------
+incdir+../../yuu_common/include
+incdir+../../yuu_common/src/sv
+incdir+../include
+incdir+../src/sv
+incdir+../test

//-------------------------------------
// SV Filelist
//-------------------------------------
../../yuu_common/include/yuu_common_pkg.sv
../include/yuu_ahb_pkg.sv

//-------------------------------------
// Case List
//-------------------------------------
../test/yuu_ahb_base_case.sv
../test/yuu_ahb_direct_case.sv
../test/yuu_ahb_ral_case.sv
../test/yuu_ahb_single_base_case.sv
../test/yuu_ahb_single_direct_case.sv

//-------------------------------------
// Top Module
//-------------------------------------
single_top.sv

ext: f


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This file contains interfaces for the maxwell fmm3d and
! direct calculation. It contains the two subroutines:
!
!   emfmm3d - Maxwell fmm3d with nd densities
!
!   em3ddirect - Maxwell direct calculation with nd
!    densities
!
!   (input sources and output requests are controleld by flags)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine emfmm3d(nd,eps,zk,ns,source,ifa_vect,a_vect,&
 ifb_vect,b_vect,iflambda,lambda,nt,targets,ifE,E,ifcurlE,curlE,&
 ifdivE,divE,ier)
implicit none

!jmc_edc

!  This function computes
!
!    E = curl S_{k}[a_vect] + S_{k}[b_vect] + grad S_{k}[lambda]  -- (1)
!  using the vector Helmholtz fmm.
!
!  The subroutine also computes divE, curlE
!  with appropriate flags
!
!  Remark: the subroutine uses a stabilized representation
!    for computing the divergence by using integration by parts
!    wherever possible. If the divergence is not requested, then the
!    helmholtz fmm is called with 3*nd densities, while if the divergence
!    is requested, then the helmholtz fmm is calld with 4*nd densities
!
!
!  input:
!    nd - integer
!      number of densities (a_vect,b_vect,lambda)
!
!    eps - real *8
!      epsilon for the fmm call
!
!    zk - complex *16
!      Helmholtz parameter
!
!    ns - integer
!      number of sources
!
!    source - real *8 (3,ns)
!      location of the sources
!
!    ifa_vect - integer
!     Contribution due to curl S_{k}[a_vect] is included if
!     ifa_vect = 1
!
!    a_vect - complex *16(nd,3,ns)
!      a vector source
!
!    ifb_vect - integer
!      Contribution due to S_{k}[b_vect] is included if
!      ifb_vect = 1
!
!    b_vect - complex *16(nd,3,ns)
!      b vector source
!
!    iflambda - integer
!      Contribution due to grad S_{k} [lambda] is included
!      if iflambda = 1
!
!    lambda - complex *16(nd,ns)
!      lambda source
!
!    nt - integer
!      number of targets
!
!    targets - real *8 (3,nt)
!      location of the targets
!
!    ifE - integer
!      E is returned at the target locations if ifE = 1
!
!    ifcurlE - integer
!      curl E is returned at the target locations if
!      ifcurlE = 1
!
!    ifdivE - integer
!      div E is returned at the target locations if
!      ifdivE = 1
!
!  output:
!
!    E - complex *16(nd,3,nt)
!      E field defined in (1) above
!
!    curlE - complex *16(nd,3,nt)
!      curl of E field
!
!    divE - complex *16(nd,nt)
!      divergence of E
!

  !List of calling arguments
  double precision, intent(in) :: eps
  double complex, intent(in) :: zk
  integer, intent(in) :: nd,ns,nt
  double precision, intent(in) :: source(3,ns)
  integer, intent(in) :: ifa_vect
  double complex, intent(in) :: a_vect(nd,3,*)
  integer, intent(in) :: ifb_vect
  double complex, intent(in) :: b_vect(nd,3,*)
  integer, intent(in) :: iflambda
  double complex, intent(in) :: lambda(nd,*)
  integer, intent(in) :: ifE
  double complex, intent(out) :: E(nd,3,*)
  integer, intent(in) :: ifcurlE
  double complex, intent(out) :: curlE(nd,3,*)
  integer, intent(in) :: ifdivE
  double complex, intent(out) :: divE(nd,*)
  double precision, intent(in) :: targets(3,nt)

  !List of local variables
  double complex, allocatable :: sigma_vect(:,:,:)
  double complex, allocatable :: dipvect_vect(:,:,:,:)
  double complex, allocatable :: Etmp(:,:,:)
  double complex, allocatable :: gradE_vect(:,:,:,:)

  double complex, allocatable :: pot(:),grad(:,:),hess(:,:)
  double complex, allocatable :: hesstarg(:,:)
  integer ndens
  integer ier,iper

  integer i,j,nd0,l,m
  integer ifcharge,ifdipole,ifpgh,ifpghtarg

  ndens = 3
  if(ifdivE.eq.1) ndens = 4

  !!Initialize sources
  allocate(sigma_vect(nd,ndens,ns))
  allocate(dipvect_vect(nd,ndens,3,ns))
  allocate(Etmp(nd,ndens,nt))
  allocate(gradE_vect(nd,ndens,3,nt))

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,l,j,m)
  do i=1,nt
    do l=1,ndens
      do j=1,nd
        Etmp(j,l,i)=0.0d0
      enddo
    enddo

    do l=1,3
      do m=1,ndens
        do j=1,nd
          gradE_vect(j,m,l,i)=0.0d0
        enddo
      enddo
    enddo
  enddo
!$OMP END PARALLEL DO

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,l,j,m)
  do i=1,ns
    do l=1,ndens
      do j=1,nd
        sigma_vect(j,l,i)=0.0d0
      enddo
    enddo

    do l=1,3
      do m=1,ndens
        do j=1,nd
          dipvect_vect(j,m,l,i)=0.0d0
        enddo
      enddo
    enddo
  enddo
!$OMP END PARALLEL DO

  if (ifb_vect.eq.1) then

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,l,j)
    do i=1,ns
      do l=1,3
        do j=1,nd
          sigma_vect(j,l,i)=sigma_vect(j,l,i)+b_vect(j,l,i)
        enddo
      enddo
    enddo
!$OMP END PARALLEL DO

    if(ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,l,j)
      do i=1,ns
        do l=1,3
          do j=1,nd
            dipvect_vect(j,4,l,i) = dipvect_vect(j,4,l,i) - b_vect(j,l,i)
          enddo
        enddo
      enddo
!$OMP END PARALLEL DO
    endif
  endif

  if (iflambda.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,ns
      do j=1,nd
        dipvect_vect(j,1,1,i)=dipvect_vect(j,1,1,i)-lambda(j,i)
        dipvect_vect(j,2,2,i)=dipvect_vect(j,2,2,i)-lambda(j,i)
        dipvect_vect(j,3,3,i)=dipvect_vect(j,3,3,i)-lambda(j,i)
      enddo
    enddo
!$OMP END PARALLEL DO

!
!  Add in the contribution corresponding to the divergence
!  Here we use the identity that
!  \nabla \cdot \nabla S_{k} [\lambda] = k^2 S_{k} [\lambda]
!
    if(ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
      do i=1,ns
        do j=1,nd
          sigma_vect(j,4,i) = sigma_vect(j,4,i) - zk**2*lambda(j,i)
        enddo
      enddo
!$OMP END PARALLEL DO
    endif
  endif

  if (ifa_vect.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,ns
      do j=1,nd
        dipvect_vect(j,1,2,i)=dipvect_vect(j,1,2,i)-a_vect(j,3,i)
        dipvect_vect(j,1,3,i)=dipvect_vect(j,1,3,i)+a_vect(j,2,i)

        dipvect_vect(j,2,1,i)=dipvect_vect(j,2,1,i)+a_vect(j,3,i)
        dipvect_vect(j,2,3,i)=dipvect_vect(j,2,3,i)-a_vect(j,1,i)

        dipvect_vect(j,3,1,i)=dipvect_vect(j,3,1,i)-a_vect(j,2,i)
        dipvect_vect(j,3,2,i)=dipvect_vect(j,3,2,i)+a_vect(j,1,i)
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  ifcharge=1
  ifdipole=1
  ifpgh = 0
  ifpghtarg=2

  ifcharge=1
  ifdipole=1
  ifpgh = 0
  ifpghtarg=2

  if ((ifcurlE.ne.1).and.(ifdivE.ne.1)) then
    ifpghtarg=1
  endif

  if(ifdivE.ne.1) then
    if ((ifa_vect.ne.1).and.(iflambda.ne.1)) then
      ifdipole=0
    endif
    if (ifb_vect.ne.1) then
      ifcharge=0
    endif
  else
    if (ifb_vect.ne.1.and.iflambda.ne.1) then
      ifcharge = 0
    endif
  endif

  nd0 = ndens*nd

  allocate(pot(nd0),grad(nd0,3),hess(nd0,6),hesstarg(nd0,6))

  call hfmm3d(nd0,eps,zk,ns,source,ifcharge,sigma_vect,ifdipole, &
   dipvect_vect,iper,ifpgh,pot,grad,hess,nt,targets,ifpghtarg,Etmp, &
   gradE_vect,hesstarg,ier)

  if (ifE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,nt
      do l=1,3
        do j=1,nd
          E(j,l,i) = Etmp(j,l,i)
        enddo
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  if (ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,nt
      do j=1,nd
        divE(j,i)=Etmp(j,4,i)
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  if (ifcurlE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,nt
      do j=1,nd
        curlE(j,1,i)=gradE_vect(j,3,2,i)-gradE_vect(j,2,3,i)
        curlE(j,2,i)=gradE_vect(j,1,3,i)-gradE_vect(j,3,1,i)
        curlE(j,3,i)=gradE_vect(j,2,1,i)-gradE_vect(j,1,2,i)
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  deallocate(sigma_vect)
  deallocate(dipvect_vect)
  deallocate(gradE_vect)
  deallocate(Etmp)

return
end subroutine emfmm3d
!
!
!
!
!

subroutine em3ddirect(nd,zk,ns,source,ifa_vect,a_vect,&
 ifb_vect,b_vect,iflambda,lambda,nt,targets,ifE,E,ifcurlE,curlE,&
 ifdivE,divE,thresh)
implicit none

!  This function computes
!    E = curl S_{k}[a_vect] + S_{k}[b_vect] + grad S_{k}[lambda]  -- (1)
!  via direct computation.
!  The subroutine also returns divE, curlE with appropriate flags
!
!  input:
!    nd - integer
!      number of densities (a_vect,b_vect,lambda)
!
!    eps - real *8
!      epsilon for the fmm call
!
!    zk - complex *16
!      Helmholtz parameter
!
!    ns - integer
!      number of sources
!
!    source - real *8 (3,ns)
!      location of the sources
!
!    ifa_vect - integer
!     Contribution due to curl S_{k}[a_vect] is included if
!     ifa_vect = 1
!
!    a_vect - complex *16(nd,3,ns)
!      a vector source
!
!    ifb_vect - integer
!      Contribution due to S_{k}[b_vect] is included if
!      ifb_vect = 1
!
!    b_vect - complex *16(nd,3,ns)
!      b vector source
!
!    iflambda - integer
!      Contribution due to grad S_{k} [lambda] is included
!      if iflambda = 1
!
!    lambda - complex *16(nd,ns)
!      lambda source
!
!    nt - integer
!      number of targets
!
!    targets - real *8 (3,nt)
!      location of the targets
!
!    ifE - integer
!      E is returned at the target locations if ifE = 1
!
!    ifcurlE - integer
!      curl E is returned at the target locations if
!      ifcurlE = 1
!
!    ifdivE - integer
!      div E is returned at the target locations if
!      ifdivE = 1
!
!
!    thresh - real *8
!      threshold for self interaction term
!
!
!  output:
!
!    E - complex *16(nd,3,nt)
!      E field
!
!    curlE - complex *16(nd,3,nt)
!      curlE curl of E field
!
!    divE - complex *16(nd,nt)
!      divE divergence of E
!

  !List of calling arguments
  double complex, intent(in) :: zk
  integer, intent(in) :: nd,ns,nt
  double precision, intent(in) :: source(3,ns)
  integer, intent(in) :: ifa_vect
  double complex, intent(in) :: a_vect(nd,3,*)
  integer, intent(in) :: ifb_vect
  double complex, intent(in) :: b_vect(nd,3,*)
  integer, intent(in) :: iflambda
  double complex, intent(in) :: lambda(nd,*)
  integer, intent(in) :: ifE
  double complex, intent(out) :: E(nd,3,*)
  integer, intent(in) :: ifcurlE
  double complex, intent(out) :: curlE(nd,3,*)
  integer, intent(in) :: ifdivE
  double complex, intent(out) :: divE(nd,*)
  double precision, intent(in) :: targets(3,nt)
  double precision, intent(in) :: thresh

  !List of local variables
  double complex, allocatable :: sigma_vect(:,:,:)
  double complex, allocatable :: dipvect_vect(:,:,:,:)
  double complex, allocatable :: gradE_vect(:,:,:,:)
  double complex, allocatable :: Etmp(:,:,:)
  integer i,j,nd0,l,m
  integer ifcharge,ifdipole,ifpot,ifgrad,ifpghtarg
  integer ndens

  !!Initialize sources

  ndens = 3
  if(ifdivE.eq.1) ndens = 4

  allocate(sigma_vect(nd,ndens,ns))
  allocate(dipvect_vect(nd,ndens,3,ns))
  allocate(Etmp(nd,ndens,nt))
  allocate(gradE_vect(nd,ndens,3,nt))

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,m)
  do i=1,nt
    do l=1,ndens
      do j=1,nd
        Etmp(j,l,i)=0.0d0
      enddo
    enddo

    do l=1,3
      do m=1,ndens
        do j=1,nd
          gradE_vect(j,m,l,i) = 0.0d0
        enddo
      enddo
    enddo
  enddo
!$OMP END PARALLEL DO

!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,m)
  do i=1,ns
    do l=1,ndens
      do j=1,nd
        sigma_vect(j,l,i)=0.0d0
      enddo
    enddo

    do l=1,3
      do m=1,ndens
        do j=1,nd
          dipvect_vect(j,m,l,i)=0.0d0
        enddo
      enddo
    enddo
  enddo
!$OMP END PARALLEL DO

  if (ifb_vect.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l)
    do i=1,ns
      do l=1,3
        do j=1,nd
          sigma_vect(j,l,i)=sigma_vect(j,l,i)+b_vect(j,l,i)
        enddo
      enddo
    enddo
!$OMP END PARALLEL DO

    if(ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l)
      do i=1,ns
        do l=1,3
          do j=1,nd
            dipvect_vect(j,4,l,i) = dipvect_vect(j,4,l,i) - b_vect(j,l,i)
          enddo
        enddo
      enddo
!$OMP END PARALLEL DO
    endif
  endif

  if (iflambda.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,ns
      do j=1,nd
        dipvect_vect(j,1,1,i)=dipvect_vect(j,1,1,i)-lambda(j,i)
        dipvect_vect(j,2,2,i)=dipvect_vect(j,2,2,i)-lambda(j,i)
        dipvect_vect(j,3,3,i)=dipvect_vect(j,3,3,i)-lambda(j,i)
      enddo
    enddo
!$OMP END PARALLEL DO

!
!  Add in the contribution corresponding to the divergence
!  Here we use the identity that
!  \nabla \cdot \nabla S_{k} [\lambda] = k^2 S_{k} [\lambda]
!
    if(ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
      do i=1,ns
        do j=1,nd
          sigma_vect(j,4,i) = sigma_vect(j,4,i) - zk**2*lambda(j,i)
        enddo
      enddo
!$OMP END PARALLEL DO
    endif
  endif

  if (ifa_vect.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,ns
      do j=1,nd
        dipvect_vect(j,1,2,i)=dipvect_vect(j,1,2,i)-a_vect(j,3,i)
        dipvect_vect(j,1,3,i)=dipvect_vect(j,1,3,i)+a_vect(j,2,i)

        dipvect_vect(j,2,1,i)=dipvect_vect(j,2,1,i)+a_vect(j,3,i)
        dipvect_vect(j,2,3,i)=dipvect_vect(j,2,3,i)-a_vect(j,1,i)

        dipvect_vect(j,3,1,i)=dipvect_vect(j,3,1,i)-a_vect(j,2,i)
        dipvect_vect(j,3,2,i)=dipvect_vect(j,3,2,i)+a_vect(j,1,i)
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  ifcharge=1
  ifdipole=1
  ifpghtarg=2

  if ((ifcurlE.ne.1).and.(ifdivE.ne.1)) then
    ifpghtarg=1
  endif

  if(ifdivE.ne.1) then
    if ((ifa_vect.ne.1).and.(iflambda.ne.1)) then
      ifdipole=0
    endif
    if (ifb_vect.ne.1) then
      ifcharge=0
    endif
  else
    if (ifb_vect.ne.1.and.iflambda.ne.1) then
      ifcharge = 0
    endif
  endif

  nd0=ndens*nd
  if(ifpghtarg.eq.1) then
    if(ifcharge.eq.1.and.ifdipole.eq.0) then
      call h3ddirectcp(nd0,zk,source,sigma_vect,ns, &
         targets,nt,Etmp,thresh)
    endif
    if(ifcharge.eq.0.and.ifdipole.eq.1) then
      call h3ddirectdp(nd0,zk,source,dipvect_vect,ns, &
         targets,nt,Etmp,thresh)
    endif
    if(ifcharge.eq.1.and.ifdipole.eq.1) then
      call h3ddirectcdp(nd0,zk,source,sigma_vect,dipvect_vect,ns, &
         targets,nt,Etmp,thresh)
    endif

  else if(ifpghtarg.eq.2) then
    if(ifcharge.eq.1.and.ifdipole.eq.0) then
      call h3ddirectcg(nd0,zk,source,sigma_vect,ns, &
         targets,nt,Etmp,gradE_vect,thresh)
    endif
    if(ifcharge.eq.0.and.ifdipole.eq.1) then
      call h3ddirectdg(nd0,zk,source,dipvect_vect,ns, &
         targets,nt,Etmp,gradE_vect,thresh)
    endif
    if(ifcharge.eq.1.and.ifdipole.eq.1) then
      call h3ddirectcdg(nd0,zk,source,sigma_vect,dipvect_vect,ns, &
         targets,nt,Etmp,gradE_vect,thresh)
    endif
  endif

  if(ifE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l)
    do i=1,nt
      do l=1,3
        do j=1,nd
          E(j,l,i) = E(j,l,i) + Etmp(j,l,i)
        enddo
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  if (ifdivE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,nt
      do j=1,nd
        divE(j,i)=divE(j,i) + Etmp(j,4,i)
       enddo
    enddo
!$OMP END PARALLEL DO
  endif

  if (ifcurlE.eq.1) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
    do i=1,nt
      do j=1,nd
        curlE(j,1,i)=curlE(j,1,i)+gradE_vect(j,3,2,i)-gradE_vect(j,2,3,i)
        curlE(j,2,i)=curlE(j,2,i)+gradE_vect(j,1,3,i)-gradE_vect(j,3,1,i)
        curlE(j,3,i)=curlE(j,3,i)+gradE_vect(j,2,1,i)-gradE_vect(j,1,2,i)
      enddo
    enddo
!$OMP END PARALLEL DO
  endif

  deallocate(Etmp)
  deallocate(sigma_vect)
  deallocate(dipvect_vect)
  deallocate(gradE_vect)

return
end subroutine em3ddirect


text:
program useClass

! Object-oriented programming (OOP) is a programming paradigm based on
! the concept of "objects", which can contain data and code: data in the
! form of fields (often known as attributes or properties), and code,
! in the form of procedures (often known as methods).

use iso_fortran_env
use particles_objs

  implicit none
  integer :: ii
  character(4) :: nn
  character(1) :: c
  real(real64) :: r, m, q, x, carga
  type(Particle) :: p
  type(Electron) :: eneg, epos
  type(Electron), dimension(4) :: W

  r = 2.4e-10
  m = 9.109e-31
  q = 1.602e-19

  x = 4.0_real64
  p = Particle("stone", r, m)
  eneg = Electron("negatron", r, m, -q)
  epos = Electron("positron", r, m, q)

  do ii = 1, 4
    write(c,'(I1)') ii
    nn = "w"//c
    carga = ((-1.0)**ii) * q/(real(ii))
    W(ii) = Electron(nn, r/real(ii), m*real(ii)/2.0, carga)
    call w(ii)%show()
  end do

  !W(1) = Electron("w1", r/1.0, m*2, q/4.0)
  !W(2) = Electron("w2", r/2.0, m*2, q/3.0)
  !W(3) = Electron("w3", r/3.0, m*2, q/2.0)
  !W(4) = Electron("w4", r/4.0, m*2, q/1.0)

  print *, p%radius, p%mass
  print *, eneg%radius, eneg%mass, eneg% charge
  print *, epos%radius, epos%mass, epos% charge

  call p%cross_area()
  call eneg%cross_area()
  call epos%cross_area()

  print *, eneg%electric_potential(x)
  print *, epos%electric_potential(x)

print *, "-------------------------------------------"
print*, " "
  call p%show()
  call eneg%show()
  call epos%show()

  do ii = 1, 4
    call w(ii)%show()
  end do

end program useClass

max_stars_repo_name: EdgardoBonzi/Fortran-Examples

max_issues_repo_name: EdgardoBonzi/Fortran-Examples

max_forks_repo_name: EdgardoBonzi/Fortran-Examples


text:
	SUBROUTINE GH_KGLB ( timstr, type, sname, sped, adnm, wind,
     +                       alat, alon, iflags, rmnlat, rmnlon, rmxlat,
     +                       rmxlon, tzone, iwht, iblk, idisp, namstr,
     +                       poslb, xkey, ykey, nknt, wdir, jflags,
     +			     disstr, zone2, iret )
C************************************************************************
C* GH_KGLB								*
C*									*
C* This subroutine draws the legend box (key) and plots the storm       *
C* type/name string.  It also creates and plots the date, time and day, *
C* advisory number, current  location, maximum wind speed, current      *
C* movement, marker keys, and track label strings (if any) for the      *
C* legend box, and plots the logos.                                     *
C*									*
C* GH_KGLB ( TIMSTR, TYPE, SNAME, SPED, ADNM, WIND, ALAT, ALON, IFLAGS, *
C*	     RMNLAT, RMNLON, RMXLAT, RMXLON, TZONE, IWHT, IBLK, IDISP,  *
C*           NAMSTR, POSLB, XKEY, YKEY, NKNT, WDIR, JFLAGS, DISSTR,     *
C*	     ZONE2, IRET )						*
C*									*
C* Input parameters:							*
C*	TIMSTR		CHAR*		Ending time for watches		*
C*	TYPE		CHAR*		Storm type			*
C*	SNAME		CHAR*		Storm name			*
C*	SPED		CHAR*		Current movement in kts		*
C*	ADNM		CHAR*		Storm advisory number		*
C*	WIND		CHAR*		Max sustained Wind in kts	*
C*	ALAT (*)	REAL		Current latitude		*
C*	ALON (*)	REAL		Current longitude		*
C*      IFLAGS(4)	INTEGER		Flags for storm watch/warning   *
C*	RMNLAT 		REAL		Minimum y axis map value	*
C*	RMNLON 		REAL		Minimum x axis map value	*
C*	RMXLAT 		REAL		Maximum y axis map value	*
C*	RMXLON 		REAL		Maximum x axis map value	*
C*      TZONE		CHAR*		Storm's current time zone	*
C*	IWHT		INTEGER  	White color value		*
C*	IBLK		INTEGER		Black color value		*
C*	IDISP		INTEGER		Storm dissipation flag		*
C*	NAMSTR		CHAR*		Storm type/name string          *
C*	POSLB		CHAR*		Flag for label/logo position    *
C*	XKEY (*)	REAL		Norm. x coords of legend box    *
C*	YKEY (*)	REAL		Norm. y coords of legend box    *
C*      NKNT		INTEGER         Number of forecast points	*
C*      WDIR		CHAR*           Direction of present movement   *
C*      JFLAGS(3)	INTEGER		Flags for forecast intensity    *
C*									*
C* Output parameters:							*
C*	DISSTR		CHAR*		Storm dissipation time string	*
C*	ZONE2		CHAR*		Local time zone                 *
C*	IRET		INTEGER		Return code			*
C*									*
C**									*
C* Log:									*
C* A. Hardy/GSC		 2/01   					*
C* A. Hardy/GSC		 6/01	Added diss. check and label pos. parms  *
C* A. Hardy/GSC		 6/01	Added color tag; GH_SAVE, GH_REST	*
C* A. Hardy/SAIC	 8/01	Added TD type;widen box;add text,lines  *
C* A. Hardy/SAIC	10/01	Added current lat/lon to GH_KGIP	*
C* D. Kidwell/NCEP	 4/02	Removed code for label box placement;   *
C*			        default marker color 20->17             *
C* A. Hardy/NCEP	 8/02   Set 115 kts wind to 135 mph		*
C* J. Wu/SAIC		 8/03   Added legend for day 4-5 track area	*
C* D. Kidwell/NCEP	11/03   Added output argument zone2             *
C* D. Kidwell/NCEP	 1/04   Added wdir, 'TPC' line & intensity lines*
C* D. Kidwell/NCEP	 4/04   CSC - added jflags processing           *
C* S. Jacobs/NCEP	11/10	Removed "TPC" from the label string	*
C* M. Sardi/NHC	        10/02   Set 115 kts wind to 130 mph		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	PARAMETER	( NUMH2 = 6 )
C*
	CHARACTER*(*)	timstr, type, sname, sped, wind, tzone, disstr,
     +                  adnm, namstr, poslb, zone2, wdir
        REAL		alat(*), alon(*), xkey(*), ykey(*)
        INTEGER		iflags(*), jflags (*)
C*
        CHARACTER       dtstr*40, strtm*20, clat*5, clon*6,
     +			nslab*1, ewlab*1, curloc*48, maxwin*35,
     +			curmov*40, currlc*24, forpos*28, pttkar*28,
     +			smwd*4, cmov*3, coltag*33, hscolr(4)*12,
     +                  advstr*12, pttkar5*28, attstr*50, cmpdir*3,
     +			fintns(3)*27
      	CHARACTER       fmonth*9, fmon*9, fday*9, cyear*4, cdate*2,
     +                  ctime*2, ampm*2
	INTEGER		jtarr(5), lenstp(4), ihs(4)
        CHARACTER       flgstr*240, stype(4)*36
	LOGICAL 	baddir
	REAL            xx(5), yy(5), xlin(6), ylin(6), xl(9), yl(9)
C*
	DATA		stype / 'Hurricane Warning',
     +                          'Hurricane Watch',
     +                          'Tropical Storm Warning' ,
     +                          'Tropical Storm Watch' /
	DATA		lenstp  / 17, 15, 22, 20 /
	DATA		hscolr  / 'w_hur_warn', 'w_hur_watch',
     +                            'w_ts_warn', 'w_ts_watch' /
	DATA		ihs     / 2, 11, 4, 20 /
	DATA		attstr  / 'NWS National Hurricane Center' /
	DATA		fintns  / 'H  Sustained wind > 73 mph ',
     +				  'S  Sustained wind 39-73 mph',
     +				  'D  Sustained wind < 39 mph ' /
C-----------------------------------------------------------------------
      	iret = 0
C
C*	Query attributes.
C
        CALL GH_SAVE ( ier )
C
C*	Create the date and time strings.
C
 	CALL GH_TIME ( timstr, tzone, jtarr, hours, zone2,
     +                 fmonth, fmon, fday, ampm, ier )
C
C*	Create date string.
C
	CALL ST_INCH ( jtarr(1), cyear, ier )
	CALL ST_INCH ( jtarr(3), cdate, ier )
	CALL ST_INCH ( jtarr(4), ctime, ier )
	CALL ST_LSTR( fmonth, lenm, ier )
	CALL ST_LSTR( cyear, leny, ier )
	CALL ST_LSTR( cdate, lend, ier )
	CALL ST_LSTR( ctime, lent, ier )
	CALL ST_LSTR( ampm, lena, ier )
	CALL ST_LSTR( zone2, lenz, ier )
C
	dtstr = fmonth(:lenm) //' '// cdate(:lend) //', '// cyear(:leny)
	strtm = ctime(:lent) // ' ' // ampm(:lena) // ' '
     +          // zone2(:lenz) // ' '// fday
        disstr = ctime(:lent) // ' ' // ampm(:lena) // ' '//
     +		 zone2(:lenz) // ' ' // fmonth(:lenm) //' '//
     +		 cdate(:lend)
C
        IF ( ( idisp .eq. 0 ) .or. ( idisp .eq. 2 ) ) THEN
C
C*	    Create the advisory number string.
C
	    CALL ST_LSTR( adnm, lenv, ier )
            advstr = 'Advisory '// adnm (:lenv)
C
C*	    Create the current location string.
C
	    IF ( alat(1) .gt. 0.0) THEN
                nslab = 'N'
              ELSE IF ( alat (1) .lt. 0.0 ) THEN
                nslab = 'S'
              ELSE
                nslab = ' '
            END IF
C
	    IF ( alon(1) .lt. 0.0) THEN
                ewlab = 'W'
              ELSE IF ( alon (1) .gt. 0.0 ) THEN
                ewlab = 'E'
              ELSE
                ewlab = ' '
            END IF
            CALL ST_RLCH ( alat(1), 1, clat, ier )
            CALL ST_RLCH ( alon(1), 1, clon, ier )
            IF ( clat(1:1) .eq. '-') clat = clat(2:)
            IF ( clon(1:1) .eq. '-') clon = clon(2:)
            CALL ST_LSTR (clat, leng, ier )
            CALL ST_LSTR (clon, lens, ier )
            CALL ST_LSTR (ewlab, lenw, ier )
            curloc = 'Current Center Location '//
     +               clat(:leng)// ' ' // nslab //
     +                '  '//clon(:lens) // ' '//ewlab(:lenw)
C
C*	    Create the max sustained wind string. Round wind speed
C*          to the 5 mph. If the wind speed is equal to 115 knots then
C*          it will be automatically set to 130 mph.
C
            CALL ST_LSTR ( wind, lens, ier )
            CALL ST_CRNM ( wind(:lens), fwnd, ier )
            wnd = PR_KNMH ( fwnd )
            IF ( ( wnd .gt. 131.18 ) .and. ( wnd .lt. 133.0 ) ) THEN
                imswnd = 130
              ELSE
                mswnd = NINT ( wnd )
                imswnd = ( mswnd / 5 ) * 5
                IF ( (MOD (mswnd, 5) .ge. 2.5 ) ) imswnd  = imswnd + 5
            END IF
            CALL ST_INCH ( imswnd ,  smwd, ier )
            CALL ST_LSTR ( smwd, lenx, ier )
	    maxwin = 'Max Sustained Wind  ' // smwd(:lenx) // ' mph'
C
C*	    Create the current movement string.
C
	    baddir = .true.
	    IF ( wdir .ne. ' ' ) THEN
	        CALL ST_CRNM ( wdir, dirmov, ier )
	        CALL CLO_COMPASS ( dirmov, cmpdir, leni, ier )
	        IF ( ier .eq. 0 ) baddir = .false.
	    END IF
	    IF ( baddir ) THEN
	        cmpdir = ' '
	        leni   = 1
	    END IF
C
            CALL ST_CRNM ( sped, fspd, ier )
            spd = PR_KNMH ( fspd )
            mspd = NINT ( spd )
            CALL ST_INCH ( mspd , cmov, ier )
            CALL ST_LSTR ( cmov, lens, ier )
            IF (cmov(:lens) .ne. '0' ) THEN
	        curmov = 'Current Movement ' // cmpdir (:leni) //
     +			 ' at  ' // cmov(:lens) // ' mph'
              ELSE IF (cmov(:lens) .eq. '0' ) THEN
	        curmov = 'Current Movement  Stationary'
            END IF
C
C*	    Create the current location, forecast positions and
C*          potential track area strings.
C
	    currlc = 'Current Center Location'
	    forpos = 'Forecast Center Positions'
	    pttkar = 'Potential Day 1-3 Track Area'
	    pttkar5 = 'Potential Day 4-5 Track Area'
C
	    CALL ST_LSTR ( attstr, lenr, ier )
            CALL ST_LSTR ( namstr, lenn, ier )
            CALL ST_LSTR ( dtstr,  lend, ier )
            CALL ST_LSTR ( strtm,  lens, ier )
            CALL ST_LSTR ( advstr, lenv, ier )
            CALL ST_LSTR ( curloc, lenc, ier )
            CALL ST_LSTR ( maxwin, lenm, ier )
            CALL ST_LSTR ( curmov, lenu, ier )
            CALL ST_LSTR ( currlc, lenl, ier )
            CALL ST_LSTR ( forpos, leno, ier )
            CALL ST_LSTR ( pttkar, lenp, ier )
C
	    IF ( nknt .gt. NUMH2 ) THEN
	        CALL ST_LSTR ( pttkar5, lenq, ier )
	    END IF
C
	    CALL GQSYSZ ( rxszmk,ryszmk,rxsztx,rysztx,rxszwb,ryszwb,ier)
C
	    CALL GSCOLR ( iwht, ier )
            CALL GSSMTH ( 0, 1.5, ier )
	    CALL GSFILL ( 1.0, 1, ier )
	    CALL GFILL ( 'N', 5, xkey, ykey, ier )
	    CALL GSCOLR ( iblk, ier )
            CALL GSLINE ( 1, 1, 2, 1, ier )
            CALL GLINE ( 'N', 5, xkey, ykey, ier )
C
	    CALL GSCOLR ( iblk, ier )
C
            CALL GSTEXT ( 22, 2, 1.2, 1, 111, 1, 1, ier )
C
C*          Plot storm name line.
C
            xxz   = xkey (1) + .01
            yyz   = ykey (1) - .02
            ixoff = 0
            CALL GTEXT ( 'N', xxz, yyz, namstr(:lenn), 0.0, ixoff,
     +                   0, ier )
C
C*          Plot date line.
C
            yyz = yyz - .022
            CALL GTEXT ( 'N', xxz, yyz, dtstr(:lend), 0.0, ixoff,
     +                   0, ier )
C
C*          Plot time line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, strtm(:lens), 0.0, ixoff,
     +                   0, ier )
C
            CALL GSTEXT ( 22, 2, 1.1, 1, 111, 1, 1, ier )
	    CALL GQSYSZ ( rxszmk,ryszmk,rxsztx,rysztx,rxszwb,ryszwb,ier)
C
C*          Plot attribution line
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, attstr(:lenr), 0.0, ixoff,
     +                   0, ier )
C
C*          Plot advisory number line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, advstr(:lenv), 0.0, ixoff,
     +               	 0, ier )
C
C*          Plot lat. lon. line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, curloc(:lenc), 0.0, ixoff,
     +               	 0, ier )
C
C*          Plot wind line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, maxwin(:lenm), 0.0, ixoff,
     +                   0, ier )
C
C*          Plot movement line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, curmov(:lenu), 0.0, ixoff,
     +                   0, ier )
C
C*          Plot current location line.
C
            yyz   = yyz - .02
            ixoff = 12
            CALL GTEXT ( 'N', xxz, yyz, currlc(:lenl), 0.0, ixoff,
     +                   0, ier )
            sym = 3.
            xzz = xxz + .03
	    CALL GSCOLR ( iblk, ier )
	    CALL GSSPCL ( 1.1, 1, ier )
	    CALL GSPCL ( 'N', 1, sym, xzz, yyz, 0, 0, ier )
	    coltag = 'w_lbl_curr_pos_mkr'
            CALL ST_LSTR ( coltag, lens, ier )
 	    CALL GH_COLR ( coltag(:lens), 17, ier )
	    CALL GSSPCL ( .9, 1, ier )
	    CALL GSPCL ( 'N', 1, sym, xzz, yyz, 0, 0, ier )
            CALL GSCOLR ( iblk, ier )
            CALL GSMRKR ( 17, 0, 0.8, 2, ier )
            CALL GMARK ( 'N', 1, xzz, yyz, ier )
	    CALL GSCOLR ( iblk, ier )
C
C*          Plot forecast position line.
C
            yyz = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, forpos(:leno), 0.0, ixoff,
     +                   0, ier )
            CALL GSCOLR ( iblk, ier )
            CALL GSMRKR ( 17, 0, 1.4, 2, ier )
            CALL GMARK ( 'N', 1, xzz, yyz, ier )
	    CALL GSCOLR ( iblk, ier )
C
C*          Plot forecast intensity lines.
C
	    ixoff = 15
	    DO ii = 1, 3
		IF ( jflags ( ii ) .gt. 0 ) THEN
                    yyz = yyz - .02
                    CALL GTEXT ( 'N', xxz, yyz, fintns (ii) (:27), 0.0,
     +			         ixoff, 0, ier )
		END IF
	    END DO
C
C*          Plot day 1-3 potential track area line.
C
            ixoff = 12
            yyz   = yyz - .02
            CALL GTEXT ( 'N', xxz, yyz, pttkar(:lenp), 0.0, ixoff,
     +                   0, ier )
            xlin(1) = xxz + .05
            xlin(2) = xlin(1) - .033
            xlin(3) = xxz
            xlin(4) = xxz
            xlin(5) = xlin(4) + .010
            xlin(6) = xlin(1)
C
            ylin(1) = yyz - .006
            ylin(2) = ylin(1)
            ylin(3) = ylin(2) + .005
            ylin(4) = ylin(3) + .01
            ylin(5) = ylin(4) + .005
            ylin(6) = ylin(1)
C
            CALL GSLINE ( 1, 1, 2, 1, ier )
            CALL GLINE ( 'N', 6, xlin, ylin, ier )
C
C*          Plot and fill day 4-5 potential track area line.
C
	    IF ( nknt .gt. NUMH2 ) THEN
		yyz = yyz - .02
                CALL GTEXT ( 'N', xxz, yyz, pttkar5(:lenq), 0.0, ixoff,
     +                       0, ier )
                xl(1) = xxz + .04
                xl(2) = xl(1) - .025
                xl(3) = xxz + 0.005
                xl(4) = xl(3)
                xl(5) = xl(4) + .01
                xl(6) = xl(5) + .025
                xl(7) = xl(6) - .01
                xl(8) = xl(7)
                xl(9) = xl(1)
C
                yl(1) = yyz - .008
                yl(2) = yl(1)
                yl(3) = yl(2) + .004
                yl(4) = yl(3) + .008
                yl(5) = yl(4) + .004
                yl(6) = yl(5)
                yl(7) = yl(6) - .006
                yl(8) = yl(7) - .004
                yl(9) = yl(1)
C
                CALL GQFILL ( szfil, iftyp, ier )
C
                CALL GSFILL ( 1.0, 2, ier )
                CALL GSLINE ( 1, 1, 2, 1, ier )
                CALL GLINE ( 'N', 9, xl, yl, ier )
     	        CALL GFILL ( 'N', 9, xl, yl, ier )
C
                CALL GSFILL ( szfil, iftyp, ier )
C
	    END IF
C
C*	    Plot active hurricane /T.S. watches and warnings text lines.
C
            DO is = 1,4
 	        IF ( iflags (is) .ne. 0 ) THEN
                    flgstr =   stype(is)(:lenstp(is))
                    CALL ST_LSTR (flgstr, lens, ier )
                    yyz = yyz - .02
                    CALL GTEXT ( 'N', xxz, yyz, flgstr(:lens), 0.0,
     +                            ixoff, 0, ier )
C
C*	            Plot color legend for tropical storm.
C
                    xx (1) = xxz +.04
                    xx ( 2 ) = xx ( 1 ) - .033
	            xx ( 3 ) = xx ( 2 )
	            xx ( 4 ) = xx ( 1 )
	            xx ( 5 ) = xx ( 1 )
 	            yy ( 1 ) = yyz -.008
	            yy ( 2 ) = yy ( 1 )
                    IF ( is .lt. 3 ) yy ( 3 ) = yy ( 1 ) + .015
                    IF ( is .ge. 3 ) yy ( 3 ) = yy ( 1 ) + .010
	            yy ( 4 ) = yy ( 3 )
	            yy ( 5 ) = yy ( 1 )
C
                    coltag = hscolr (is)
                    CALL ST_LSTR ( coltag, lens, ier )
                    CALL GH_COLR ( coltag(:lens), ihs(is), ier )
C
	            CALL GSFILL ( 1.0, 1, ier )
	            CALL GFILL ( 'N', 5, xx, yy, ier )
	            CALL GSCOLR ( iblk, ier )
                END IF
            END DO
        END IF
C
C*      Plot logos opposite of text box label
C
        IF ( poslb .eq. 'B' ) THEN
             xnwx = rmnlon + .09
             ynwx = rmxlat - .28
             xnoaa = rmxlon + .13
             ynoaa = rmxlat - .28

          ELSE
             xnwx = rmnlon + .09
             ynwx = rmnlat + .07
             xnoaa = rmxlon + .13
             ynoaa = rmnlat + .07
        END IF
C
C*      Plot NWS logo.
C
	CALL GLOGO ( 'N', xnwx, ynwx, 2.5, 2, 2, ier )
C
C*      Plot NOAA logo.
C
        CALL GLOGO ( 'N', xnoaa, ynoaa, 2.5, 2, 1, ier )
C
C*	Restore attributes.
C
        CALL GH_REST ( ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


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

      SUBROUTINE LCOLOR (COLOR)
C***********************************************************************

C  SUBROUTINE LCOLOR = SETS THE LINE COLOR

C***********************************************************************

      CHARACTER*5 COLOR

      IF (COLOR .EQ. 'WHITE') THEN
         CALL PLTSTD(1, 7.)
      ELSEIF (COLOR .EQ. 'BLACK') THEN
         CALL PLTSTD(1, 0.)
      ELSEIF (COLOR .EQ. 'YELOW') THEN
         CALL PLTSTD(1, 3.)
      ELSEIF (COLOR .EQ. 'PINK ') THEN
         CALL PLTSTD(1, 5.)
      ELSEIF (COLOR .EQ. 'RED  ') THEN
         CALL PLTSTD(1, 1.)
      ELSEIF (COLOR .EQ. 'BLUE ') THEN
         CALL PLTSTD(1, 4.)
      ENDIF
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
module Mod_supm_TemperatureCoupling
   use Mod_supm_BaseElmope
   use Mod_supm_ExternalForces
   use Mod_nsi_BoussinesqForces
   use typre
   implicit none
   private
   public :: SetPointersTemperatureCoupling
   real(rp), allocatable :: eltem(:)
   real(rp) :: gptem(1), lambdaT
   integer(ip), allocatable :: kfl_IsSet

contains

   subroutine SetPointersTemperatureCoupling(itask)
      implicit none
      integer(ip) :: itask
      select case (itask)

      case(0)
         allocate(kfl_IsSet)
         call a%Memor%allocObj(0,'kfl_IsSet','InitProcedurePointer',1)
         kfl_IsSet = -1

      case(1)
         if (kfl_IsSet == -1) then
            !Pointers are already set
            kfl_IsSet = 1

            if (a%kfl_cotem == 1) then
               call SetPointersExternalForces(1)
               call ConcatenateProcedures(ProcHook_Initializations,AllocTem)
               call ConcatenateProcedures(ProcHook_Gathers,GatherTem)
               call ConcatenateProcedures(ProcHook_Interpolates,InterpolateTem)
               if (a%kfl_ExternalTemperatureSGS == 1) then
                  call ConcatenateProcedures(ProcHook_Interpolates,InterpolateTemSGS)
               endif

               if (a%kfl_cotem_Boussinesq == 1) call ConcatenateProcedures(ProcPointer%ExternalForces_sup,BoussinesqForces)
               if (a%kfl_cotem_WLF==1)          call ConcatenateProcedures(ProcHook_Physicalprop,ViscoelasticTempParametersWLFmodel)
               if (a%kfl_cotem_Arrhenius==1)    call ConcatenateProcedures(ProcHook_Physicalprop,ViscoelasticTempParameters_ArrheniusModel)

               if (a%kfl_cotem_WLF==1 .or. a%kfl_cotem_Arrhenius==1) then
                  !Update physical parameters
                  call ConcatenateProcedures(ProcHook_Physicalprop,UpdateViscoelasticAuxiliarParameter)
                  if (a%LogFormulation==1)call ConcatenateProcedures(ProcHook_Physicalprop,UpdateLogarithmicAuxiliarParameters)
                  !Print viscosity and lambda
                  if (a%npp_stepi(5)>=1)  call ConcatenateProcedures(ProcHook_Physicalprop,PlotViscosity)
                  if (a%npp_stepi(23)>=1) call ConcatenateProcedures(ProcHook_Physicalprop,PlotLambda)
               end if

               call ConcatenateProcedures(ProcHook_Finalizations,DeallocTem)
            endif
         endif

      case(100)
         deallocate(kfl_IsSet)
         call a%Memor%deallocObj(0,'kfl_IsSet','InitProcedurePointer',1)

      end select
   end subroutine

   !-------------------------------------------------------------------
   !Computation Subroutines
   subroutine AllocTem
      implicit none
      call a%Memor%alloc(e%mnode,eltem,'eltem','nsm_elmope')
   end subroutine

   subroutine DeallocTem
      implicit none
      call a%Memor%dealloc(e%mnode,eltem,'eltem','nsm_elmope')
   end subroutine

   subroutine GatherTem
      implicit none
      call e%gather(1,eltem,a%tempe)
   end subroutine

   subroutine InterpolateTem
      implicit none
      call e%interpg(1,eltem,gptem(1))

   end subroutine

   subroutine InterpolateTemSGS
      implicit none
      gptem = gptem + a%tesgs(ielem)%a(1,e%igaus)
   end subroutine

   subroutine BoussinesqForces
      implicit none
      call nsi_BoussinesqForces(e,acden,a%bougr,a%boube,a%boutr,gptem(1),elext)
   end subroutine

   subroutine ViscoelasticTempParametersWLFmodel
      implicit none
      call sup_templaw_WLF(e, gptem, a%ReferenceTemp, a%c1_WLF, a%c2_WLF, a%MatProp(imat)%LawViParam, acvis, lambdaT)
   end subroutine

   subroutine ViscoelasticTempParameters_ArrheniusModel
      implicit none
      call sup_templaw_Arrhenius(e, gptem, a%ReferenceTemp, a%alpha_Arrhenius, a%MatProp(imat)%LawViParam, acvis, lambdaT)
   end subroutine

   subroutine UpdateViscoelasticAuxiliarParameter
      implicit none
      call a%IncrementalLambda(imat,lambdaT,lambda)
      auxVE  = lambda/(2.0_rp*acvis)
   end subroutine

   subroutine UpdateLogarithmicAuxiliarParameters
      implicit none
      if(a%kfl_timei==0) then
        lambda0=lambda*a%MatProp(1)%LawviParam(5)
      else
        lambda0=a%MatProp(imat)%LawViParam(3)*a%MatProp(1)%LawviParam(5)
      end if

      auxL = ((1.0_rp-beta)*acvis)/lambda0
      auxG   = a%MatProp(imat)%LawviParam(4)*lambda/(2.0_rp*(lambda0**2_ip) + 0.00001_rp) !GIESEKUSLCR
   end subroutine

   subroutine PlotViscosity
      implicit none
      a%viscarray(ielem)%a(e%igaus) = acvis
   end subroutine

    subroutine PlotLambda
      implicit none
      a%lambdarray(ielem)%a(e%igaus) = lambda
    end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM