text:
PROGRAM porefem
!-------------------------------------------------------------------------
! Program porefem
!
!       Three-dimensional strain of an elastic solid using
!       8-, 14- or 20-node brick hexahedra. Mesh numbered in x-y
!       planes then in the z-direction. No global stiffness matrix
!       assembly. Diagonally preconditioned conjugate gradient solver.
!-------------------------------------------------------------------------
 USE main
 USE geom
 USE gaf95
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15),npri=1
 INTEGER::cg_iters,cg_limit,fixed_freedoms,i,iel,k,loaded_nodes,ndim=3,   &
   ndof,nels,neq,nip,nn,nprops=2,np_types,nod,nodof=3,nr,nst=6,nxe,nye,   &
   nze,nlen,nlen_out,ii,iii,s
 INTEGER::ip,iq,maxfld,j,ipore
 REAL(iwp)::alpha,beta,cg_tol,det,one=1.0_iwp,penalty=1.0e20_iwp,up,      &
   zero=0.0_iwp,xl,yl,zl,dx1,dy1,dz1,small=1.0e-6_iwp,dtim=1.0_iwp
 REAL(iwp)::ee,ei,vi,meanqe,sdqe,meanqv,meanqv1,sdqv,sdqv1,pernom,porosity,ll,ul, &
            three=3.0_iwp,two=2.0_iwp,meanpo,sdpo,start_time,end_time,    &
            pcg_start_time,pcg_end_time
 CHARACTER(LEN=15)   :: element='hexahedron',argv,trash
 CHARACTER(LEN=15)   :: filename_out, format_string
 LOGICAL::cg_converged,solid=.true.
 LOGICAL::debug2 = .false.
!-----------------------dynamic arrays------------------------------------
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),nf(:,:),no(:),    &
   node(:),num(:),sense(:),imap(:)
 REAL(iwp),ALLOCATABLE::bee(:,:),coord(:,:),d(:),dee(:,:),der(:,:),       &
   deriv(:,:),diag_precon(:),eld(:),fun(:),gc(:),g_coord(:,:),jac(:,:),   &
   km(:,:),loads(:),p(:),points(:,:),prop(:,:),sigma(:),store(:),         &
   storkm(:,:,:),u(:),value(:),weights(:),x(:),xnew(:),x_coords(:),       &
   y_coords(:),z_coords(:),val(:,:),temp1(:),temp2(:)
 REAL(iwp),ALLOCATABLE::enom(:),qe(:),qv(:),perc(:),qv1(:),kmi(:,:)
!-----------------------input and initialisation--------------------------
 integer ieplt(3)
 real, allocatable:: cfld(:),phifld(:),psifld(:),gamfld(:),efld(:),vfld(:)
 real  c(7),phi(7),psi(7),gam(7),e(7),v(7),R(6,6)
 integer kseed,istat,nsim,ns,nseed,iterm,nrfx,nrfy,nrfz,nfld,jfld,ifld
 real dx,dy,dz,thx,thy,thz,ti,ts
 logical verbos,debug,lxfld,dmpfld
 character sub2*128, job*80, sub1*80,varfnc*6
 common/dbgrfl/ istat, debug
 common/lastym/ ti, ts
 istat = 13
 iterm = 14
 dmpfld=.false.
 debug=.false.
 verbos=.false.
 lxfld=.false.
 sub1='none'
 sub2='none'
 job='equivalent_3D'
 varfnc='dlavx3'
 nfld=1
 jfld=2
 ifld=19
 call flush(6)
 CALL getname(argv,nlen)
 OPEN(10,FILE=argv(1:nlen)//'.dat')
 OPEN(11,FILE=argv(1:nlen)//'.res')
 CALL cpu_time(start_time)
 READ(10,*)nod,nip
 READ(10,*)xl,yl,zl
 READ(10,*)nxe,nye,nze
 READ(10,*)ei,vi
 READ(10,*)cg_tol,cg_limit
 READ(10,*)trash
 read(10,*)c
 read(10,*)phi
 read(10,*)psi
 read(10,*)gam
 read(10,*)e
 read(10,*)v
 READ(10,*)trash
 READ(10,*)r
 READ(10,*)thx,thy,thz
 READ(10,*)kseed,nsim
 READ(10,*)porosity
 kseed = iseed(kseed)
 WRITE(11,'(A)')job
 WRITE(11,'(A,I2,A)')'Analysis uses', nod,' noded elements'
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Overall element dimensions'
 WRITE(11,'(A,F5.2)')'xl =',xl
 WRITE(11,'(A,F5.2)')'yl =',yl
 WRITE(11,'(A,F5.2)')'zl =',zl
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Number of elements used'
 WRITE(11,'(A,I3)')'x-direction =',nxe
 WRITE(11,'(A,I3)')'y-direction =',nye
 WRITE(11,'(A,I3)')'z-direction =',nze
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Stiffness parameters of intact rock'
 WRITE(11,'(A,F8.2)')"Young's modulus =",ei
 WRITE(11,'(A,F5.2)')"Poisson's ratio =",vi
 WRITE(11,*)
! WRITE(11,'(A)')"c phi psi gam e v prop erties of embankment (mean, SD, dist, L, U, m, s)"
 WRITE(11,'(A)')"Normal distribution (voids)"
 write(11,*)c
! write(11,*)phi
! write(11,*)psi
! write(11,*)gam
! write(11,*)e
! write(11,*)v
! WRITE(11,*)
! WRITE(11,'(A)')"Correlation matrix"
! WRITE(11,'(6f8.2)')r
! WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Spatial correlation length of voids'
 WRITE(11,'(A,F8.2)')"x-direction =",thx
 WRITE(11,'(A,F8.2)')"y-direction =",thy
 WRITE(11,'(A,F8.2)')"z-direction =",thz
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Starting seed'
 WRITE(11,'(I5)')kseed
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Number of simulations'
 WRITE(11,'(I5)')nsim
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Target porosity'
 WRITE(11,'(F5.2)')porosity
 ll=dphinv(0.5_iwp-porosity*0.5_iwp)
 ul=-ll
 !  dphinv(x) gives the abscissa that has the area x to its left in a normal distibution
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Limits of porosity zone in standard space'
 WRITE(11,'(A,F7.4)')"lower =",ll
 WRITE(11,'(A,F7.4)')"upper =",ul
 CALL chkne3(istat,iterm,verbos,nxe,nye,nze,nrfx,nrfy,nrfz)
 maxfld=9*nrfx*nrfy*nrfz/8
 WRITE(11,'(A,I2,A)')'nrfx nrfy nrfz'
 WRITE(11,*)nrfx,nrfy,nrfz
 WRITE(11,*)
 CALL mesh_size(element,nod,nels,nn,nxe,nye,nze)
 ndof=nod*nodof

! Allocation of large arrays that dominate memory usage

 ALLOCATE(nf(nodof,nn))
 ALLOCATE(g_coord(ndim,nn))
 ALLOCATE(g_num(nod,nels),g_g(ndof,nels))
 ALLOCATE(imap(nels),etype(nels))

 ! ALLOCATE(enom(nxe*nye*nze)) ! - LM - allocated but not used

! Allocation of small arrays

 ALLOCATE(weights(nip),num(nod))
 ALLOCATE(points(nip,ndim),dee(nst,nst),coord(nod,ndim))
 ALLOCATE(jac(ndim,ndim),der(ndim,nod),deriv(ndim,nod),fun(nod),gc(ndim))
 ALLOCATE(bee(nst,ndof),km(ndof,ndof),eld(ndof),sigma(nst))
 ALLOCATE(g(ndof),x_coords(nxe+1),y_coords(nye+1),z_coords(nze+1),kmi(ndof,ndof))
 ALLOCATE(storkm(ndof,ndof,2))

! ALLOCATE(g(ndof),y_coords(nye+1),z_coords(nze+1),storkm(ndof,ndof,nels))
! ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),coord(nod,ndim),     &
!   jac(ndim,ndim),der(ndim,nod),deriv(ndim,nod),fun(nod),gc(ndim),        &
!   bee(nst,ndof),km(ndof,ndof),eld(ndof),sigma(nst),g_coord(ndim,nn),     &
!   g_num(nod,nels),weights(nip),num(nod),g_g(ndof,nels),x_coords(nxe+1),  &
!   g(ndof),y_coords(nye+1),z_coords(nze+1),storkm(ndof,ndof,nels),        &
!   cfld(maxfld),phifld(maxfld),psifld(maxfld),gamfld(maxfld),efld(maxfld),&
!   vfld(maxfld),enom(nxe*nye*nze),imap(nels),etype(nels))

 x_coords=zero
 y_coords=zero
 z_coords=zero
 dx1=xl/dble(nxe)
 do i=2,nxe+1
 x_coords(i)=x_coords(i-1)+dx1
 enddo
 dy1=yl/dble(nye)
 do i=2,nye+1
 y_coords(i)=y_coords(i-1)+dy1
 enddo
 dz1=yl/dble(nze)
 do i=2,nze+1
 z_coords(i)=z_coords(i-1)-dz1
 enddo
 dx=dx1*dble(nxe)/dble(nrfx)
 dy=dy1*dble(nye)/dble(nrfy)
 dz=dz1*dble(nze)/dble(nrfz)
! call exca(dx,dy,thx,thy,varfnc,vf)
!
! prop(1)=c, prop(2)=phi, prop(3)=psi, prop(4)=gamma, prop(5)=e, prop(6)=v
!
 call propst(c)
 call propst(phi)
 call propst(psi)
 call propst(gam)
 call propst(e)
 call propst(v)
 IF(SUM(ABS(r))>6)lxfld=.TRUE.
 nf=1
  !-----------------------loop the elements to find global arrays sizes-----
 DO iel=1,nels
   CALL hexahedron_xz(iel,x_coords,y_coords,z_coords,coord,num)
   g_coord(:,num)=TRANSPOSE(coord)
 END DO
do i=1,nn
  if(abs(zl+g_coord(3,i))<=small)nf(3,i)=0
  if(abs(g_coord(1,i))<=small)nf(1,i)=0
  if(abs(g_coord(2,i)-yl)<=small)nf(2,i)=0
  if(abs(g_coord(3,i))<=small)nf(3,i)=0
  if(abs(g_coord(2,i))<=small)nf(2,i)=0
  if(abs(g_coord(1,i)-xl)<=small)nf(1,i)=0
enddo
 nf(2,1)=1
 nf(3,1)=1
 nf(1,nxe+1)=1
CALL formnf(nf)
do i=1,nn
 if(abs(g_coord(2,i))<=small)nf(2,i)=nf(2,1)
 if(abs(g_coord(3,i))<=small)nf(3,i)=nf(3,1)
 if(abs(g_coord(1,i)-xl)<=small)nf(1,i)=nf(1,nxe+1)
enddo
 neq=MAXVAL(nf)
 loaded_nodes=(nxe+1)*(nye+1)
 ALLOCATE(no(loaded_nodes),val(loaded_nodes,ndim))
 ALLOCATE(loads(0:neq),qe(nsim),qv(nsim),qv1(nsim),perc(nsim))
 WRITE(11,'(A,I5,A)')" There are",neq," equations"

!----------element stiffness integration, storage and preconditioner------
elements_1: DO iel=1,nels
  CALL hexahedron_xz(iel,x_coords,y_coords,z_coords,coord,num)
  CALL num_to_g(num,nf,g)
  g_num(:,iel)=num
  g_coord(:,num)=TRANSPOSE(coord)
  g_g(:,iel)=g
 END DO elements_1
 loads=zero
!-----------------------element stiffness integration and assembly--------
 CALL sample(element,points,weights)

   CALL deemat(dee,ei,vi)
   num=g_num(:,1)
   coord=TRANSPOSE(g_coord(:,num))
   kmi=zero
   eld=zero
   gauss_pts_1: DO i=1,nip
     CALL shape_der(der,points,i)
     jac=MATMUL(der,coord)
     det=determinant(jac)
     CALL invert(jac)
     deriv=MATMUL(jac,der)
     CALL beemat(bee,deriv)
     kmi=kmi+MATMUL(matmul(transpose(bee),dee),bee)*det*weights(i)
   END DO gauss_pts_1

   storkm(:,:,2) = kmi          ! medium
   storkm(:,:,1) = kmi/100._iwp ! void

   DEALLOCATE(g_num,g_coord)

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

 do ns =1, nsim

 ALLOCATE(cfld(maxfld))
!ALLOCATE(phifld(maxfld),psifld(maxfld),gamfld(maxfld))
!ALLOCATE(efld(maxfld),vfld(maxfld))

! write(*,*)ns
 if( ns/=1)nseed = iseed(kseed + ns - 1)
!CALL simpl3(istat,iterm,verbos,cfld,phifld,psifld,gamfld,efld,           &
!                   vfld,c,phi,psi,gam,e,v,R,lxfld,thx,thy,thz,nrfx,      &
!                      nrfy,nrfz,dx,dy,dz,dmpfld,nfld,jfld,ieplt,         &
!                      ifld,job,sub1,sub2,varfnc,kseed,debug)

 CALL simpl3cfld(istat,iterm,verbos,cfld,c,phi,psi,gam,e,v,R,lxfld,thx,   &
                       thy,thz,nrfx,                                      &
                       nrfy,nrfz,dx,dy,dz,dmpfld,nfld,jfld,ieplt,         &
                       ifld,job,sub1,sub2,varfnc,kseed,debug)

!DEALLOCATE(phifld,psifld,gamfld,efld,vfld)

 ! phifld, psifld, gamfld, efld & vfld are not used in the main program
 ! perhaps these should be allocated and deallocated in subroutine simpl3

      k = 0
      do iii=1, nze
        do s = 1, nxe
          do i = 1, nye
            k = k + 1
            imap(k) = s + (nrfy-i)*nrfx+(nrfz-iii)*nrfx*nrfy
          enddo
        enddo
      enddo

 ALLOCATE(diag_precon(0:neq))
 diag_precon=zero

 j=0
 elements_2: DO iel=1,nels
   g=g_g(:,iel)
   ii=imap(iel)
   if(cfld(ii)>=ll.AND.cfld(ii)<=ul)then
    etype(iel)=1    !  type 1   (voids)
    km=kmi/100.0_iwp
    j=j+1
   else
    etype(iel)=2    !  type 2   (rock)
    km=kmi
   endif
   DO k=1,ndof
     diag_precon(g(k))=diag_precon(g(k))+km(k,k)
   END DO
 END DO elements_2

 IF(debug2) CALL sortvoid(etype)

 DEALLOCATE(cfld) ! not used again in the program until the next iteration

 ALLOCATE(p(0:neq))
 ALLOCATE(x(0:neq))
 ALLOCATE(xnew(0:neq))
 ALLOCATE(u(0:neq))
 ALLOCATE(d(0:neq))
 ALLOCATE(temp1(ndof),temp2(ndof))

 !-----------------------invert the preconditioner and get starting loads--
 perc(ns)= dble(j)/dble(nxe*nye*nze)
 loads=zero
 loads(nf(3,1))=xl*yl
 loads(0)=0.0
 diag_precon(1:)=one/diag_precon(1:)
 diag_precon(0)=zero
 d=diag_precon*loads
 p=d
!-----------------------pcg equation solution-----------------------------
 x=zero
 cg_iters=0
 CALL cpu_time(pcg_start_time)
 pcg: DO
   cg_iters=cg_iters+1
   WRITE(*,*)"cg_iters = ", cg_iters
   u=zero
   elements_3: DO iel=1,nels
     g    = g_g(:,iel)
     u(g) = u(g)+MATMUL(storkm(:,:,etype(iel)),p(g))
   END DO elements_3
   up    = DOT_PRODUCT(loads,d)
   alpha = up/DOT_PRODUCT(p,u)
   xnew  = x+p*alpha
   loads = loads-u*alpha
   d     = diag_precon*loads
   beta  = DOT_PRODUCT(loads,d)/up
   p     = d+p*beta
   CALL checon(xnew,x,cg_tol,cg_converged)
   IF(cg_converged.OR.cg_iters==cg_limit)EXIT
 END DO pcg
 CALL cpu_time(pcg_end_time)
 xnew(0)=zero
 !write(11,*)xnew(nf(3,1)),loads(nf(3,1)),xl
 qe(ns)=(xl*yl)/(xl*abs(xnew(nf(3,1))))
 qv(ns)=abs(xnew(nf(2,1)))/abs(xnew(nf(3,1)))
 qv1(ns)=abs(xnew(nf(1,nxe+1)))/abs(xnew(nf(3,1)))
 write(*,*)ns,cg_iters
 write(*,*) "Time in PCG = ", pcg_end_time - pcg_start_time
 write(11,'(2I4)')ns,cg_iters

   DEALLOCATE(p,x,xnew,u,temp1,temp2,diag_precon,d)

!  IF(ns==1) THEN

     IF(ns < 10) THEN
       format_string = "(2A,I1)"
     ELSE IF(ns < 100) THEN
       format_string = "(2A,I2)"
     ELSE IF(ns < 1000) THEN
       format_string = "(2A,I3)"
     ELSE IF(ns < 10000) THEN
       format_string = "(2A,I4)"
     ELSE
       WRITE(*,*) "Limit of 9999 realisations in format_string"
       STOP
     END IF

     WRITE(filename_out,format_string) trim(argv),"_r",ns
     nlen_out = len(trim(filename_out))

     ALLOCATE(g_coord(ndim,nn),g_num(nod,nels))
     DO iel=1,nels
       CALL hexahedron_xz(iel,x_coords,y_coords,z_coords,coord,num)
       ! CALL num_to_g(num,nf,g)
       g_num(:,iel)=num
       g_coord(:,num)=TRANSPOSE(coord)
     END DO
!    CALL mesh_ensi(argv,nlen,g_coord,g_num,element,etype,                    &
!                   loads(1:),1,npri,dtim,solid)
     CALL mesh_ensi(filename_out,nlen_out,g_coord,g_num,element,etype,nf,     &
                    loads(1:),1,npri,dtim,solid)
     DEALLOCATE(g_coord,g_num)
!  END IF

  ENDDO
 !!End simulation
 meanqe=sum(qe)/nsim
 sdqe=sqrt(dot_product(qe-meanqe,qe-meanqe)/(nsim-1))
 meanqv=sum(qv)/nsim
 sdqv=sqrt(dot_product(qv-meanqv,qv-meanqv)/(nsim-1))
 meanqv1=sum(qv1)/nsim
 sdqv1=sqrt(dot_product(qv1-meanqv1,qv1-meanqv1)/(nsim-1))
 meanpo=sum(perc)/nsim
 sdpo=sqrt(dot_product(perc-meanpo,perc-meanpo)/(nsim-1))
 write(11,*)
 WRITE(11,'(A)')"Simulation   Young's Modulus  Poisson's ratio(y)   Poisson's ratio(x)  Porosity"
 write(11,'(I5,8X,F10.4,8X,F10.3,10X,F10.3,5X,f10.3)')(i,qe(i),qv(i),qv1(i),perc(i),i=1,nsim)
 WRITE(11,*)
 WRITE(11,'(A,8X,F10.4,8X,F10.3,10X,F10.3,5X,f10.3)')'Means',meanqe,meanqv,meanqv1,meanpo
 WRITE(11,'(A,8X,F10.4,8X,F10.3,10X,F10.3,5X,f10.3)')'SDs  ',sdqe,sdqv,sdqv1,sdpo
! write(11,*)
! WRITE(11,'(A,I5)')" Number of cg iterations to convergence was",cg_iters
! WRITE(11,'(/A)')"  Node   x-disp      y-disp      z-disp"
! DO k=1,nn
!   WRITE(11,'(I5,3f12.4)')k,xnew(nf(:,k))
! END DO
CALL cpu_time(end_time)
WRITE(11,'(/,A,f12.4)')"time taken=",end_time-start_time
WRITE(*,*)"Time taken = ", end_time-start_time
STOP

CONTAINS

SUBROUTINE sortvoid(etype)

INTEGER,INTENT(INOUT) :: etype(:)
INTEGER               :: iel, nels, npores, count
REAL(iwp)             :: start_sort, end_sort
LOGICAL               :: debug = .false.
LOGICAL               :: found = .false.

CALL cpu_time(start_sort)

nels = UBOUND(etype,1)

IF(debug) THEN
  DO iel = 1,nels
    WRITE(*,*) "Element ", iel, " = type ", etype(iel)
  END DO
END IF

ipore = 0
DO iel = 1,nels
  IF(etype(iel)==1) npores = npores + 1
END DO

DO j = 1, npores
  count = 0
  IF(found) THEN
    found = .false.
  END IF
  DO iel = 1, nels
    count = count + 1
    IF(etype(iel)==1) THEN
      found = .true.
      IF(debug) WRITE(*,*) "Position = ",count
      CYCLE
    END IF
  END DO
END DO

IF(debug) THEN
  WRITE(*,*) "Number of pore elements  = ",npores
  WRITE(*,*) "Number of solid elements = ",nels-npores
END IF

CALL cpu_time(end_sort)
WRITE(*,*)"Time to sort voids = ", end_sort - start_sort

END SUBROUTINE sortvoid

END PROGRAM porefem

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
A.J. “Bud” Henle is a retired Teachers in the Davis Joint Unified School District teacher and coach of numerous sports at Davis Senior High School. Coach Henle is the fourth generation of his family to grow up in Davis and Graduates of the Davis Joint Unified School District graduated from DHS in 1942.

During his years as a student at Davis High he was an outstanding athlete in football, basketball and track. His long jump record stood at Davis High for 35 years. Henle is a decorated Veterans veteran of World War II, having served with the 102nd medical battalion of the 27th infantry division in the Pacific Theatre where he earned the Bronze Star and Combat Medal Badge for Gallantry and Service. After the war, he attended UC Davis where he played quarterback for the football team. He transferred to College of the Pacific in Stockton where he received his degree in 1950 and his teaching credential in 1952. Although he began his teaching career elsewhere, he returned to Davis in 1956 at the request of his former coach, Dewey Halden. Initially upon his return, he coached Varsity Football, Basketball, JV and Varsity Baseball while he also taught social studies and physical education. Over the next 29 years, Coach Henle served as athletic director (22 years), coached football (29 years), baseball (15 years), and track (10 years), served as the Block “D” advisor (16 years) and meet director of the Halden Relays (9 years) as he continued to teach in the classroom.

Coach Henle is a member of the Davis Senior High School Hall of Fame, the SacJoaquin Section Football Coaches Association Hall of Fame, the Northern California Old Timers Sports Association Hall of Fame and the La Salle Coaches Hall of Fame.

 Links
  http://www.dhsblueandwhite.org/hofmembers.html DHS Hall of Fame Members

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! ewald_module.f90
! r-space and k-space parts of Ewald sum for ions
MODULE ewald_module

  !------------------------------------------------------------------------------------------------!
  ! This software was written in 2016/17                                                           !
  ! by Michael P. Allen <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk>                        !
  ! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("the authors"),                             !
  ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"),     !
  ! published by Oxford University Press ("the publishers").                                       !
  !                                                                                !
  ! LICENCE                                                                                !
  ! Creative Commons CC0 Public Domain Dedication.                                                 !
  ! To the extent possible under law, the authors have dedicated all copyright and related         !
  ! and neighboring rights to this software to the PUBLIC domain worldwide.                        !
  ! This software is distributed without any warranty.                                             !
  ! You should have received a copy of the CC0 Public Domain Dedication along with this software.  !
  ! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.                               !
  !                                                                                !
  ! DISCLAIMER                                                                                !
  ! The authors and publishers make no warranties about the software, and disclaim liability       !
  ! for all uses of the software, to the fullest extent permitted by applicable law.               !
  ! The authors and publishers do not recommend use of this software for any purpose.              !
  ! It is made freely available, solely to clarify points made in the text. When using or citing   !
  ! the software, you should not imply endorsement by the authors or publishers.                   !
  !------------------------------------------------------------------------------------------------!

  ! References:
  ! Woodcock and Singer, Trans. Faraday Soc. 67, 12, 1971.
  ! de Leeuw et al., Proc. Roy. Soc. A 373, 27, 1980.
  ! Heyes, J. Chem. Phys. 74, 1924, 1981.
  ! see also Fincham, mdions, CCP5 program library.

  ! The self term is subtracted from the k-space contribution
  ! The surface term for simulations in vacuum is not included
  ! A cubic box and unit box length are assumed throughout
  ! No special lists are used

  USE, INTRINSIC :: iso_fortran_env, ONLY : error_unit

  IMPLICIT NONE
  PRIVATE

  ! Public routines
  PUBLIC :: pot_r_ewald, pot_k_ewald, pot_k_pm_ewald

  ! Consider k-vectors within a sphere bounded by the cube (-nk,+nk) in each component
  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: kfac     ! k-dept quantities (k_sq_max)
  INTEGER,                         SAVE :: k_sq_max ! =nk**2

CONTAINS

  FUNCTION pot_r_ewald ( n, r, q, kappa ) RESULT ( pot )
    IMPLICIT NONE
    REAL                              :: pot   ! Returns r-space part of potential energy
    INTEGER,              INTENT(in)  :: n     ! Number of atoms
    REAL, DIMENSION(3,n), INTENT(in)  :: r     ! Positions
    REAL, DIMENSION(n),   INTENT(in)  :: q     ! Charges
    REAL,                 INTENT(in)  :: kappa ! Ewald width parameter

    INTEGER            :: i, j
    REAL, DIMENSION(3) :: rij
    REAL               :: rij_mag, vij

    pot = 0.0

    DO i = 1, n - 1 ! Begin outer loop
       DO j = i + 1, n ! Begin inner loop

          rij(:)  = r(:,i) - r(:,j)           ! Separation vector (box=1 units)
          rij(:)  = rij(:) - ANINT ( rij(:) ) ! Periodic boundaries
          rij_mag = SQRT ( SUM ( rij**2 ) )   ! Magnitude of separation

          vij = q(i) * q(j) * ERFC ( kappa * rij_mag ) / rij_mag ! Screened Coulomb term
          pot = pot + vij

       END DO ! End inner loop
    END DO ! End outer loop

  END FUNCTION pot_r_ewald

  FUNCTION pot_k_ewald ( nk, n, r, q, kappa ) RESULT ( pot )
    IMPLICIT NONE
    REAL                              :: pot   ! Returns k-space part of potential energy
    INTEGER,              INTENT(in)  :: nk    ! Determines number of wave-vectors
    INTEGER,              INTENT(in)  :: n     ! Number of atoms
    REAL, DIMENSION(3,n), INTENT(in)  :: r     ! Positions
    REAL, DIMENSION(n),   INTENT(in)  :: q     ! Charges
    REAL,                 INTENT(in)  :: kappa ! Ewald width parameter

    INTEGER                      :: k, kx, ky, kz, k_sq
    REAL                         :: b, factor, kr_sq
    COMPLEX, DIMENSION(n,  0:nk) :: eikx
    COMPLEX, DIMENSION(n,-nk:nk) :: eiky, eikz
    COMPLEX                      :: term

    LOGICAL, SAVE :: first_call = .TRUE.

    REAL, PARAMETER :: pi = 4.0*ATAN(1.0), twopi = 2.0*pi, twopi_sq = twopi**2

    IF ( first_call ) THEN ! Precalculation of expressions at the start

       b = 1.0 / 4.0 / kappa / kappa
       k_sq_max = nk**2            ! Store this in module data
       ALLOCATE ( kfac(k_sq_max) ) ! Allocate module data array

       ! Lazy triple loop, which over-writes the same values of ksq several times
       ! and leaves some locations empty (those which are not the sum of three squared integers)
       ! We are only interested in the value of k**2, so skip all negative components
       DO kx = 0, nk
          DO ky = 0, nk
             DO kz = 0, nk

                k_sq = kx**2 + ky**2 + kz**2

                IF ( ( k_sq <= k_sq_max ) .AND. ( k_sq /= 0 ) ) THEN ! Test to ensure within range
                   kr_sq      = twopi_sq * REAL ( k_sq )           ! k**2 in real units
                   kfac(k_sq) = twopi * EXP ( -b * kr_sq ) / kr_sq ! Stored expression for later use
                END IF ! End test to ensure within range

             END DO
          END DO
       END DO
       ! End lazy triple loop

       first_call = .FALSE.

    END IF ! End of precalculation of expressions at the start

    ! Double-check value on later calls
    IF ( k_sq_max /= nk**2 ) THEN
       WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'nk error ', k_sq_max, nk**2
       STOP 'Error in pot_k_ewald'
    END IF

    ! Construct exp(ik.r) for all ions and k-vectors

    ! Calculate kx, ky, kz = 0, 1 explicitly
    eikx(:,0) = (1.0, 0.0)
    eiky(:,0) = (1.0, 0.0)
    eikz(:,0) = (1.0, 0.0)

    eikx(:,1) = CMPLX ( COS ( twopi * r(1,:) ), SIN ( twopi * r(1,:) ) )
    eiky(:,1) = CMPLX ( COS ( twopi * r(2,:) ), SIN ( twopi * r(2,:) ) )
    eikz(:,1) = CMPLX ( COS ( twopi * r(3,:) ), SIN ( twopi * r(3,:) ) )

    ! Calculate remaining positive kx, ky and kz by recurrence
    DO k = 2, nk
       eikx(:,k) = eikx(:,k-1) * eikx(:,1)
       eiky(:,k) = eiky(:,k-1) * eiky(:,1)
       eikz(:,k) = eikz(:,k-1) * eikz(:,1)
    END DO

    ! Negative k values are complex conjugates of positive ones
    ! We do not need negative values of kx
    eiky(:,-nk:-1) = CONJG ( eiky(:,nk:1:-1) )
    eikz(:,-nk:-1) = CONJG ( eikz(:,nk:1:-1) )

    pot = 0.0

    DO kx = 0, nk ! Begin outer loop over non-negative kx

       IF ( kx == 0 ) THEN
          factor = 1.0
       ELSE
          factor = 2.0 ! Accounts for skipping negative kx
       END IF

       DO ky = -nk, nk ! Begin inner loop over all ky

          DO kz = -nk, nk ! Begin inner loop over all kz

             k_sq = kx**2 + ky**2 + kz**2

             IF ( ( k_sq <= k_sq_max ) .AND. ( k_sq /= 0 ) ) THEN ! Test to ensure within range

                term = SUM ( q(:) * eikx(:,kx) * eiky(:,ky) * eikz(:,kz) ) ! Sum over all ions
                pot  = pot + factor * kfac(k_sq) * REAL ( CONJG ( term ) * term )

             END IF ! End test to ensure within range

          END DO ! End inner loop over all kz

       END DO ! End inner loop over all ky

    END DO ! End outer loop over non-negative kx

    ! Subtract self part of k-space sum
    pot = pot - kappa * SUM ( q(:)**2 ) / SQRT(pi)

  END FUNCTION pot_k_ewald

  FUNCTION pot_k_pm_ewald ( nk, n, r, q, kappa ) RESULT ( pot )
    USE, INTRINSIC :: iso_c_binding
    USE mesh_module, ONLY : mesh_function, sharpen

    IMPLICIT NONE
    INCLUDE 'fftw3.f03'
    REAL                              :: pot   ! Returns k-space part of potential energy
    INTEGER,              INTENT(in)  :: nk    ! Determines number of wave-vectors
    INTEGER,              INTENT(in)  :: n     ! Number of atoms
    REAL, DIMENSION(3,n), INTENT(in)  :: r     ! Positions
    REAL, DIMENSION(n),   INTENT(in)  :: q     ! Charges
    REAL,                 INTENT(in)  :: kappa ! Ewald width parameter

    INTEGER(C_INT)                                           :: sc      ! Number of points for FFT
    COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:,:,:), ALLOCATABLE :: fft_inp ! Data to be transformed (0:sc-1,0:sc-1,0:sc-1)
    COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:,:,:), ALLOCATABLE :: fft_out ! Output data (0:sc-1,0:sc-1,0:sc-1)
    TYPE(C_PTR)                                              :: fft_plan! Plan needed for FFTW

    INTEGER                             :: ix, iy, iz
    REAL                                :: kx, ky, kz, fx, fxy, fxyz, g, k_sq, dr
    real, DIMENSION(:,:,:), ALLOCATABLE :: rho_sq ! Square modulus of charge density (0:sc-1,0:sc-1,0:sc-1)
    real, dimension(:),     allocatable :: kmesh  ! k-values in wraparound convention (0:sc-1)

    REAL, PARAMETER    :: pi = 4.0*ATAN(1.0), twopi = 2.0*pi, fourpi = 4.0*pi, dk = twopi

    dr = 1.0 / REAL(2*nk) ! r-spacing

    ! Use nk to determine mesh dimension
    sc = 2*nk
    ALLOCATE ( fft_inp(0:sc-1,0:sc-1,0:sc-1), fft_out(0:sc-1,0:sc-1,0:sc-1) )
    allocate ( rho_sq(0:sc-1,0:sc-1,0:sc-1), kmesh(0:sc-1) )

    ! Assign charge density to complex array ready for FFT
    ! Assume r in unit box with range (-0.5,0.5)
    ! Mesh function expects coordinates in range 0..1 so we add 0.5
    fft_inp = CMPLX ( mesh_function ( r+0.5, q, 2*nk ), 0.0 )
    fft_plan = fftw_plan_dft_3d ( sc, sc, sc, fft_inp, fft_out, FFTW_FORWARD, FFTW_ESTIMATE) ! Set up plan for the FFT
    CALL fftw_execute_dft  ( fft_plan, fft_inp, fft_out )                                    ! Execute FFT
    CALL fftw_destroy_plan ( fft_plan )                                                      ! Release plan

    fft_out = fft_out / REAL(sc**3)           ! Incorporate scaling by number of grid points
    rho_sq  = real ( fft_out*conjg(fft_out) ) ! Convert to square modulus of charge density

    kmesh      = [ ( real(ix)*dk, ix = 0, sc-1 ) ] ! Set up k-components
    kmesh(nk:) = kmesh(nk:) - real(sc)*dk          ! in wraparound convention
    pot = 0.0

    ! Triple loop over xyz grid points (uses wraparound indexing)
    DO ix = 0, sc-1
       kx = kmesh(ix)
       fx = sharpen( 0.5*kx*dr )
       DO iy = 0, sc-1
          ky  = kmesh(iy)
          fxy = fx * sharpen( 0.5*ky*dr )
          DO iz = 0, sc-1
             kz   = kmesh(iz)
             fxyz = fxy * sharpen( 0.5*kz*dr )

             IF ( ix==0 .AND. iy==0 .AND. iz==0 ) CYCLE ! Skip zero wave vector

             k_sq = kx**2 + ky**2 + kz**2                          ! Squared magnitude of wave vector
             g    = (fourpi/k_sq) * EXP ( -k_sq / (4.0*kappa**2) ) ! Uncorrected influence function
             g    = g * fxyz                                       ! Apply simple correction
             pot  = pot + g * rho_sq(ix,iy,iz)

          END DO
       END DO
    END DO

    ! Divide by 2; box volume is 1
    pot = pot / 2.0

    ! Subtract self part of k-space sum
    pot = pot - kappa * SUM ( q(:)**2 ) / SQRT(pi)

    DEALLOCATE ( fft_inp, fft_out, rho_sq )

  END FUNCTION pot_k_pm_ewald

END MODULE ewald_module


text:
      SUBROUTINE XERABT(MESSG,NMESSG)                                   XER   10
C
C     ABSTRACT
C        ***NOTE*** MACHINE DEPENDENT ROUTINE
C        XERABT ABORTS THE EXECUTION OF THE PROGRAM.
C        THE ERROR MESSAGE CAUSING THE ABORT IS GIVEN IN THE CALLING
C        SEQUENCE IN CASE ONE NEEDS IT FOR PRINTING ON A DAYFILE,
C        FOR EXAMPLE.
C
C     DESCRIPTION OF PARAMETERS
C        MESSG AND NMESSG ARE AS IN XERROR, EXCEPT THAT NMESSG MAY
C        BE ZERO, IN WHICH CASE NO MESSAGE IS BEING SUPPLIED.
C
C     WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE
C END OF ABSTRACT
C     LATEST REVISION ---  7 JUNE 1978
C
      DIMENSION MESSG(NMESSG)
      IF (.TRUE.) STOP
      RETURN
      END

max_stars_repo_name: jhu-cisst/cisstNetlib

max_issues_repo_name: jhu-cisst/cisstNetlib

max_forks_repo_name: jhu-cisst/cisstNetlib


text:
program adapter_main
    use adapter_module, only: client_t, computer_t, mac_t, windows_t, windows_adapter_t
    implicit none
    type(client_t) :: client
    type(mac_t) :: mac
    type(windows_t), target :: windows
    type(windows_adapter_t) :: windows_adapter

    call client%insert_lightning_connector_into_computer(mac)
    windows_adapter%windows_machine => windows
    call client%insert_lightning_connector_into_computer(windows_adapter)

end program adapter_main

!> Results shall be:

!  Client inserts Lightning connector into computer.
!  Lightning connector is plugged into mac machine.
!  Client inserts Lightning connector into computer.
!  Adapter converts Lightning signal to USB.
!  USB connector is plugged into windows machine.

max_stars_repo_name: zoziha/Fortran-Design-Patterns

max_issues_repo_name: zoziha/Fortran-Design-Patterns

max_forks_repo_name: zoziha/Fortran-Design-Patterns


text:
  subroutine statSDCM (sDCM, u, v, partin, ipatne, mv, np, results)

! ======================================================================
!
!     Determining secondary particles characteristics for
!     gamma + N --> 2 pi interaction with statistical model.
!
!     See Barashenkov, et al. Nucl. Phys. A231 462 (1974).
!
!  Definition of partin:
!                       partin(4); sin(theta), direction of momentum
!                       partin(5); cos(theta), direction of momentum
!                       partin(6); sin(phi), direction of momentum
!                       partin(7); cos(phi), direction of momentum
!                       partin(8); kineti! energy of particle
!                       partin(9); rest mass of particle
!
!  Definition of ipatne:
!                       ipatne(1); charge of target nucleon
!
!    CEM95 written by S. G. Mashnik
!
!    Edited by A. J. Sierk,  LANL  T-2  February, 1996.
!    Edited by AJS,  August, 1997.
!    Modification: 22-oct-1998 by NVMokhov
!   "Last" change: 13-AUG-2003 by NVMokhov
!    Edited by A. J. Sierk, LANL T-16, October, 2003.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection)
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use standardDCMParams, only: zro, one, two, thrd, twthrd, pi, twpi, &
         & emneut, emprot, massPiPM, massPi0
    use standardDCMData,   only: photonEG

    implicit none
    class(StandardDCM), intent(inout) :: sDCm
    real(real64),   intent(in   ) :: u
    real(real64),   intent(in   ) :: v(3)
    real(real64),   intent(in   ) :: partin(9)
    integer(int32), intent(in   ) :: ipatne(5)
    integer(int32), intent(in   ) :: mv
    integer(int32), intent(  out) :: np
    class(StandardDCMResults), intent(inout) :: results

    real(real64) :: aux, b1, ct1, ct2, ct3, e1, e2, emnu, emnupi, empi, &
         & empi2, epim, f1, fi1, fi2, fi3, p1, p2, p3, st1, st2, t1, &
         & t2, t3, temp, temp1, temp2, temp3, temp4, tpim
    real(real64), dimension(3) :: pv1=zro, pv2=zro, pv3=zro, ps1=zro, &
         & ps2=zro, ps3=zro, pin=zro, pinst=zro

! ======================================================================

    results%imemo(2,mv+1) = 0
    results%imemo(3,mv+1) = 0
    results%imemo(4,mv+1) = 0
    results%imemo(2,mv+2) = 0
    results%imemo(3,mv+2) = 0
    results%imemo(4,mv+2) = 0
    results%imemo(2,mv+3) = 0
    results%imemo(3,mv+3) = 0
    results%imemo(4,mv+3) = 1
!  Determine randomly the resulting particle types
    temp4 = sDCM%rang()
    if (ipatne(1) <= 0) then
!  I.  gamma + n case:
       results%imemo(1,mv+3) = 0
       results%pmemo(9,mv+3) = emneut
       if (temp4 <= thrd) then
!   n + 2 pi0
          results%imemo(1,mv+1) = 0
          results%pmemo(9,mv+1) = massPi0
          results%imemo(1,mv+2) = 0
          results%pmemo(9,mv+2) = massPi0
       elseif (temp4 < twthrd) then
!   p + pi0 + pi-
          results%imemo(1,mv+1) = 0
          results%pmemo(9,mv+1) = massPi0
          results%imemo(1,mv+2) = -1
          results%pmemo(9,mv+2) = massPiPM
          results%imemo(1,mv+3) = 1
          results%pmemo(9,mv+3) = emprot
       else
!   n + pi- + pi+
          results%imemo(1,mv+1) = -1
          results%pmemo(9,mv+1) = massPiPM
          results%imemo(1,mv+2) = 1
          results%pmemo(9,mv+2) = massPiPM
       endif
    else
!  II. gamma + p case:
       results%imemo(1,mv+3) = 1
       results%pmemo(9,mv+3) = emprot
       if (temp4 <= thrd) then
!   p + pi0 + pi0
          results%imemo(1,mv+1) = 0
          results%pmemo(9,mv+1) = massPi0
          results%imemo(1,mv+2) = 0
          results%pmemo(9,mv+2) = massPi0
       elseif (temp4 <= twthrd) then
!    n + pi+ + pi0
          results%imemo(1,mv+1) = 0
          results%pmemo(9,mv+1) = massPi0
          results%imemo(1,mv+2) = 1
          results%pmemo(9,mv+2) = massPiPM
          results%imemo(1,mv+3) = 0
          results%pmemo(9,mv+3) = emneut
       else
!    p + pi+ + pi-
          results%imemo(1,mv+1) = 1
          results%pmemo(9,mv+1) = massPiPM
          results%imemo(1,mv+2) = -1
          results%pmemo(9,mv+2) = massPiPM
       endif
    endif
    np = 3

!  Other combinations might be possible!?
!  HOW, if we restrict to N + 2 pi?  AJS  10/2/03
    empi = max(results%pmemo(9,mv+1), results%pmemo(9,mv+2))
    empi2 = empi**2
    emnu = results%pmemo(9,mv+3)
    emnupi = emnu + min(results%pmemo(9,mv+1), results%pmemo(9,mv+2))
    temp = u
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(sDCM%io%message,1000) "130 and 137"
       call sDCM%io%print(4, 3, sDCM%io%message)
    end if
    epim = (temp**2 + empi2 - emnupi**2)/(two*temp)
    tpim = epim - empi
10  t1 = sDCM%rang()*tpim
    t2 = sDCM%rang()*tpim
    e1 = t1 + results%pmemo(9,mv+1)
    e2 = t2 + results%pmemo(9,mv+2)
    aux = u - e1 - e2
    f1 = 27.d0*e1*e2*aux/(temp**3)
    b1 = sDCM%rang()
    if (b1 < f1)then
       t3 = aux - results%pmemo(9,mv+3)
    else
       go to 10
    endif
    if (t3 <= zro) go to 10
    p1 = sqrt(abs(t1*(t1 + two*results%pmemo(9,mv+1))))
    p2 = sqrt(abs(t2*(t2 + two*results%pmemo(9,mv+2))))
    p3 = sqrt(abs(t3*(t3 + two*results%pmemo(9,mv+3))))
    temp1 = (p1 + p2 - p3)*(p1 - p2 + p3)*(p2 + p3 - p1)
    if (temp1 <= zro) go to 10
    ct3 = one - two*sDCM%rang()
    fi3 = twpi*sDCM%rang()
    temp2 = sqrt(abs(one - ct3**2))
    pv3(1) = p3*temp2*cos(fi3)
    pv3(2) = p3*temp2*sin(fi3)
    pv3(3) = p3*ct3
    temp3 = sqrt(abs(partin(8)*(partin(8) + two*partin(9))))
    pin(1) = temp3*partin(4)*partin(7)
    pin(2) = temp3*partin(4)*partin(6)
    pin(3) = temp3*partin(5)
    call photonEG%lorentzTransform (pin, v, pinst, partin(8), partin(9))
    call sDCM%rotor (pinst, v, pv3, ps3)
    temp = two*p3*p1
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(sDCM%io%message,1000) "167"
       call sDCM%io%print(4, 3, sDCM%io%message)
    end if
    ct1 = (p2**2 - p1**2 - p3**2)/(temp)
    temp = two*p3*p2
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(sDCM%io%message,1000) "173"
       call sDCM%io%print(4, 3, sDCM%io%message)
    end if
    ct2 = (p1**2 - p2**2 - p3**2)/(temp)
    fi1 = twpi*sDCM%rang()
    fi2 = pi + fi1
    st1 = sqrt(abs(one - ct1**2))
    st2 = sqrt(abs(one - ct2**2))
    pv1(1) = p1*st1*cos(fi1)
    pv1(2) = p1*st1*sin(fi1)
    pv1(3) = p1*ct1
    call sDCM%rotor (ps3, v, pv1, ps1)
    pv2(1) = p2*st2*cos(fi2)
    pv2(2) = p2*st2*sin(fi2)
    pv2(3) = p2*ct2
    call sDCM%rotor (ps3, v, pv2, ps2)
    results%pmemo(1,mv+1) = ps1(1)
    results%pmemo(2,mv+1) = ps1(2)
    results%pmemo(3,mv+1) = ps1(3)
    results%pmemo(1,mv+2) = ps2(1)
    results%pmemo(2,mv+2) = ps2(2)
    results%pmemo(3,mv+2) = ps2(3)
    results%pmemo(1,mv+3) = ps3(1)
    results%pmemo(2,mv+3) = ps3(2)
    results%pmemo(3,mv+3) = ps3(3)
    return

! ======================================================================
1000 format("Divide by zero error prevented in 'stat.f90' line(s) ", A)
! ======================================================================
  end subroutine statSDCM

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
module Source_IO
   use Environment

   implicit none

   ! Структура данных для хранения строки исходного текста
   type TextLine
      character(:, CH_), allocatable :: String
      type(TextLine), pointer        :: Next  => Null()
   end type TextLine

contains

   ! Получение значений из файла конфигурации
   subroutine Read_Config_File(ConfigFile, First, Last, K)
      character(*), intent(in)  :: ConfigFile
      integer,      intent(out) :: First, Last, K

      integer                   :: In

      open (file=ConfigFile, newunit=In)
         read(In, *) First, Last, K
      close (In)

      if (.not. (First <= Last .and. (K < First .or. K > Last))) then
         write(*,*) 'Incorrect settings'
         STOP
      endif
   end subroutine Read_Config_File

   ! Чтение исходного кода
   function Read_Source(InputFile) result (Source)
      character(*), intent(in) :: InputFile

      type(TextLine), pointer :: Source
      integer                 :: In

      open (file=InputFile, encoding=E_, newunit=In)
         Source => Read_Source_Line(In)
      close (In)
   end function Read_Source

   ! Чтение строки исходного кода
   recursive function Read_Source_Line(in) result(Code)
      type(TextLine), pointer :: Code
      integer, intent(in)       :: In
      integer, parameter        :: max_len = 1024
      character(max_len, CH_)   :: string
      integer                   :: IO

      ! Чтение строки во временную строку бОльшей длины
      read (In, "(a)", iostat=IO) string
      call Handle_IO_Status(IO, "reading line from source code")
      if (IO == 0) then
         allocate (Code)
         ! Хранение в размещаемом поле символов без завершающих пробелов
         Code%String = Trim(string)
         Code%Next => Read_Source_Line(In)
      else
         Code => Null()
      end if
   end function Read_Source_Line

   ! Вывод исходного кода
   subroutine Output_Source_Code(OutputFile, Code)
      character(*), intent(in)      :: OutputFile
      type(TextLine), intent(in)  :: Code
      integer  :: Out

      open (file=OutputFile, encoding=E_, newunit=Out)
         call Output_Source_Line(Out, Code)
      close (Out)
   end subroutine Output_Source_Code

   ! Вывод строки исходного кода
   recursive subroutine Output_Source_Line(Out, Code)
      integer, intent(in)           :: Out
      type(TextLine), intent(in)  :: Code
      integer  :: IO

      write (Out, "(a)", iostat=IO) Code%String
      call Handle_IO_Status(IO, "writing line to file")
      if (Associated(Code%next)) &
         call Output_Source_Line(Out, Code%next)
   end subroutine Output_Source_Line
end module Source_IO

max_stars_repo_name: taxnuke/algorithm-exercises

max_issues_repo_name: taxnuke/algorithm-exercises

max_forks_repo_name: taxnuke/algorithm-exercises


text:
      SUBROUTINE ZEROIN(F,B,C,RE,AE,IFLAG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     SANDIA MATHEMATICAL PROGRAM LIBRARY
C     APPLIED MATHEMATICS DIVISION 2646
C     SANDIA LABORATORIES
C     ALBUQUERQUE, NEW MEXICO  87185
C     CONTROL DATA 6600/7600  VERSION 8.1  AUGUST 1980
C                   *************************
C                   *       ISSUED BY       *
C                   *  SANDIA LABORATORIES, *
C                   *   A PRIME CONTRACTOR  *
C                   ********     TO THE     *
C                          *  UNITED STATES *
C                          *   DEPARTMENT   *
C                          *       OF       *
C                          *     ENERGY     *
C      *********************  ---NOTICE---  *********************
C      *THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED*
C      *  BY THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED  *
C      *   STATES NOR THE UNITED STATES DEPARTMENT OF ENERGY,   *
C      *               NOR ANY OF THEIR EMPLOYEES,              *
C      * NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR *
C      * EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR  *
C      * ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE  *
C      *          **********    ACCURACY,   **********          *
C      *          *        *  COMPLETENESS  *        *          *
C      *          *        *  OR USEFULNESS *        *          *
C      *          *        *     OF ANY     *        *          *
C      *          *        *  INFORMATION,  *        *          *
C      *          *        *   APPARATUS,   *        *          *
C      *       ****        *     PRODUCT    *        ****       *
C      *       *           *   OR PROCESS   *           *       *
C      *       *           *   DISCLOSED,   *           *       *
C      *       *           *  OR REPRESENTS *           *       *
C      *       *          **    THAT ITS    **          *       *
C      *       *          **  USE WOULD NOT **          *       *
C      *********          **    INFRINGE    **          *********
C                         **    PRIVATELY   **
C                         **      OWNED     **
C                         **     RIGHTS.    **
C                         **                **
C                         **                **
C                         **                **
C                         ********************
C
C     BASED ON A METHOD BY T J DEKKER
C     WRITTEN BY L F SHAMPINE AND H A WATTS
C     MODIFIED FOR THE MATH LIBRARY BY C B BAILEY
C
C     ABSTRACT
C        ZEROIN SEARCHES FOR A ZERO OF A FUNCTION F(X) BETWEEN
C        THE GIVEN VALUES B AND C UNTIL THE WIDTH OF THE INTERVAL
C        (B,C) HAS COLLAPSED TO WITHIN A TOLERANCE SPECIFIED BY
C        THE STOPPING CRITERION, ABS(B-C) .LE. 2.*(RW*ABS(B)+AE).
C        THE METHOD USED IS AN EFFICIENT COMBINATION OF BISECTION AND
C        THE SECANT RULE.  IN ORDER TO INSURE THAT ZEROIN WILL CONVERGE
C        TO A ZERO, THE USER SHOULD PICK VALUES FOR B AND C AT WHICH
C        THE FUNCTION DIFFERS IN SIGN.
C
C     DESCRIPTION OF ARGUMENTS
C     F,B,C,RE AND AE ARE INPUT PARAMETERS
C     B,C AND IFLAG ARE OUTPUT PARAMETERS
C        F     - NAME OF THE REAL VALUED EXTERNAL FUNCTION.  THIS NAME
C                MUST BE IN AN EXTERNAL STATEMENT IN THE CALLING
C                PROGRAM.  F MUST BE A FUNCTION OF ONE REAL ARGUMENT.
C        B     - ONE END OF THE INTERVAL (B,C).  THE VALUE RETURNED FOR
C                B USUALLY IS THE BETTER APPROXIMATION TO A ZERO OF F.
C        C     - THE OTHER END OF THE INTERVAL (B,C)
C        RE    - RELATIVE ERROR USED FOR RW IN THE STOPPING CRITERION.
C                IF THE REQUESTED RE IS LESS THAN MACHINE PRECISION,
C                THEN RW IS SET TO APPROXIMATELY MACHINE PRECISION.
C        AE    - ABSOLUTE ERROR USED IN THE STOPPING CRITERION.  IF THE
C                GIVEN INTERVAL (B,C) CONTAINS THE ORIGIN, THEN A
C                NONZERO VALUE SHOULD BE CHOSEN FOR AE.
C        IFLAG - A STATUS CODE.  USER MUST CHECK IFLAG AFTER EACH CALL.
C                CONTROL RETURNS TO THE USER FROM ZEROIN IN ALL CASES.
C                XERROR DOES NOT PROCESS DIAGNOSTICS IN THESE CASES.
C                 1 B IS WITHIN THE REQUESTED TOLERANCE OF A ZERO.
C                   THE INTERVAL (B,C) COLLAPSED TO THE REQUESTED
C                   TOLERANCE, THE FUNCTION CHANGES SIGN IN (B,C), AND
C                   F(X) DECREASED IN MAGNITUDE AS (B,C) COLLAPSED.
C                 2 F(B) = 0.  HOWEVER, THE INTERVAL (B,C) MAY NOT HAVE
C                   COLLAPSED TO THE REQUESTED TOLERANCE.
C                 3 B MAY BE NEAR A SINGULAR POINT OF F(X).
C                   THE INTERVAL (B,C) COLLAPSED TO THE REQUESTED
C                   TOLERANCE AND THE FUNCTION CHANGES SIGN IN (B,C) BUT
C                   F(X) INCREASED IN MAGNITUDE AS (B,C) COLLAPSED,I.E.
C                     ABS(F(B OUT)) .GT. MAX(ABS(F(B IN)),ABS(F(C IN)))
C                 4 NO CHANGE IN SIGN OF F(X) WAS FOUND ALTHOUGH THE
C                   INTERVAL (B,C) COLLAPSED TO THE REQUESTED TOLERANCE.
C                   THE USER MUST EXAMINE THIS CASE AND DECIDE WHETHER
C                   B IS NEAR A LOCAL MINIMUM OF F(X), OR B IS NEAR A
C                   ZERO OF EVEN MULTIPLICITY, OR NEITHER OF THESE.
C                 5 TOO MANY (.GT. 500) FUNCTION EVALUATIONS USED.
C
C     REFERENCES
C       1.  L F SHAMPINE AND H A WATTS, ZEROIN, A ROOT-SOLVING CODE,
C           SC-TM-70-631, SEPT 1970.
C       2.  T J DEKKER, FINDING A ZERO BY MEANS OF SUCCESSIVE LINEAR
C           INTERPOLATION, *CONSTRUCTIVE ASPECTS OF THE FUNDAMENTAL
C           THEOREM OF ALGEBRA*, EDITED BY B DEJON AND P HENRICI, 1969.
C
C
C     ER IS TWO TIMES THE COMPUTER UNIT ROUNDOFF VALUE WHICH IS
C     DEFINED HERE TO BE THE VALUE FOR THE IBM PC DOUBLE PRECISION
C
      ER = 2.0D-13
C
C     INITIALIZE
      RW=DMAX1(RE,ER)
      AW=DMAX1(AE,0.0D0)
      IC=0
      ACBS=DABS(B-C)
      A=C
      T=A
      FA=F(T)
      T=B
      FB=F(T)
      FC=FA
      KOUNT=2
      FX=DMAX1(DABS(FB),DABS(FC))
C
    1 IF (DABS(FC) .GE. DABS(FB)) GO TO 2
C     PERFORM INTERCHANGE
      A=B
      FA=FB
      B=C
      FB=FC
      C=A
      FC=FA
C
    2 IF (FB .EQ. 0.0D0) GO TO 11
      CMB=0.5D0*(C-B)
      ACMB=DABS(CMB)
      TOL=RW*DABS(B)+AW
C
C     TEST STOPPING CRITERION
      IF (ACMB .LE. TOL) GO TO 10
C
C     CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q
C     WHERE WE ARRANGE P .GE. 0.
C     THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW.
      P=(B-A)*FB
      Q=FA-FB
      IF (P .GE. 0.0D0) GO TO 3
      P=-P
      Q=-Q
C
C     UPDATE A AND CHECK FOR SATISFACTORY REDUCTION
C     IN THE SIZE OF OUR BOUNDING INTERVAL.
    3 A=B
      FA=FB
      IC=IC+1
      IF (IC .LT. 4) GO TO 4
      IF (8.0D0*ACMB .GE. ACBS) GO TO 6
      IC=0
      ACBS=ACMB
C
C     TEST FOR TOO SMALL A CHANGE
    4 IF (P .GT. DABS(Q)*TOL) GO TO 5
C
C     INCREMENT BY TOLERANCE
      B=B+DSIGN(TOL,CMB)
      GO TO 7
C
C     ROOT OUGHT TO BE BETWEEN B AND (C+B)/2.0D0
    5 IF (P .GE. CMB*Q) GO TO 6
C
C     INTERPOLATE
      B=B+P/Q
      GO TO 7
C
    6 B=0.5D0*(C+B)
C     BISECT
C
C     HAVE COMPLETED COMPUTATION FOR NEW ITERATE B
    7 T=B
      FB=F(T)
      IF (FB .EQ. 0.0D0) GO TO 11
C
C     DECIDE WHETHER NEXT STEP IS INTERPOLATION OR EXTRAPOLATION
      IF (DSIGN(1.0D0,FB) .NE. DSIGN(1.0D0,FC)) GO TO 8
      C=A
      FC=FA
    8 KOUNT=KOUNT+1
      IF (KOUNT .GT. 500) GO TO 15
      GO TO 1
C
C
C     FINISHED. PROCESS RESULTS FOR PROPER SETTING OF IFLAG
C
   10 IF (DSIGN(1.0D0,FB) .EQ. DSIGN(1.0D0,FC)) GO TO 13
      IF (DABS(FB) .GT. FX) GO TO 12
      IFLAG = 1
      RETURN
   11 IFLAG = 2
      RETURN
   12 IFLAG = 3
      RETURN
   13 IFLAG = 4
      RETURN
   15 IFLAG = 5
      RETURN
      END

max_stars_repo_name: marient/PelePhysics

max_issues_repo_name: marient/PelePhysics

max_forks_repo_name: marient/PelePhysics


text:
*=======================================================================
* File:    rdistm.f
* Author:  Alex Stivala
* Created: August 2008
*
* Read numeric SSE distance matrix.
*
* Each row is a row of the distance matrix, lower triangle
* only (since it is symmetric).
* The diagonal entries are meaningless (self-distance)
* and are included instead to specify the SSE type, with
* the following codes:
*
* 0.000 beta strand
* 1.000 alpha helix
* 2.000 pi helix
* 3.000 3_10 helix
*
* Each entry in matrix is in Angstroms format
* F6.3 with a space between each on a line, and one line
* per row of matrix.
*
* E.g.:
*
*   0.000
*   4.501  0.000
*  11.662 10.386  1.000
*  16.932 17.644  9.779  3.000
*  10.588 13.738 11.815 10.527  0.000
*  15.025 18.692 17.143 15.341  6.466  0.000
*  15.298 17.276 16.276 20.075 13.264 11.610  3.000
*   7.549 11.072 12.248 12.446  4.583  9.903 15.689  0.000
*
*
* This subroutine is used to read the distance matrix from the
* database, where the identifier and order and tableau have already
* been read by subroutine RDTABD.
*
* $Id: rdistm.F 2969 2009-11-22 01:29:39Z astivala $
*=======================================================================

      subroutine rdistm(nunit, dmat, ldo, n)

      implicit none
*
* rdistm - read a numeric SSE distance matrix
*
*
*     .. Scalar Arguments ..
      integer nunit,ldo,n
*     ..
*     .. Array Arguments
      double precision dmat(ldo, *)
*     ..
*
*
* Arguments
* =========
*
* nunit  (input) INTEGER
*        Unit number to read numeric distance matrix from
*
* dmat  (output) DOUBLE PRECISION array, dimension (ldo, ldo)
*        Distance matrix read from the file. Must be dimensioned large
*        enough to hold dimension (n,n) data.
*        The full symmetric matrix is stored.
*
* ldo   (input) INTEGER
*       Leading dimension of distance matrix array.
*
* n     (input) INTEGER
*       Order of distance matrix (dimension of matrix).
*       n <= ldo
*
*=======================================================================
*

*     .. Parameters ..
*     ..
*     .. Local Scalars ..
      integer i,j
      real dist
      character line*2048
*     ..
*     .. Local Arrays ..
*     ..
*     .. Intrinsic Functions ..
      intrinsic dble
*     ..
*     .. External Subroutines and Functions ..
*     ..
*     .. Executable Statements ..
*
      do 50 i = 1, n
         read(nunit,'(a2048)') line
         do 40 j = 1, i
            read(line((j-1)*7+1:(j-1)*7+7), '(f6.3)') dist
            dmat(i,j) = dble(dist)
            dmat(j,i) = dmat(i,j)
 40      continue
 50   continue
      return
      end

max_stars_repo_name: stivalaa/qptabsearch

max_issues_repo_name: stivalaa/qptabsearch

max_forks_repo_name: stivalaa/qptabsearch


text:
      program t
* test for e, en, es, and p in format-item-list
      real::a=3.141592654
 5678 format (3p,E15.7,-3p,E15.7e4,0p,EN15.5,EN10.3e5,ES10.3,ES10.3e1)
      write (*,5678) a,a,a,a,a,a
      endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
MODULE FUN_RandomColor
!!#### PURPOSE
!! Get a random color in RGBA real format, R,G,B,A in [0,1].

!!#### EXTERNAL PARAMETERS
USE KND_IntrinsicTypes,ONLY: KIND_Rsp,KIND_Rdp !!((01-A-KND_IntrinsicTypes.f90))

!!#### EXTERNAL PROCEDURES
USE SUB_Randomize                              !!((03-A-SUB_Randomize.f90))

!!#### DEFAULT IMPLICIT
IMPLICIT NONE

!!#### DEFAULT ACCESS
PRIVATE

!!#### PROCEDURE OVERLOADING
INTERFACE RandomColor
 MODULE PROCEDURE RandomColor_Rsp
 MODULE PROCEDURE RandomColor_Rdp
END INTERFACE

!!#### PUBLIC ACCESS LIST
PUBLIC :: RandomColor

CONTAINS

FUNCTION RandomColor_Rsp( A ) RESULT(RGBA)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp  !!((01-A-KND_IntrinsicTypes.f90))
REAL(KIND_R),INTENT(IN) :: A
REAL(KIND_R) :: RGBA(1:4)
!!--begin--
CALL Randomize( RGBA(1:3) )
RGBA(1) = A
!!--end--
END FUNCTION

FUNCTION RandomColor_Rdp( A ) RESULT(RGBA)
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp  !!((01-A-KND_IntrinsicTypes.f90))
REAL(KIND_R),INTENT(IN) :: A
REAL(KIND_R) :: RGBA(1:4)
!!--begin--
CALL Randomize( RGBA(1:3) )
RGBA(1) = A
!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
 SUBROUTINE contol(actio)
 !******************************************************************
 !*** READ control DATA
 !******************************************************************
 USE param_db,ONLY: mnam, mttl
 USE ctrl_db,ONLY: ndime, ndofn, nrotd, neulr, ntype, nacce, ifunc, ifixd, dtrec, endtm,  &
                   begtm, ncdlt, mscal, dtscal, dtuser, text, dtime,ttime, vefac,  &
                   lumped, itemp, therm, ndoft,                  &
                   tdtuse, tdtime, tdtsca, tmscal, tref, tscal, ascal
 USE c_input,ONLY: listen, exists, ludat, lures, get_name, getrea, getint,  &
                   card, backs, words
 USE outp_db, ONLY: energ, postype, toutp, toutd, ncdis, iener, cname, iwrit
 USE meshmo_db,ONLY: r_nrfset
 USE npo_db, ONLY : acceg
 USE esets_db, ONLY : nrenu
 USE gvar_db, ONLY : actchk,static
 USE static_db, ONLY : nides,miter,ittol,tol,sstep,pcont,nperi,lovec,nss,nds, &
                       sstop,autod,pmin,pmax,ipred,tolm,fdamp,spback
 IMPLICIT NONE

  !--- Dummy variables
  CHARACTER(len=*),INTENT(IN):: actio
  !--- Local variables
  INTEGER(kind=4):: i, j, point, dof
  LOGICAL done(5), dout,tout,cout

  INTERFACE
    SUBROUTINE readtp(toutp,t_d,task)
      IMPLICIT NONE
      CHARACTER(len=1), INTENT(IN) :: t_d
      INTEGER (kind=4), INTENT(IN) :: task
      REAL (kind=8), POINTER :: toutp(:)
    END SUBROUTINE readtp
  END INTERFACE

  !*** READ the first DATA card, and echo it immediately.
  IF (TRIM(actio) == 'NEW') THEN
    done = .FALSE.      !initializes all tasks to FALSE.
    !  assign default values for some parameters
    dtuser = 0d0        !time increment
    ncdlt = 0           !do not recompute
    dtscal = 0.75d0     !time increment scale factor
    mscal  = 1d0        !mass scale factor
    iwrit  = 1          !flag to echo data and results writing
    iener  = 0          !flat to compute energy parameters
    !        default values for heat transfer
    tscal  = 1.0d0      !time scale factor for thermal analysis
    therm  = .FALSE.    !thermo-mechanical analysis
    tdtuse = 0d0        !time increment for heat transfer
    tdtsca = 0.75d0     !time increment scale factor for heat transfer
    tmscal = 1d0        !mass scale factor for heat transfer
    ascal = .FALSE.     !don't auto-compute mass scale factor for thermo-mech coupling
    tref   = 293d0      !reference temperature (20 Celsius)
    itemp  = .FALSE.    !there are temperatures in the model
    ndoft  = 0          !no temperature DOFs
  ELSE
    done = .TRUE.       !initializes all flags to TRUE
    cout = .FALSE.      !only OUTPUT control is compulsory
  END IF

  CALL listen('CONTOL')            !read a line
  DO
    IF (exists('ENDCON')) EXIT     !end line found, exit loop
    ! C O N T R O L    D A T A   (new problem only)
    IF (exists('CONTRO')) THEN     !Set of control variables
      IF (TRIM(actio) /= 'NEW') CALL runend('CONTOL: CONTROL Params cannot chang')  !check
      WRITE (lures,"(/,'  C O N T R O L   P A R A M E T E R S ',/)",ERR=9999)
      ndime=getint('NDIME ',3,' Number of coordinate Dimension ...')
      IF( ndime == 2 )ntype=getint('NTYPE ',3,' Problem type for 2 Dimensional ...')
      neulr=getint('NEULR ',0,' Nodal system code 0:NO 1:YES .....')
      ndofn=getint('NDOFN ',ndime,' Number of Degrees Of Freedom .....')
      nrotd = MIN(ndofn,6)
      IF( exists('CONSIS') )THEN
        lumped = .FALSE.  !consistent mass matrix
        WRITE (lures, "(/,'      Consistent Mass Matrix will be used',/)",ERR=9999)
        nrenu=getint('NRENU ',0,' Node Renumbering flag ............')
      END IF
      WRITE (lures, "(/,'Type of analysis:')",ERR=9999)
      IF (exists('COUPLE')) THEN
        WRITE (lures, "(/,'      coupled thermo-mechanical',/)",ERR=9999)
        therm = .TRUE.
      ELSE
        therm = .FALSE.
        WRITE (lures, "(/,'      mechanical only',/)",ERR=9999)
      END IF
      ndoft=getint('NDOFT ',1,' Number of Temperature DOFs..........')
      done(1)=.TRUE.    ! set first task as done
    END IF
    !  **** S T E P   T I M E   D A T A
    IF (exists('TIME  ') .OR. exists('STATIC')) THEN     !time control parameters
      WRITE (lures,"(/,'  T I M E   P A R A M E T E R S ',/)",ERR=9999)
      IF (TRIM(actio) == 'NEW') THEN     !compulsory variable for NEW problems
        begtm = 0d0 ! initial time of the new strategy
        endtm =getrea('ENDTM ',0d0,'!Final Time .......................')
      ELSE                         !optional for regular strategies
        begtm = ttime ! initial time of the new strategy (=final of the last
        endtm =getrea('ENDTM ',endtm,' Final Time .......................')
      END IF
      dtscal=getrea('DTSCAL',dtscal,' SCAle fact. for DT in auto mode...')
      IF (dtscal == 0d0) dtscal=0.75d0
      IF (dtscal > 0.8d0) THEN
        WRITE(*,"(/,' Warning! Scaling factor DTSCAL might be too HIGH')")
        WRITE(55,"(/,' WARNING! Scaling factor DTSCAL might be too HIGH',/)",ERR=9999)
      END IF
      ncdlt = getint('NCDLT ',ncdlt,' No of steps between DT computat..')
      IF (exists('STATIC')) THEN     !static control parameters
        static = .TRUE.
        spback= exists('SPRING')      !elastic springback type
        miter =getint('MITER ',10000,' Maximum number of iterations .....')
        ittol =getint('ITTOL ',1    ,' Tolerance type ...................')
        tol   =getrea('TOLER ',0.01d0,' Convergence Tolerance ............')
        tolm  =getrea('TOLMAX',5d0*tol,' Maximum Convergence Tolerance ...')
        IF( spback ) THEN
          sstep = 1
          fdamp = 1d0
          sstop = .TRUE.
          autod = .FALSE.
          nides = miter/100
          nperi = miter/10
          ipred = 1 ; lovec = 1 ; nss = 1 ; nds = 1
        ELSE
          sstep =getint('SSTEP ',10   ,' Number of steps for the stage ....')
          fdamp =getrea('FACTOR',1.4d0,' Damping factor between steps......')
          sstop = exists('STOP')        !stop if no convergence
          autod = exists('AUTODA')      !use automatic damping
          nides =getint('NIDES ',100  ,' Number of steps to apply displ....')
          nperi =getint('NPERI ',miter/10,' Minimum Number of iteration   ....')
          ipred =getint('IPRED',1,     ' Prediction order .................')
          lovec =getint('LVSET ',0    ,' Load (+) or Velocity (-) control .')
          nss   =getint('NSS   ',2    ,' Number of subdivisions of veloc. .')
          nds   =getint('NDS   ',5    ,' Number of displacements steps ....')
          pmin  =getint('PMIN  ',   1 ,' Minimum Damping period ...........')
          pmax  =getint('PMAX  ',miter,' Maximum Damping period ...........')
        END IF
      ELSE
        dtuser=getrea('DTUSER',dtuser,' Time Increment = 0.-> Auto calc...')
        IF( dtuser == 0d0 .AND. ncdlt == 0 ) ncdlt = 100    !frequency of dtime re-evaluation
        mscal =getrea('MSCALE',mscal,' Mass SCALing factor  .............')
        tscal =getrea('TSCALE',1d0,' Time SCALing factor for therm.ana.')
        IF (exists('DTUSER') .AND. dtuser > 0d0 )THEN
          IF( dtime < dtuser .AND. TRIM(actio) /= 'NEW')THEN
             WRITE (lures,"(/,'  W A R N I N G   DTIME modified by user',/ &
                        &     '  DTUSER:',e12.4,' Old_DTIME:',e12.4)",ERR=9999)dtuser,dtime
             WRITE (*,"(/,'  W A R N I N G   DTIME modified by user',/  &
                        &     '  DTUSER:',e12.4,' Old_DTIME:',e12.4)",ERR=9999)dtuser,dtime
          END IF
          dtime = dtuser
        END IF
        vefac=getrea('VEFAC',vefac,' Velocity smoothing factor ........')
        static = .FALSE.
      END IF
      done(2) = .TRUE.  ! set second task as done
    END IF
    !  **** O U T P U T   C O N T R O L
    IF (exists('OUTPUT')) THEN   ! control parameters for output results
      WRITE (lures,"(/,'  O U T P U T   P A R A M E T E R S ',/)",ERR=9999)
      point = ncdis/10           ! node point
      dof   = MOD(ncdis,10)      ! DOF
      IF (dof == 0) dof=ndime   ! arbitrary set to Y (2D) or Z (3D)
      point = getint('POINT ',point,' Node for Output Control ..........')
      dof   = getint('DOF   ',dof  ,' Dof  for Output Control ..........')
      ncdis = point*10+dof       ! keep in one variable
      ! see if new output control exists
      j = 0   !initializes number of pos-types
      IF (exists('COUTD ') .OR. exists('COUTP ')) THEN   !curve controlled
        cout = .TRUE.                                    !set to true
        j = j + 1                                        !increase counter
        IF (exists('CURVE ',i )) THEN                    !name of the curve exists
          cname = get_name(posin=i,stype='CURV')         !get curve name
        ELSE
          CALL runend(' CURVE key word is compulsory')   !unnamed curve
        END IF
      ELSE
        cout = .FALSE.                                   !no curve data
      END IF
      IF (exists('TOUTD ') .OR. exists('TOUTP ')) THEN   !time controlled
        tout = .TRUE.                                    !set to true
        j = j + 1                                        !increase counter
      ELSE
        tout = .FALSE.                                   !no time  data
      END IF
      IF (exists('DOUTD ') .OR. exists ('DOUTP ')) THEN  !displacement controlled
        dout = .TRUE.                                    !set to true
        j = j + 1                                        !increase counter
        IF (ncdis < 10) CALL runend('CONTOL: check POINT and DOF')
      ELSE
        dout = .FALSE.                                   !no displacement data
      END IF

      IF (j > 1) CALL runend('CONTOL: check output control pars ')

      IF (TRIM(actio) == 'NEW') THEN  !for a new problem
        IF (cout) THEN             ! if a curve data will be used as control
          ncdis = -1               !use negative as a flag
          toutd = getrea('COUTD ',0d0,'!History Output function Interval ..')  !frequency measured as curve incr.
          postype = 'C'                   !set to pos-type CURVE
          CALL readtp(toutp,postype,1)    !read frequency for FULL output
        ELSE IF (dout) THEN      ! if a node displ. will be used as control
          toutd = getrea('DOUTD ',0d0,'!History Output Travel Interval ...')   ! frequency measured as displ. incr.
          postype = 'D'                   ! set to pos-type DISPLACEMENT
          CALL readtp(toutp,postype,1)
        ELSE                     ! if TIME will be used as control
          toutd =getrea('TOUTD ',0d0,'!History Output Interval  .........')  ! time frequency for selected output
          postype = 'T'                   ! set to pos-type TIME
          CALL readtp(toutp,postype,1)    ! read frequency for FULL output
          ncdis = 0
        END IF
      ELSE                       ! for an existing problem
        IF (cout) THEN           ! if a curve data will be used as control
          ncdis = -1
          toutd = getrea('COUTD ',toutd,' History Output Function Interval .')  !frequency measured as curve incr.
          postype = 'C'                   !set to pos-type CURVE
          CALL readtp(toutp,postype,2)    ! read frequency for FULL output
        ELSE IF (dout) THEN      ! if a node displ. will be used as control
          toutd = getrea('DOUTD ',toutd,' History Output Travel Interval ...')  ! frequency measured as displ. incr.
          postype = 'D'                   ! set to pos-type DISPLACEMENT
          CALL readtp(toutp,postype,2)    ! read frequency for FULL output
        ELSE                     ! if TIME will be used as control
          ncdis = 0
          toutd = getrea('TOUTD ',toutd,' History Output Interval  .........')  ! time frequency for selected output
          postype = 'T'                   ! set to pos-type TIME
          CALL readtp(toutp,postype,2)    ! read frequency for FULL output
        END IF
      END IF

      IF( toutd == 0d0 )toutd = 1d0
      iwrit = getint('IWRIT ',iwrit,' Print-out code -      0:NO 1:Yes .')   !write flag
      iener = getint('IENER ',iener,' Calculation of Energy 0:NO 1:Yes .')   !energy computation flag
      done(3) = .TRUE.  ! set third task as done
      IF( static )THEN
        point = getint('CPOINT',point,' Node for Convergence Output.......')
        dof   = getint('CDOF  ',dof  ,' Dof  for Convergence Output.......')
        pcont = point*10+dof       ! keep in one variable
      END IF
      ! curve data position is set in NCDIS, this value is assigned in EQNUMS
      IF( cout ) CALL rdcurv('COUTP', cname )
    END IF
    ! G R O U N D   A C C E L E R A T I O N   D A T A
    IF (exists('ACCELE')) THEN
      !STOP 'ACCELEROGRAMs not active presently'
      WRITE (lures,"(/,'  A C C E L E R O G R A M  PARAMETERS ',/)",ERR=9999)
      dtrec = getrea('DTREC ',0d0,'!Time increment of Accelerogram....')
      nacce = getint('NACCE ',1,'!Points in the Accelerogram .......')
      ifixd = getint('IFIXD ',0,' Direction of excitation ......... ')
      ifunc = getint('IFUNC ',1,' Activate  0:Yes  /=0:No ......... ')
      done(4) = .TRUE.  ! set fourth task as done
    END IF
    ! C O U P L E D   T H E R M O - M E C H   D A T A
    IF (exists('THERM ')) THEN     !heat transfer control parameters
      WRITE (lures,"(/,'  H E A T   T R A N S F E R   P A R A M E T E R S',/)",ERR=9999)
      tdtuse=getrea('TDTUSE',tdtuse,' Time Increment = 0.-> Auto calc...')
      tdtsca=getrea('TDTSCA',tdtsca,' SCAle fact. for DT in auto mode...')
      tmscal=getrea('TMSCAL',tmscal,' Capacity SCALing factor ..........')
      IF(exists('AUTOSC')) &
        ascal = .TRUE. !automatic evaluation of MASS MATRIX SCALING FACTOR
      tref  =getrea('TREF  ',tref,' Reference temperature ............')
      IF(tdtsca == 0.) tdtsca = 0.75d0
      IF(tdtsca > 0.8d0)THEN
        WRITE(* ,"(/,' Warning! Scaling factor TDTSCA might be too HIGH')")
        WRITE(55,"(/,' WARNING! Scaling factor TDTSCA might be too HIGH',/)",ERR=9999)
      END IF
      ! curve data position is set in TCURV, this value is assigned in TEQNUMS
      itemp = .TRUE.
      done(5) = .TRUE.  ! set five task as done
    END IF

    CALL listen('CONTOL')
    IF (ALL(done(1:5)) .AND. TRIM(actio) == 'NEW') THEN
      IF (.NOT.exists('ENDCON')) backs = .TRUE.
      EXIT    !exit if all tasks done and no END_CONTROL card exists
    END IF
  END DO
  !-------       Check all information exists
  IF (.NOT.done(1)) CALL runend('CONTOL:  Specify CONTROL Parameters')
  IF (.NOT.done(2)) CALL runend('CONTOL:  Specify TIME Parameters   ')
  IF (.NOT.done(3)) CALL runend('CONTOL:  Specify OUTPUT Parameters ')
  IF (.NOT.done(4)) THEN  !ground accelerations not read
    dtrec = 0d0         !
    nacce = 1
    ifixd = 0
    ifunc = 1
  END IF
  IF (therm .AND. .NOT.done(5) )THEN
    CALL runend('CONTOL:  Specify THERMAL Parameters')
  ELSE IF (.NOT.therm )THEN
    WRITE (lures, "(/,'Type of analysis: mechanical',/)",ERR=9999)
  END IF
  !-- deal with arrays ACCEG and ENERG
  ALLOCATE(acceg(nacce,ndime), energ(6*iener))
  acceg = 0d0  ! to avoid problems in the debugger version
  energ = 0d0  ! to avoid problems in the debugger version
  !-------
  CALL listen('CONTOL')     !refinament control data
  IF (TRIM(actio) == 'NEW') THEN      !for a new problem
    r_nrfset = 0   !initializes remesh modifications flag
  END IF
  ! M E S H   M O D I F I C A T I O N S   D A T A
  IF (exists('MESHMO')) THEN   ! mesh modifications data exist
    CALL rdmeshmo()
  ELSE                  ! no data
    backs = .TRUE.
  END IF

  IF (ndime == 3 .AND. neulr > 0) neulr=9    ! Local system in 3D defined by a matrix
  IF(actchk) iwrit = 0

RETURN
 9999 CALL runen2('')
END SUBROUTINE contol

SUBROUTINE readtp(toutp,t_d,task)
!-----------------------------------------------------------------------
! read control parameters for full output
!-----------------------------------------------------------------------
USE param_db,ONLY: midn, mlin
USE c_input,ONLY: exists, words, lures, param, getrea
IMPLICIT NONE

  !--- Dummy variables
  CHARACTER(len=1),INTENT(IN):: t_d    ! Time or Displacement control
  INTEGER(kind=4),INTENT(IN):: task   ! 1=NEW  2=OLD
  REAL(kind=8),POINTER:: toutp(:)     ! INTENT(IN OUT)

  !--- Local variables
  CHARACTER(len=midn):: name
  CHARACTER(len=mlin):: label
  INTEGER (kind=4) :: i,j,k
  REAL (kind=8) :: defau, first

  name = t_d // 'OUTP '           !TOUTP or DOUTP to read from data file
  IF (t_d == 'C') THEN            ! curve controlled output
    label = ' Complete Output Curve IntervaL ...'
  ELSE IF (t_d == 'T') THEN       ! time controlled output
    label = ' Complete Output IntervaL .........'
  ELSE                            ! displacement controlled output
    label = ' Complete Output Travel IntervaL ..'
  END IF

  IF (task == 1) THEN            ! for a NEW problem
    label(1:1) = '!'             ! parameters are COMPULSORY
    defau = 1d0                  ! no value as default
  ELSE                           ! for old problems
    defau = toutp(2)             ! use previous value as default
  END IF

  IF (exists(TRIM(name),i)) THEN  ! if ?OUTP exists
    j = i
    DO
      IF ((LEN_TRIM(words(j+1)) > 0) .OR. (param(j+1) == 0d0)) EXIT
      j = j + 1
    END DO
    first = getrea(name,defau,label)
    k = j - i + 2                      !array dimension
    IF (task == 2) DEALLOCATE(toutp)   !for an OLD problem release space
    ALLOCATE(toutp(k))                 !get space for data
    toutp(1) = REAL(k-1,8)             !number of prescribed values
    toutp(2:k) = param(i:j)            !prescribed values for output
    IF (k > 2) WRITE(lures,"(5e15.5)",ERR=9999) toutp(3:k)
  END IF
  IF( toutp(2) == 0d0 ) toutp(2) = 1d0
 RETURN
 9999 CALL runen2('')
 END SUBROUTINE readtp

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


module hs_forcing_mod

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

use     constants_mod, only: KAPPA, CP_air, GRAV, rdgas, rvgas, hlv, tfreeze,pi

use           fms_mod, only: error_mesg, FATAL, file_exist,       &
                             open_namelist_file, check_nml_error, &
                             mpp_pe, mpp_root_pe, close_file,     &
                             write_version_number, stdlog,        &
                             uppercase

use  time_manager_mod, only: time_type

use  diag_manager_mod, only: register_diag_field, send_data

use  field_manager_mod, only: MODEL_ATMOS, parse
use tracer_manager_mod, only: query_method, get_tracer_index

implicit none
private

!-----------------------------------------------------------------------
!---------- interfaces ------------

   public :: hs_forcing, hs_forcing_init

!-----------------------------------------------------------------------
!-------------------- namelist -----------------------------------------

   real :: t_zero=315., t_strat=200., delh=60.0, delv=10., eps=0., sigma_b=0.7, CHI=0.0
   real :: P00 = 1.e5

!added by cw 9/16/04
   real :: phi0 = 0.0
!end cw addition

   real :: ka = -40. !  negative values are damping time in days
   real :: ks =  -4., kf = -1.

   logical :: do_conserve_energy = .true.

   real :: trflux = 1.e-5   !  surface flux for optional tracer
   real :: trsink = -4.     !  damping time for tracer

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

!modified by cw 9/16/04
   namelist /hs_forcing_nml/  t_zero, t_strat, delh, delv, eps, CHI, &
                              sigma_b, ka, ks, kf, do_conserve_energy, &
                              trflux, trsink, phi0
! end cw modification
!-----------------------------------------------------------------------

   character(len=128) :: version='$Id: hs_forcing.f90,v 1.2 2002/07/16 22:33:19 fms Exp $'
   character(len=128) :: tag='$Name: havana $'

   real :: tka, tks, vkf
   real :: trdamp

   integer :: id_teq, id_tdt, id_udt, id_vdt,  &
              id_tdt_diss, id_diss_heat
   real    :: missing_value = -1.e10
   character(len=10) :: mod_name = 'hs_forcing'

   logical :: do_init = .true.

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

contains

!#######################################################################

 subroutine hs_forcing ( is, ie, js, je, dt, Time, lat, p_half, p_full, &
                         u, v, t, r, um, vm, tm, rm, udt, vdt, tdt, rdt,&
                         mask, kbot )

!-----------------------------------------------------------------------
   integer, intent(in)                        :: is, ie, js, je
      real, intent(in)                        :: dt
 type(time_type), intent(in)                  :: Time
      real, intent(in),    dimension(:,:)     :: lat
      real, intent(in),    dimension(:,:,:)   :: p_half, p_full
      real, intent(in),    dimension(:,:,:)   :: u, v, t, um, vm, tm
      real, intent(in),    dimension(:,:,:,:) :: r, rm
      real, intent(inout), dimension(:,:,:)   :: udt, vdt, tdt
      real, intent(inout), dimension(:,:,:,:) :: rdt

      real, intent(in),    dimension(:,:,:), optional :: mask
   integer, intent(in),    dimension(:,:)  , optional :: kbot
!-----------------------------------------------------------------------
   real, dimension(size(t,1),size(t,2))           :: ps, diss_heat
   real, dimension(size(t,1),size(t,2),size(t,3)) :: ttnd, utnd, vtnd, teq, pmass
   real, dimension(size(r,1),size(r,2),size(r,3)) :: rst, rtnd
   integer :: i, j, k, kb, n
   logical :: used
   real    :: flux, sink, value
   character(len=128) :: scheme, params

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

     if (do_init) call error_mesg ('hs_forcing','hs_forcing_init has not been called', FATAL)

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

!     surface pressure

     if (present(kbot)) then
         do j=1,size(p_half,2)
         do i=1,size(p_half,1)
            kb = kbot(i,j)
            ps(i,j) = p_half(i,j,kb+1)
         enddo
         enddo
     else
            ps(:,:) = p_half(:,:,size(p_half,3))
     endif

!-----------------------------------------------------------------------
!     rayleigh damping of wind components near the surface

      call rayleigh_damping ( ps, p_full, u, v, utnd, vtnd, mask=mask )

      if (do_conserve_energy) then
         ttnd = -((um+.5*utnd*dt)*utnd + (vm+.5*vtnd*dt)*vtnd)/CP_air
         tdt = tdt + ttnd
         if (id_tdt_diss > 0) used = send_data ( id_tdt_diss, ttnd, Time, is, js)
       ! vertical integral of ke dissipation
         if ( id_diss_heat > 0 ) then
          do k = 1, size(t,3)
            pmass(:,:,k) = p_half(:,:,k+1)-p_half(:,:,k)
          enddo
          diss_heat = CP_air/GRAV * sum( ttnd*pmass, 3)
          used = send_data ( id_diss_heat, diss_heat, Time, is, js)
         endif
      endif

      udt = udt + utnd
      vdt = vdt + vtnd

      if (id_udt > 0) used = send_data ( id_udt, utnd, Time, is, js)
      if (id_vdt > 0) used = send_data ( id_vdt, vtnd, Time, is, js)

!-----------------------------------------------------------------------
!     thermal forcing for held & suarez (1994) benchmark calculation

      call newtonian_damping ( lat, ps, p_full, t, ttnd, teq, mask )

      tdt = tdt + ttnd

      if (id_tdt > 0) used = send_data ( id_tdt, ttnd, Time, is, js)
      if (id_teq > 0) used = send_data ( id_teq, teq,  Time, is, js)

!-----------------------------------------------------------------------
!     -------- tracers -------
      do n = 1, size(rdt,4)
         rtnd = 0.0
         flux = trflux
         sink = trsink
         if (query_method('tracer_sms', MODEL_ATMOS, n, scheme, params)) then
            if (uppercase(trim(scheme)) == 'NONE') cycle
            if (parse(params,'flux',value) == 1) flux = value
            if (parse(params,'sink',value) == 1) sink = value
         endif
         rst = rm(:,:,:,n) + dt*rdt(:,:,:,n)
         call tracer_source_sink ( flux, sink, p_half, rst, rtnd, kbot )
         rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd
      enddo

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

 end subroutine hs_forcing

!#######################################################################

 subroutine hs_forcing_init ( axes, Time )

!-----------------------------------------------------------------------
!
!           routine for initializing the model with an
!              initial condition at rest (u & v = 0)
!
!-----------------------------------------------------------------------

           integer, intent(in) :: axes(4)
   type(time_type), intent(in) :: Time

!-----------------------------------------------------------------------
   integer  unit, io, ierr

!     ----- read namelist -----

      if (file_exist('input.nml')) then
         unit = open_namelist_file ( )
         ierr=1; do while (ierr /= 0)
            read  (unit, nml=hs_forcing_nml, iostat=io, end=10)
            ierr = check_nml_error (io, 'hs_forcing_nml')
         enddo
  10     call close_file (unit)
      endif

!     ----- write version info and namelist to log file -----

      call write_version_number (version,tag)
      if (mpp_pe() == mpp_root_pe()) write (stdlog(),nml=hs_forcing_nml)

!     ----- compute coefficients -----

      if (ka < 0.) ka = -86400.*ka
      if (ks < 0.) ks = -86400.*ks
      if (kf < 0.) kf = -86400.*kf

      tka = 0.; if (ka > 0.) tka = 1./ka
      tks = 0.; if (ks > 0.) tks = 1./ks
      vkf = 0.; if (kf > 0.) vkf = 1./kf

!     ----- for tracers -----

      if (trsink < 0.) trsink = -86400.*trsink
      trdamp = 0.; if (trsink > 0.) trdamp = 1./trsink

!     ----- register diagnostic fields -----

      id_teq = register_diag_field ( mod_name, 'teq', axes(1:3), Time, &
                      'equilibrium temperature', 'deg_K'   , &
                      missing_value=missing_value, range=(/100.,400./) )

      id_tdt = register_diag_field ( mod_name, 'tdt_ndamp', axes(1:3), Time, &
                      'newtonian damping', 'deg_K/sec' ,    &
                       missing_value=missing_value     )

      id_udt = register_diag_field ( mod_name, 'udt_rdamp', axes(1:3), Time, &
                      'rayleigh damping (zonal wind)', 'm/s2',       &
                       missing_value=missing_value     )

      id_vdt = register_diag_field ( mod_name, 'vdt_rdamp', axes(1:3), Time, &
                      'rayleigh damping (meridional wind)', 'm/s2',  &
                       missing_value=missing_value     )

      if (do_conserve_energy) then
         id_tdt_diss = register_diag_field ( mod_name, 'tdt_diss_rdamp', axes(1:3), &
                   Time, 'Dissipative heating from Rayleigh damping', 'deg_K/sec',&
                   missing_value=missing_value     )

         id_diss_heat = register_diag_field ( mod_name, 'diss_heat_rdamp', axes(1:2), &
                   Time, 'Integrated dissipative heating for Rayleigh damping', 'W/m2')
      endif

      do_init  = .false.

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

!added by cw 9/16/04
      phi0 = phi0*pi/180.0
!end cw addition

 end subroutine hs_forcing_init

!#######################################################################

 subroutine hs_forcing_end

!-----------------------------------------------------------------------
!
!       routine for terminating held-suarez benchmark module
!             (this routine currently does nothing)
!
!-----------------------------------------------------------------------

 end subroutine hs_forcing_end

!#######################################################################

 subroutine newtonian_damping ( lat, ps, p_full, t, tdt, teq, mask )

!-----------------------------------------------------------------------
!
!   routine to compute thermal forcing for held & suarez (1994)
!   benchmark calculation.
!
!-----------------------------------------------------------------------

real, intent(in),  dimension(:,:)   :: lat, ps
real, intent(in),  dimension(:,:,:) :: p_full, t
real, intent(out), dimension(:,:,:) :: tdt, teq
real, intent(in),  dimension(:,:,:), optional :: mask

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

          real, dimension(size(t,1),size(t,2)) :: &
     sin_lat, sin_lat_2, cos_lat_2, t_star, cos_lat_4, &
     tstr, sigma, the, tfactr, rps, p_norm

       real, dimension(size(t,1),size(t,2),size(t,3)) :: tdamp

       integer :: k
       real    :: tcoeff, pref

!-----------------------------------------------------------------------
!------------latitudinal constants--------------------------------------

      sin_lat  (:,:) = sin(lat(:,:))
      sin_lat_2(:,:) = sin_lat(:,:)*sin_lat(:,:)
      cos_lat_2(:,:) = 1.0-sin_lat_2(:,:)
      cos_lat_4(:,:) = cos_lat_2(:,:)*cos_lat_2(:,:)

!modified by cw 9/16/04
! original:
!      t_star(:,:) = t_zero - delh*(sin_lat - sin(phi0))**2 + sin(phi0)**2
! modified by cw 1/29/12
      t_star(:,:) = t_zero - delh*(sin_lat - sin(phi0))**2 + delh*sin(phi0)**2 - 20.0*CHI*sin_lat
! end cw modification
!      t_star(:,:) = t_zero - delh*sin_lat_2(:,:) - eps*sin_lat(:,:)

      tstr  (:,:) = t_strat - eps*sin_lat(:,:)

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

      tcoeff = (tks-tka)/(1.0-sigma_b)
      pref = P00
      rps  = 1./ps

      do k = 1, size(t,3)

!  ----- compute equilibrium temperature (teq) -----
         p_norm(:,:) = p_full(:,:,k)/pref
!modified by cw 9/16/04
         the   (:,:) = t_star(:,:) - delv*(cos(lat-phi0))**2*log(p_norm(:,:))
!end cw modification
!         the   (:,:) = t_star(:,:) - delv*cos_lat_2(:,:)*log(p_norm(:,:))
         teq(:,:,k) = the(:,:)*(p_norm(:,:))**KAPPA
         teq(:,:,k) = max( teq(:,:,k), tstr(:,:) )

!  ----- compute damping -----
         sigma(:,:) = p_full(:,:,k)*rps(:,:)
         where (sigma(:,:) <= 1.0 .and. sigma(:,:) > sigma_b)
           tfactr(:,:) = tcoeff*(sigma(:,:)-sigma_b)
!modified by cw 9/16/04
!           tdamp(:,:,k) = tka + cos(lat - phi0)**8*tfactr(:,:)
!end cw modification
           tdamp(:,:,k) = tka + cos_lat_4(:,:)*tfactr(:,:)
         elsewhere
           tdamp(:,:,k) = tka
         endwhere

      enddo

!*** note: if the following loop uses vector notation for all indices
!          then the code will not run ??????

      do k=1,size(t,3)
         tdt(:,:,k) = -tdamp(:,:,k)*(t(:,:,k)-teq(:,:,k))
      enddo

      if (present(mask)) then
         tdt = tdt * mask
         teq = teq * mask
      endif

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

 end subroutine newtonian_damping

!#######################################################################

 subroutine rayleigh_damping ( ps, p_full, u, v, udt, vdt, mask )

!-----------------------------------------------------------------------
!
!           rayleigh damping of wind components near surface
!
!-----------------------------------------------------------------------

real, intent(in),  dimension(:,:)   :: ps
real, intent(in),  dimension(:,:,:) :: p_full, u, v
real, intent(out), dimension(:,:,:) :: udt, vdt
real, intent(in),  dimension(:,:,:), optional :: mask

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

real, dimension(size(u,1),size(u,2)) :: sigma, vfactr, rps

integer :: i,j,k
real    :: vcoeff

!-----------------------------------------------------------------------
!----------------compute damping----------------------------------------

      vcoeff = -vkf/(1.0-sigma_b)
      rps = 1./ps

      do k = 1, size(u,3)

         sigma(:,:) = p_full(:,:,k)*rps(:,:)

         where (sigma(:,:) <= 1.0 .and. sigma(:,:) > sigma_b)
            vfactr(:,:) = vcoeff*(sigma(:,:)-sigma_b)
            udt(:,:,k)  = vfactr(:,:)*u(:,:,k)
            vdt(:,:,k)  = vfactr(:,:)*v(:,:,k)
         elsewhere
            udt(:,:,k) = 0.0
            vdt(:,:,k) = 0.0
         endwhere
      enddo

      if (present(mask)) then
          udt = udt * mask
          vdt = vdt * mask
      endif

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

 end subroutine rayleigh_damping

!#######################################################################

 subroutine tracer_source_sink ( flux, damp, p_half, r, rdt, kbot )

!-----------------------------------------------------------------------
      real, intent(in)  :: flux, damp, p_half(:,:,:), r(:,:,:)
      real, intent(out) :: rdt(:,:,:)
   integer, intent(in), optional :: kbot(:,:)
!-----------------------------------------------------------------------
      real, dimension(size(r,1),size(r,2),size(r,3)) :: source, sink
      real, dimension(size(r,1),size(r,2))           :: pmass

      integer :: i, j, kb
      real    :: rdamp
!-----------------------------------------------------------------------

      rdamp = damp
      if (rdamp < 0.) rdamp = -86400.*rdamp   ! convert days to seconds
      if (rdamp > 0.) rdamp = 1./rdamp

!------------ simple surface source and global sink --------------------

      source(:,:,:)=0.0

   if (present(kbot)) then
      do j=1,size(r,2)
      do i=1,size(r,1)
         kb = kbot(i,j)
         pmass (i,j)    = p_half(i,j,kb+1) - p_half(i,j,kb)
         source(i,j,kb) = flux/pmass(i,j)
      enddo
      enddo
   else
         kb = size(r,3)
         pmass (:,:)    = p_half(:,:,kb+1) - p_half(:,:,kb)
         source(:,:,kb) = flux/pmass(:,:)
   endif

     sink(:,:,:) = rdamp*r(:,:,:)
     rdt(:,:,:) = source(:,:,:)-sink(:,:,:)

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

 end subroutine tracer_source_sink

!#######################################################################

end module hs_forcing_mod


program main
    use parser
    implicit none
    character (len=100) :: filepath, memory
    integer(kind = 8) :: imemory
    if (command_argument_count() /= 2) then
        print *, "Error: Two arguments needed; ./output <filepath> <memory size>"
        call EXIT(1)
    end if
    call get_command_argument(1, filepath)
    call get_command_argument(2, memory)
    read(memory, *) imemory
    call parse(filepath, imemory)
    print *, "Hello World"
end program main