text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> This module provide functions to generate Femap neutral output
module m_heat_lib_neutral
contains

  !----------------------------------------------------------------------
  !     heat_put_neutral_601 ( hecMESH )                                l
  !     heat_put_neutral_402 ( hecMESH )                                l
  !     heat_put_neutral_403 ( hecMESH )                                l
  !     heat_put_neutral_404 ( hecMESH )                                l
  !     heat_put_neutral_409 ( hecMESH )                                l
  !     heat_put_neutral_450 ( hecMESH )                                l
  !     heat_put_neutral_451 ( hecMESH )                                l
  !----------------------------------------------------------------------
  subroutine heat_put_neutral_601 ( hecMESH )

    use m_fstr

    implicit real(kind=kreal) (A - H, O - Z)
    dimension rdum(0:9)
    character(len=80) :: line

    type(hecmwST_local_mesh) :: hecMESH

    CP = 0.0d0
    CD = 0.0d0
    RHO = 0.0d0

    !C
    !C==put MATERIAL : ( BLOCK NO. = 601 )
    !C
    write(INEU,'(a)') '   -1'
    write(INEU,*) '  601'

    icou = 0
    do im1 = 1, hecMESH%material%n_mat
      line = ' '
      line(7:23) = ',-601,55,0,0,1,0,'
      write(line(1:6),'(i6)') im1
      write(INEU,'(a23)') line(1:23)
      write(INEU,*) '<NULL>'
      write(INEU,*) '10,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '25,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,'
      write(INEU,*) '200,'

      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'

      do im2= 1, 3
        icou = icou  + 1
        is= hecMESH%material%mat_TABLE_index(icou -1) + 1
        val  = hecMESH%material%mat_VAL (is)

        if (im2.eq.3) CD  = val
        if (im2.eq.2) CP  = val
        if (im2.eq.1) RHO = val
      enddo

      rdum( 0) = 0.d0
      rdum( 1) = 0.d0
      rdum( 2) = CD
      rdum( 3) = CD
      rdum( 4) = CD
      rdum( 5) = CD
      rdum( 6) = CD
      rdum( 7) = CD
      rdum( 8) = CP
      rdum( 9) = RHO

      write(INEU,'(1p,10(E9.2,'',''))') ( rdum(i),i=0,9 )

      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,'
      write(INEU,*) '50,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '70,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'

    enddo
    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_601

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

  !C
  subroutine heat_put_neutral_402 ( hecMESH )

    use m_fstr

    implicit real(kind=kreal) (A - H, O - Z)
    dimension idum(0:9)

    type(hecmwST_local_mesh) :: hecMESH

    !C
    !C==put Section : ( BLOCK NO. = 402 )
    !C
    write(INEU,'(a)') '   -1'
    write(INEU,*) '  402'

    do im = 1, hecMESH%section%n_sect
      idum(0) = im
      idum(1) = 110
      idum(2) = hecMESH%section%sect_mat_ID_item(im)
      if    ( hecMESH%section%sect_type(im) == 1 ) then
        idum(3) = 25
      elseif( hecMESH%section%sect_type(im) == 2 ) then
        idum(3) = 17
      elseif( hecMESH%section%sect_type(im) == 3 ) then
        idum(3) = 5
      elseif( hecMESH%section%sect_type(im) == 4 ) then
        idum(3) = 9
      endif
      idum(4) = 1
      idum(5) = 0

      write(INEU,'(6(i6,'',''))') ( idum(i),i=0,5 )
      write(INEU,*) '<NULL>'
      write(INEU,*) '0,0,0,0,'
      write(INEU,*) '90,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,'
      write(INEU,*) '0,0,'
      write(INEU,*) '190,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0.,0.,0.,0.,0.,'
      write(INEU,*) '0,'
      write(INEU,*) '0,'

    enddo

    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_402

  !----------------------------------------------------------------------
  !C
  subroutine heat_put_neutral_403 ( hecMESH )

    use m_fstr

    implicit real(kind=kreal) (A - H, O - Z)

    type(hecmwST_local_mesh) :: hecMESH

    !C
    !C==put NODE : ( BLOCK NO. = 403 )
    !C
    write(INEU,'(a)') '   -1'
    write(INEU,*) '  403'

    do i = 1, hecMESH%n_node

      inod = hecMESH%global_node_ID(i)
      xx = hecMESH%node( 3*i-2 )
      yy = hecMESH%node( 3*i-1 )
      zz = hecMESH%node( 3*i   )

      write(INEU,'(i8,a,3(e15.7,'',''))')  &
        &          inod,',0,0,1,46,0,0,0,0,0,0,',xx,yy,zz
    enddo

    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_403

  !----------------------------------------------------------------------
  !C
  subroutine heat_put_neutral_404 ( hecMESH )

    use m_fstr

    implicit real(kind=kreal) (A - H, O - Z)
    dimension nn(20), nna(10), nnb(10)

    type(hecmwST_local_mesh) :: hecMESH

    !C
    !C==put ELEMENT : ( BLOCK NO. = 404 )
    !C

    write(INEU,'(a)') '   -1'
    write(INEU,*) '  404'

    do i = 1, hecMESH%n_elem

      ielm = hecMESH%global_elem_ID(i)
      icol = 124
      isid = hecMESH%section_ID(i)
      isop = hecMESH%section%sect_opt(isid)

      is = hecMESH%elem_node_index(i-1) + 1
      ie = hecMESH%elem_node_index(i)

      k = 0
      do j = is, ie
        k = k + 1
        jj    = hecMESH%elem_node_item(j)
        nn(k) = hecMESH%global_node_ID(jj)
      enddo

      nna = 0  ;  nnb = 0
      ietyp = hecMESH%elem_type(i)
      if( ietyp == 231 ) then
        istyp = 25
        if( isop == 1 ) istyp = 19
        if( isop == 2 ) istyp = 35
        itopo = 2
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
      elseif( ietyp == 731 ) then
        istyp = 17
        itopo = 2
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
      elseif( ietyp == 232 ) then
        istyp = 26
        if( isop == 1 ) istyp = 20
        if( isop == 2 ) istyp = 36
        itopo = 3
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(5) = nn(4)
        nna(6) = nn(5)
        nna(7) = nn(6)
      elseif( ietyp == 732 ) then
        istyp = 18
        itopo = 3
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(5) = nn(4)
        nna(6) = nn(5)
        nna(7) = nn(6)
      elseif( ietyp == 241 ) then
        istyp = 25
        if( isop == 1 ) istyp = 19
        if( isop == 2 ) istyp = 35
        itopo = 4
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(4) = nn(4)
      elseif( ietyp == 741 ) then
        istyp = 17
        itopo = 4
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(4) = nn(4)
      elseif( ietyp == 242 ) then
        istyp = 26
        if( isop == 1 ) istyp = 20
        if( isop == 2 ) istyp = 36
        itopo = 5
        do ij = 1, 8
          nna(ij) = nn(ij)
        enddo
      elseif( ietyp == 742 ) then
        istyp = 18
        itopo = 5
        do ij = 1, 8
          nna(ij) = nn(ij)
        enddo
      elseif( ietyp == 341 ) then
        istyp = 25
        itopo = 6
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(5) = nn(4)
      elseif( ietyp == 351 ) then
        istyp = 25
        itopo = 7
        nna(1) = nn(1)
        nna(2) = nn(2)
        nna(3) = nn(3)
        nna(5) = nn(4)
        nna(6) = nn(5)
        nna(7) = nn(6)
      elseif( ietyp == 361 ) then
        istyp = 25
        itopo = 8
        do ij = 1, 8
          nna(ij) = nn(ij)
        enddo
      elseif( ietyp == 342 ) then
        istyp = 26
        itopo = 10
        nna( 1) = nn( 1)
        nna( 2) = nn( 2)
        nna( 3) = nn( 3)
        nna( 5) = nn( 4)
        nna( 9) = nn( 5)
        nna(10) = nn( 6)
        nnb( 1) = nn( 7)
        nnb( 3) = nn( 8)
        nnb( 4) = nn( 9)
        nnb( 5) = nn(10)
      elseif( ietyp == 352 ) then
        istyp = 26
        itopo = 11
        nna( 1) = nn( 1)
        nna( 2) = nn( 2)
        nna( 3) = nn( 3)
        nna( 5) = nn( 4)
        nna( 6) = nn( 5)
        nna( 7) = nn( 6)
        nna( 9) = nn( 7)
        nna(10) = nn( 8)
        nnb( 1) = nn( 9)
        nnb( 3) = nn(13)
        nnb( 4) = nn(14)
        nnb( 5) = nn(15)
        nnb( 7) = nn(10)
        nnb( 8) = nn(11)
        nnb( 9) = nn(12)
      elseif( ietyp == 362 ) then
        istyp = 26
        itopo = 12
        do ij = 1, 10
          nna(ij) = nn(ij)
        enddo
        nnb( 1) = nn(11)
        nnb( 2) = nn(12)
        nnb( 3) = nn(17)
        nnb( 4) = nn(18)
        nnb( 5) = nn(19)
        nnb( 6) = nn(20)
        nnb( 7) = nn(13)
        nnb( 8) = nn(14)
        nnb( 9) = nn(15)
        nnb(10) = nn(16)
      endif

      write(INEU,'(5(i8,'',''),a)') &
        &          ielm,icol,isid,istyp,itopo,'1,0,0,0,0,0,0,0,'
      write(INEU,'(10(i8,'',''))') (nna(j),j=1,10)
      write(INEU,'(10(i8,'',''))') (nnb(j),j=1,10)
      write(INEU,*) '0,0,0,'
      write(INEU,*) '0,0,0,'
      write(INEU,*) '0,0,0,'
      write(INEU,*) '0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,'

    enddo

    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_404

  !----------------------------------------------------------------------
  !C
  subroutine heat_put_neutral_409( INEU )

    !C
    !C==put View : BLOCK NO. = 409 )
    !C
    write(INEU,'(a)') '   -1'
    write(INEU,*) ' 409'
    write(INEU,*) '1,'
    write(INEU,*) 'Default XY View'
    write(INEU,*) '2,0,1,'
    !C      write(INEU,*) '35.2644,-45.,0.,'
    write(INEU,*) '0.,0.,0.,'
    write(INEU,*) '2.5,1.25,1.5,'
    write(INEU,*) '1.,1.,0,0.,0.,0.,0.,0.,0.,'
    write(INEU,*) '1.03572,0.51035,0.,'
    write(INEU,*) '1.2574,'
    write(INEU,*) '0.,0.,1.,1.,'
    write(INEU,*) '2,0,1,1,0,'
    write(INEU,*) '-1,-1,0,1,0,1,1,60031,0,4000000,'
    write(INEU,*) '9,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '0,0,0,'
    write(INEU,*) '100.,100.,1,7,'
    write(INEU,*) '0,1,1.,'
    write(INEU,*) '0.,0.,0.,'
    write(INEU,*) '0.,0.,1.,'
    write(INEU,*) '0,1,0,0,'
    write(INEU,*) '0.,0.,0.,'
    write(INEU,*) '1.,0.,0.,'
    write(INEU,*) '0.,0.,0.,'
    write(INEU,*) '0.,1.,0.,'
    write(INEU,*) '2.,1.,70.,0.5,'
    write(INEU,*) '0.,0.,0.,0.,0.,0,100.,1000.,0.,0.,0.,100.,-100.,1,'
    write(INEU,*) '5.,90.,10.,10.,1.,'
    write(INEU,*) '4,176,0,0,0,0,0,0,0,0,0.,0.,0.,'
    write(INEU,*) '0,0,0,0,0,0,14,110,'
    write(INEU,*) '0,1,1,1,0,1,0,1,1,0,1,1,0,1,1,1,0,0,1,'
    write(INEU,*) '0,0,0.00000001,25.,100.,0.,0.,0.,20,'
    write(INEU,*) '0,1,1,0,0,1,20.,0,'
    write(INEU,*) '12,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,0.,'
    write(INEU,*) '0,5,0,0,0,0.,25.,'
    write(INEU,*) '4,16408,20,16504,100,16488,'
    write(INEU,*) '0.,0.,'
    write(INEU,*) '0.,0.,0.,0.,'
    write(INEU,*) '9,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '1.,'
    write(INEU,*) '2,'
    write(INEU,*) '<NULL>'
    write(INEU,*) '<NULL>'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '0.,0.,0.,0.,'
    write(INEU,*) '0.,0.,0.,0.,'
    write(INEU,*) '90,1,124,1,0,'
    write(INEU,*) '0,60,0,0,'
    write(INEU,*) '0,24,0,0,'
    write(INEU,*) '0,100,0,0,'
    write(INEU,*) '0,2,0,0,'
    write(INEU,*) '0,24580,0,0,'
    write(INEU,*) '0,124,0,0,'
    write(INEU,*) '0,46,0,0,'
    write(INEU,*) '0,120,0,0,'
    write(INEU,*) '0,124,0,1,'
    write(INEU,*) '0,124,0,0,'
    write(INEU,*) '0,12,0,1,'
    write(INEU,*) '0,62,0,0,'
    write(INEU,*) '0,62,0,0,'
    write(INEU,*) '0,10,0,0,'
    write(INEU,*) '0,52,0,0,'
    write(INEU,*) '0,4,0,0,'
    write(INEU,*) '0,120,0,0,'
    write(INEU,*) '0,12,0,0,'
    write(INEU,*) '0,2,0,0,'
    write(INEU,*) '0,120,0,0,'
    write(INEU,*) '0,8312,0,0,'
    write(INEU,*) '0,24600,0,0,'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '1,74,0,1,'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '3,124,0,1,'
    write(INEU,*) '0,24636,0,0,'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '0,4,0,0,'
    write(INEU,*) '0,100,0,0,'
    write(INEU,*) '0,124,0,1,'
    write(INEU,*) '0,60,0,1,'
    write(INEU,*) '0,56,0,1,'
    write(INEU,*) '0,24,0,0,'
    write(INEU,*) '0,8216,0,1,'
    write(INEU,*) '0,4,0,0,'
    write(INEU,*) '0,124,2,0,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,0,1,'
    write(INEU,*) '1,124,5,1,'
    write(INEU,*) '0,0,0,1,'
    write(INEU,*) '0,24,0,1,'
    write(INEU,*) '0,124,0,0,'
    write(INEU,*) '0,100,0,1,'
    write(INEU,*) '1,100,0,1,'
    write(INEU,*) '0,0,0,1,'
    write(INEU,*) '0,16,0,0,'
    write(INEU,*) '0,124,4,1,'
    write(INEU,*) '0,62,0,0,'
    write(INEU,*) '2,124,1,1,'
    write(INEU,*) '1,8254,0,0,'
    write(INEU,*) '0,124,1,1,'
    write(INEU,*) '1,0,5,1,'
    write(INEU,*) '0,124,0,1,'
    write(INEU,*) '0,100,0,1,'
    write(INEU,*) '0,100,0,1,'
    write(INEU,*) '1,46,0,1,'
    write(INEU,*) '1,120,0,1,'
    write(INEU,*) '1,4,0,1,'
    write(INEU,*) '1,52,0,1,'
    write(INEU,*) '1,24,0,1,'
    write(INEU,*) '1,93,0,1,'
    write(INEU,*) '1,12,0,1,'
    write(INEU,*) '1,10,0,1,'
    write(INEU,*) '1,104,0,1,'
    write(INEU,*) '0,100,0,0,'
    write(INEU,*) '0,24,0,0,'
    write(INEU,*) '0,60,0,0,'
    write(INEU,*) '0,104,0,0,'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,4,0,0,'
    write(INEU,*) '0,0,1,0,'
    write(INEU,*) '0,0,0,0,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,62,1,1,'
    write(INEU,*) '0,60,4,0,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,*) '0,0,1,1,'
    write(INEU,'(a)') '-1,'
    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_409

  !----------------------------------------------------------------------
  !C
  subroutine heat_put_neutral_450( INEU )

    !C
    !C==put RESULT : BLOCK NO. = 450 )
    !C
    write(INEU,'(a)') '   -1'
    write(INEU,*) '  450'
    write(INEU,*) '1,'
    write(INEU,*) 'hecmw_FSTR_heat_result'
    write(INEU,*) '0,0,'
    write(INEU,*) '0.,'
    write(INEU,*) '1,'
    write(INEU,*) '<NULL>'
    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_450

  !----------------------------------------------------------------------
  !C
  subroutine heat_put_neutral_451 ( hecMESH, hecHEAT )

    use m_fstr

    implicit real(kind=kreal) (A - H, O - Z)

    type(fstr_heat)          :: hecHEAT
    type(hecmwST_local_mesh) :: hecMESH

    !C
    !C==put DISPLACEMENT : ( BLOCK NO. = 451 )
    !C    & STRESS/STRAIN
    !C

    tmin = 1.0e9
    tmax =-1.0e9

    !C

    do i = 1, hecMESH%n_node

      tt = hecHEAT%TEMP(i)
      if( tt > tmax ) tmax = tt
      if( tt < tmin ) tmin = tt

    enddo

    absmax = dabs(tmax)
    absmin = dabs(tmin)
    if( absmin > absmax ) absmax = absmin

    write(INEU,'(a)') '   -1'
    write(INEU,*) '  451'

    write(INEU,*) '1,1,1,'
    write(INEU,*) 'Temperature'
    write(INEU,'(3(e15.7,'',''))') tmin, tmax, absmax
    write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
    write(INEU,*) '0,0,0,0,0,0,0,0,0,0,'
    write(INEU,*) '0,0,6,7,'
    write(INEU,*) '0,0,1,'

    do i = 1, hecMESH%n_node
      inod = hecMESH%global_node_ID(i)
      tt = hecHEAT%TEMP(i)
      write(INEU,'(i8,'','',E15.7,'','')') inod,tt
    enddo
    write(INEU,*) '-1,0.'
    !C
    !C==TERMINATION
    !C
    write(INEU,'(a)') '   -1'

    return
  end subroutine heat_put_neutral_451
end module m_heat_lib_neutral

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
!-----------------------------------------------------------------------
!     Subroutine:    torlines_write
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          06/9/2015
!     Description:   This subroutine writes out the fieldline data
!-----------------------------------------------------------------------
      SUBROUTINE torlines_write(write_type)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE torlines_fieldlines
      USE torlines_background
      USE torlines_realspace
      USE torlines_runtime
      USE torlines_write_par
      USE ez_hdf5
      USE mpi_params
      USE mpi_inc
!-----------------------------------------------------------------------
!     Input Variables
!          write_type  Type of write to preform
!-----------------------------------------------------------------------
      IMPLICIT NONE
      CHARACTER(*), INTENT(IN):: write_type
!-----------------------------------------------------------------------
!     Local Variables
!          ierr        Error flag
!-----------------------------------------------------------------------
      INTEGER :: ier
!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      ier = 0
      SELECT CASE (TRIM(write_type))
         CASE('GRID_INIT')
            IF (myworkid == master) THEN
               CALL open_hdf5('torlines_'//TRIM(id_string)//'.h5',fid,ier,LCREATE=.true.)
               IF (ier /= 0) CALL handle_err(HDF5_OPEN_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
               CALL write_scalar_hdf5(fid,'VERSION',ier,DBLVAR=TORLINES_VERSION,ATT='Version Number',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'VERSION',ier)
               CALL write_scalar_hdf5(fid,'lvmec',ier,BOOVAR=lvmec,ATT='VMEC input',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lvmec',ier)
               CALL write_scalar_hdf5(fid,'lpies',ier,BOOVAR=lpies,ATT='PIES input',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lpies',ier)
               CALL write_scalar_hdf5(fid,'lspec',ier,BOOVAR=lspec,ATT='SPEC input',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lspec',ier)
               CALL write_scalar_hdf5(fid,'lcoil',ier,BOOVAR=lcoil,ATT='Coil input',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lcoil',ier)
               CALL write_scalar_hdf5(fid,'lmgrid',ier,BOOVAR=lmgrid,ATT='MGRID input',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lmgrid',ier)
               CALL write_scalar_hdf5(fid,'lvac',ier,BOOVAR=lvac,ATT='Vacuum calc',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'lvac',ier)
               CALL write_scalar_hdf5(fid,'nrho',ier,INTVAR=nrho,ATT='Number of Rho Gridpoints',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nrho',ier)
               CALL write_scalar_hdf5(fid,'nu',ier,INTVAR=nu,ATT='Number of Poloidal Gridpoints',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nu',ier)
               CALL write_scalar_hdf5(fid,'nv',ier,INTVAR=nv,ATT='Number of Toroidal Gridpoints',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nv',ier)
               CALL write_scalar_hdf5(fid,'nfp',ier,INTVAR=nfp,ATT='Number of Field Periods',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nfp',ier)
               CALL write_var_hdf5(fid,'rreal',nrho,nu,nv,ier,DBLVAR=rreal,ATT='R Gridpoints',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'rreal',ier)
               CALL write_var_hdf5(fid,'zreal',nrho,nu,nv,ier,DBLVAR=zreal,ATT='Z Gridpoints',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'zreal',ier)
               CALL close_hdf5(fid,ier)
               IF (ier /= 0) CALL handle_err(HDF5_CLOSE_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
            END IF
         CASE('FIELD')
            IF (myworkid == master) THEN
               CALL open_hdf5('torlines_'//TRIM(id_string)//'.h5',fid,ier,LCREATE=.false.)
               IF (ier /= 0) CALL handle_err(HDF5_OPEN_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
               CALL write_var_hdf5(fid,'bsreal',nrho,nu,nv,ier,DBLVAR=bsreal,ATT='B^S Real',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'bsreal',ier)
               CALL write_var_hdf5(fid,'bureal',nrho,nu,nv,ier,DBLVAR=bureal,ATT='B^U Real',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'bureal',ier)
               CALL write_var_hdf5(fid,'bvreal',nrho,nu,nv,ier,DBLVAR=bvreal,ATT='B^V Real',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'bvreal',ier)
               CALL write_var_hdf5(fid,'breal',nrho,nu,nv,ier,DBLVAR=breal,ATT='|B| Real',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'breal',ier)
               CALL close_hdf5(fid,ier)
               IF (ier /= 0) CALL handle_err(HDF5_CLOSE_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
            END IF
         CASE('FIELDLINES')
            mystart = LBOUND(R_lines,1)
            myend   = UBOUND(R_lines,1)
            CALL MPI_ALLREDUCE(myend,nlines,1,MPI_INTEGER,MPI_MAX,MPI_COMM_FIELDLINES,ierr_mpi)
            CALL torlines_write2d_parhdf5(1, nlines, 0, nsteps, mystart, myend,   'R_lines', DBLVAR=R_lines)
            CALL torlines_write2d_parhdf5(1, nlines, 0, nsteps, mystart, myend, 'PHI_lines', DBLVAR=PHI_lines)
            CALL torlines_write2d_parhdf5(1, nlines, 0, nsteps, mystart, myend,   'Z_lines', DBLVAR=Z_lines)
            CALL torlines_write2d_parhdf5(1, nlines, 0, nsteps, mystart, myend,   'B_lines', DBLVAR=B_lines)
            CALL torlines_write2d_parhdf5(1, nlines, 0, nsteps, mystart, myend,   'U_lines', DBLVAR=U_lines)
            IF (myworkid == master) THEN
               CALL open_hdf5('torlines_'//TRIM(id_string)//'.h5',fid,ier,LCREATE=.false.)
               IF (ier /= 0) CALL handle_err(HDF5_OPEN_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
               CALL write_scalar_hdf5(fid,'nlines',ier,INTVAR=nlines,ATT='Number of Fieldlines',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nlines',ier)
               CALL write_scalar_hdf5(fid,'nsteps',ier,INTVAR=nsteps+1,ATT='Number of Steps Along Fieldline',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'nsteps',ier)
               CALL write_scalar_hdf5(fid,'npoinc',ier,INTVAR=npoinc,ATT='Number of steps per field period',ATT_NAME='description')
               IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'npoinc',ier)
               CALL close_hdf5(fid,ier)
               IF (ier /= 0) CALL handle_err(HDF5_CLOSE_ERR,'torlines_'//TRIM(id_string)//'.h5',ier)
            END IF
      END SELECT
!DEC$ IF DEFINED (MPI_OPT)
      ! Make sure we wait till we're done
      CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'torlines_write',ierr_mpi)
!DEC$ ENDIF

!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE torlines_write

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
module nf95_get_missing_m

  use netcdf, only: nf90_noerr
  use nf95_get_att_m, only: nf95_get_att

  implicit none

  interface nf95_get_missing
     module procedure nf95_get_missing_real, nf95_get_missing_dble, &
          nf95_get_missing_int
  end interface nf95_get_missing

  private
  public nf95_get_missing

contains

  subroutine nf95_get_missing_real(ncid, varid, missing)

    use netcdf, only: NF90_FILL_REAL

    integer, intent(in)::  ncid, varid
    real, intent(out):: missing ! missing or fill value

    ! Local:
    integer ncerr

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

    call nf95_get_att(ncid, varid, name = "missing_value", values = missing, &
         ncerr = ncerr)

    if (ncerr /= nf90_noerr) then
       call nf95_get_att(ncid, varid, name = "_FillValue", values = missing, &
            ncerr = ncerr)
       if (ncerr /= nf90_noerr) missing = NF90_FILL_REAL
    end if

  end subroutine nf95_get_missing_real

  !**************************************************************************

  subroutine nf95_get_missing_dble(ncid, varid, missing)

    use netcdf, only: NF90_FILL_double

    integer, intent(in)::  ncid, varid
    double precision, intent(out):: missing ! missing or fill value

    ! Local:
    integer ncerr

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

    call nf95_get_att(ncid, varid, name = "missing_value", values = missing, &
         ncerr = ncerr)

    if (ncerr /= nf90_noerr) then
       call nf95_get_att(ncid, varid, name = "_FillValue", values = missing, &
            ncerr = ncerr)
       if (ncerr /= nf90_noerr) missing = NF90_FILL_double
    end if

  end subroutine nf95_get_missing_dble

  !**************************************************************************

  subroutine nf95_get_missing_int(ncid, varid, missing)

    use netcdf, only: NF90_FILL_INT

    integer, intent(in)::  ncid, varid
    integer, intent(out):: missing ! missing or fill value

    ! Local:
    integer ncerr

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

    call nf95_get_att(ncid, varid, name = "missing_value", values = missing, &
         ncerr = ncerr)

    if (ncerr /= nf90_noerr) then
       call nf95_get_att(ncid, varid, name = "_FillValue", values = missing, &
            ncerr = ncerr)
       if (ncerr /= nf90_noerr) missing = NF90_FILL_INT
    end if

  end subroutine nf95_get_missing_int

end module nf95_get_missing_m

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM


text:
      subroutine zq(imsgin,xmsgin,cmsgin,msgtype,nwds,ierror)
      implicit none
c
c#######################################################################
c
c     purpose -
c
c        this routine sets mesh quantities
c
c             ..........................................................
c             common 1 flags.
c
c             isq  - point sequence number (fitword entry)
c             imd  - model part identifer (fitword entry)
c             itp  - point type identifer (fitword entry)
c             ihb  - hydro boundary condition flag (fitword entry)
c             irb  - radiation boundary condition flag (fitword entry)
c             nn1  - nearest neighbor search flag (no neighbors are ever
c                       changed after the initial search).
c             nn2  - indicates which neighbors have been searched in a
c                       given calculational cycle.
c             icr  - constraint index to applies to this point.
c             iah  - active hydro flag (=0 ==> active, =1 ==> inactive).
c             npf  - neighbor point type flag.
c
c        format: zq/   isq/istart/iend/istep/ sequence #
c        format: zq/   imd/istart/iend/istep/ part name
c        format: zq/   itp/istart/iend/istep/ point type
c        format: zq/   ihb/istart/iend/istep/ hydro b.c.
c        format: zq/   irb/istart/iend/istep/ radiation b.c.
c        format: zq/   nn1/istart/iend/istep/ permanent n.n. flag.
c        format: zq/   nn2/istart/iend/istep/ dynamic n.n. flag.
c        format: zq/   icr/istart/iend/istep/ constraint flag.
c        format: zq/   iah/istart/iend/istep/ active hydro flag.
c        format: zq/   npf/istart/iend/istep/ neighbor point flag.
c
c             ..........................................................
c             common 2 flags.
c
c             isn  - sequence number.
c
c             ..........................................................
c             source flags.
c
c             iss  - energy source index
c             its  - temperature source index
c             ips  - pressure source index
c             ivs  - velocity (or acceleration) source index
c
c        format: zq/   iss/istart/iend/istep/ source # / start time
c        format: zq/   its/istart/iend/istep/ source # / start time
c        format: zq/   ips/istart/iend/istep/ source # / start time
c        format: zq/   ivs/istart/iend/istep/ source # / start time
c
c             ..........................................................
c             primary variables.
c
c               x  - x-coordinate for a specific mass point
c               y  - y-coordinate for a specific mass point
c               z  - z-coordinate for a specific mass point
c               u  - x-velocity for a specific mass point
c               v  - y-velocity for a specific mass point
c               w  - z-velocity for a specific mass point
c            pmat  - material pressure
c             rho  - material density
c           einter - internal energy
c             tmat - material temperature
c            trad  - radiation temperature
c            prad  - radiation pressure
c           xmass  - mass associated with a point
c          volume  - volume associated with a point
c         velocity - radial velocity relative to an arbitrary center
c           radius - radial position relative to an arbitrary center
c
c        format: zq/     x/istart/iend/istep/  x-coordinate
c        format: zq/     y/istart/iend/istep/  y-coordinate
c        format: zq/     z/istart/iend/istep/  z-coordinate
c        format: zq/     u/istart/iend/istep/  x-velocity
c        format: zq/     v/istart/iend/istep/  y-velocity
c        format: zq/     w/istart/iend/istep/  z-velocity
c        format: zq/  pmat/istart/iend/istep/  pressure
c        format: zq/   rho/istart/iend/istep/  density
c        format: zq/einter/istart/iend/istep/  internal energy
c        format: zq/  tmat/istart/iend/istep/  material temperature
c        format: zq/  trad/istart/iend/istep/  radiation temperature
c        format: zq/  prad/istart/iend/istep/  radiation pressure
c        format: zq/ xmass/istart/iend/istep/  mass
c        format: zq/volume/istart/iend/istep/  volume
c        format: zq/velocity/istart/iend/istep/velocity/xcen/ycen/zcen
c        format: zq/  radius/istart/iend/istep/  radius/xcen/ycen/zcen
c
c
c     input arguments -
c
C        imsgin()  - Integer array of command input tokens
C        xmsgin()  - Real array of command input tokens
C        cmsgin()  - Character array of command input tokens
C        msgtype() - Integer array of command input token types
C        nwds      - Number of command input tokens
c
c     output arguments -
c
C        ierror - Error Return Code (==0 ==> OK, <>0 ==> Error)
c
c     change history -
C
C        $Log: zq.f,v $
C        Revision 2.00  2007/11/09 20:04:06  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.2   Thu Feb 03 08:42:38 2000   dcg
CPVCS
CPVCS       Rev 1.1   25 Jan 2000 09:36:20   dcg
CPVCS
CPVCS       Rev 1.47   Fri Sep 03 15:49:20 1999   dcg
CPVCS    get rid of unused options
CPVCS
CPVCS       Rev 1.46   Tue Aug 24 13:35:28 1999   llt
CPVCS    testing
CPVCS
CPVCS       Rev 1.45   Fri Mar 12 12:33:10 1999   dcg
CPVCS    if nwds < 6 print out basic info if not special case
CPVCS
CPVCS       Rev 1.44   Tue Jul 21 08:29:08 1998   dcg
CPVCS    treat ipt1 as itp, icr1 as icr
CPVCS
CPVCS       Rev 1.43   Wed Dec 17 14:33:46 1997   dcg
CPVCS    remove unused cmo_get_info calls
CPVCS
CPVCS       Rev 1.42   Thu Dec 04 17:31:08 1997   dcg
CPVCS    fix bad go to if no attribute found
CPVCS
CPVCS       Rev 1.41   Mon Dec 01 16:31:56 1997   dcg
CPVCS    remove a few obsolete variables
CPVCS
CPVCS       Rev 1.40   Mon Nov 24 16:38:28 1997   dcg
CPVCS     use geom.h and calls to get_regions, get_mregions, get_surfaces
CPVCS    to access geometry data - start to isolate integer*8 dependencies
CPVCS
CPVCS       Rev 1.38   Mon Nov 17 10:39:54 1997   dcg
CPVCS    remove obsolete code
CPVCS
CPVCS       Rev 1.37   Fri Oct 31 10:51:32 1997   dcg
CPVCS    declare ipcmoprm as a pointer
CPVCS
CPVCS       Rev 1.36   Mon Apr 14 17:06:20 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.35   Mon Apr 14 16:15:38 1997   llt
CPVCS    No change.
CPVCS
CPVCS       Rev 1.34   Fri May 24 16:10:10 1996   dcg
CPVCS    fix dimension on vels, get rid of int
CPVCS
CPVCS       Rev 1.33   Sat Mar 09 10:28:38 1996   dcg
CPVCS    replace hollerith in format statements
CPVCS
CPVCS       Rev 1.32   Tue Mar 05 12:50:34 1996   dcg
CPVCS    remove int1, icn1
CPVCS
CPVCS       Rev 1.31   Mon Feb 26 13:44:54 1996   dcg
CPVCS    fixed bug with printing velocities
CPVCS
CPVCS       Rev 1.30   Fri Feb 23 16:38:10 1996   dcg
CPVCS    use dictionary to retrieve velocity names
CPVCS
CPVCS       Rev 1.29   Mon Feb 12 16:15:56 1996   dcg
CPVCS    list option works for added attributes
CPVCS    rank is now looked at for added attributes
CPVCS
CPVCS       Rev 1.28   11/20/95 09:10:08   dcg
CPVCS    test for npoints=0 then return
CPVCS
CPVCS       Rev 1.27   11/17/95 15:21:54   dcg
CPVCS    replace literal character strings in calls
CPVCS
CPVCS       Rev 1.26   11/16/95 15:23:06   dcg
CPVCS    replace character literals in calls
CPVCS
CPVCS       Rev 1.25   11/07/95 17:29:16   dcg
CPVCS    change flag to 2 in mmgetblk calls
CPVCS
CPVCS       Rev 1.24   09/29/95 09:13:00   het
CPVCS    Put in added attributes inheritance
CPVCS
CPVCS       Rev 1.23   09/19/95 08:43:48   dcg
CPVCS    allow integer type for added attribute
CPVCS
CPVCS       Rev 1.22   09/18/95 19:43:42   dcg
CPVCS    look for mesh object added attributes
CPVCS
CPVCS       Rev 1.21   08/31/95 11:54:12   dcg
CPVCS    fix default points limits (ipointi,ipointj)
CPVCS
CPVCS       Rev 1.20   08/29/95 12:15:32   dcg
CPVCS    set length for names to 40 characters
CPVCS
CPVCS       Rev 1.19   08/23/95 06:59:40   het
CPVCS    Remove the CMO prefix from SB-ids
CPVCS
CPVCS       Rev 1.18   08/22/95 06:51:36   het
CPVCS    Split the storage block for CMO variables.
CPVCS
CPVCS       Rev 1.17   07/17/95 16:00:24   dcg
CPVCS    use names for point types
CPVCS
CPVCS       Rev 1.16   07/15/95 02:21:20   het
CPVCS    Correct an error with ipointi and ipointj
CPVCS
CPVCS       Rev 1.15   06/27/95 16:38:48   dcg
CPVCS    remove second literal argument in memory management calls
CPVCS
CPVCS       Rev 1.14   06/20/95 15:41:52   dcg
CPVCS    remove character literals from arguments list to hgetprt
CPVCS
CPVCS       Rev 1.13   06/19/95 16:43:58   dcg
CPVCS    add blank after literal in calling sequence to savpart
CPVCS
CPVCS       Rev 1.12   06/13/95 09:04:22   ejl
CPVCS    Cleaned up msgtty, calling arguments.
CPVCS
CPVCS
CPVCS       Rev 1.11   06/07/95 15:32:34   het
CPVCS    Change character*32 idsb to character*132 idsb
CPVCS
CPVCS       Rev 1.10   05/15/95 13:36:56   het
CPVCS    Make changes to the regset and surfset routines
CPVCS
CPVCS       Rev 1.9   05/01/95 08:37:00   het
CPVCS    Modifiy all the storage block calles for long names
CPVCS
CPVCS       Rev 1.8   03/31/95 09:11:30   het
CPVCS    Add the buildid calles before all storage block calls
CPVCS
CPVCS       Rev 1.7   03/30/95 05:01:02   het
CPVCS    Change the storage block id packing and preidsb to buildid for long names
CPVCS
CPVCS       Rev 1.6   03/23/95 22:59:44   het
CPVCS    Add the model routines and add the cmo name into the idsbs
CPVCS
CPVCS       Rev 1.5   03/23/95 15:08:54   dcg
CPVCS     Add mesh object name to storage block id for surface,region info.
CPVCS
CPVCS       Rev 1.4   03/10/95 17:12:50   dcg
CPVCS     put in mesh object calls
CPVCS
CPVCS       Rev 1.3   02/18/95 06:57:38   het
CPVCS    Changed the parameter list to be the same as pntlimc
CPVCS
CPVCS       Rev 1.2   01/04/95 22:06:46   llt
CPVCS    unicos changes (made by het)
CPVCS
CPVCS       Rev 1.1   12/19/94 08:27:34   het
CPVCS    Add the "comdict.h" include file.
CPVCS
CPVCS
CPVCS       Rev 1.0   11/10/94 12:20:42   pvcs
CPVCS    Original version.
C
c
c#######################################################################
c
c
c#######################################################################
c
      integer nplen
      parameter (nplen=1000000)
c
      character*132 logmess

c
      pointer (ipisetwd, isetwd)
      pointer (ipimt1, imt1)
      pointer (ipitp1, itp1)
      pointer (ipicr1, icr1)
      pointer (ipisn1, isn1)
      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      integer isetwd(nplen)
      integer imt1(nplen), itp1(nplen),
     *        icr1(nplen), isn1(nplen)
      real*8 xic(nplen), yic(nplen), zic(nplen)
      pointer ( ipmpary1 , mpary1(1) )
      integer mpary1

c
      pointer ( ipies , ies1(1) )
      pointer ( iptmes , tmes1(1) )
      pointer ( ipits , its1(1) )
      pointer ( iptmts , tmts1(1) )
      pointer ( ipips , ips1(1) )
      pointer ( iptmps , tmps1(1) )
      pointer ( ipivs , ivs1(1) )
      pointer ( iptmvs , tmvs1(1) )
      pointer (ip,out)
      real*8 out(*)
c
      integer nwds, imsgin(nwds), msgtype(nwds)
      REAL*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
C
      integer ierror
C
      character*32 isubname,cnamz
      character*32 ich1,ich2,ich3, cmo, crank, ctype

      pointer(ipxfield,xfield)
      real*8 xfield(1000000)
      pointer(ipxfield,ifield)
      pointer (ipvels,vels)
      real*8 vels(3,1000000)
      integer ifield(1000000)
C
      character*32 iword, iword1, cvelnm
      character*8 cpart
      character*132 iformat
      real*8 tmdvs,tmdss,uradius,radius,emat,rho,pmat,prad,tmat,trad,
     *  wa,va,ua,rout,xa,ya,za,tmdps,tmdts,tmvs1,tmps1,tmts1,
     * tmes1
      integer ivs,ips,its,iss,nptstmp,iout,isn,nn2,nn1,irb,ihb,itp,imt,
     * isq,ier2,ier,lin,itin,ier1,irank,ierror_return,j1,i1,len,
     * icharlnf,icr,if1,if2,i,mpno,ipt1,ipt2,ipt3,itype,icscode,ipointi,
     * ipointj,mmrel2,iunpk,iprt5,iprt4,iprt1,iprt2,iprt3,ihead1,
     * ihead2,ihead3,ihead4,ihead5,ierr,ityp,icmotype,length,npoints,
     * ierrw,ivs1,ips1,l,iatt_index,ilen,ntets,ies1,its1
      character*32 cio,cpers,cinter,clen

c
c#######################################################################
c
c
c     ******************************************************************
c
c     set the memory management path name to be the subroutine name.
c
      isubname='zq'
      cnamz='zq2'
      cpart='part'
      isubname='zq'
c
c  get mesh object
      call cmo_get_name(cmo,ierror)
      if(ierror.ne.0) then
        write(logmess,'(a)') 'ZQ found bad mesh object'
        call writloga('default',0,logmess,0,ierrw)
        go to 9999
      endif
C
      call cmo_get_info('nnodes',cmo,
     *                  npoints,length,icmotype,ierror)
c  check that there are points -- if none return
      if(npoints.eq.0) then
        write(logmess,'(a)') 'npoints = 0 in subroutine z'
        call writloga('default',0,logmess,0,ierrw)
        go to 9999
      endif
      call cmo_get_info('nelements',cmo,
     *                  ntets,length,icmotype,ierror)
      call cmo_get_info('isetwd',cmo,
     *                  ipisetwd,ilen,ityp,ierr)
      call cmo_get_info('imt1',cmo,ipimt1,ilen,ityp,ierr)
      call cmo_get_info('itp1',cmo,ipitp1,ilen,ityp,ierr)
      call cmo_get_info('icr1',cmo,ipicr1,ilen,ityp,ierr)
      call cmo_get_info('isn1',cmo,ipisn1,ilen,ityp,ierr)
      call cmo_get_info('xic',cmo,ipxic,ilen,ityp,ierr)
      call cmo_get_info('yic',cmo,ipyic,ilen,ityp,ierr)
      call cmo_get_info('zic',cmo,ipzic,ilen,ityp,ierr)
c     initialize some variables.
c
      ihead1=0
      ihead2=0
      ihead3=0
      ihead4=0
      ihead5=0
      iprt1=0
      iprt2=0
      iprt3=0
      iprt4=0
      iprt5=0
      iunpk=0
      mmrel2=0
c
c     ******************************************************************
c
c     set the point index boundaries.
c
      ipointi=0
      ipointj=0
      call cmo_get_info('ipointi',cmo,ipointi,ilen,itype,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'cmo_get_info')
      call cmo_get_info('ipointj',cmo,ipointj,ilen,itype,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'cmo_get_info')
C
      if(msgtype(3).eq.1.and.imsgin(3).eq.0) then
         imsgin(3)=max(1,ipointi)
      endif
      if(msgtype(4).eq.1.and.imsgin(4).eq.0) then
         if(ipointj.le.0.or.imsgin(3).eq.1) then
            imsgin(4)=npoints
         else
            imsgin(4)=ipointj
         endif
      endif
      if(msgtype(5).eq.1.and.imsgin(5).eq.0) then
         imsgin(5)=1
      endif
C
      ich1=' '
      ich2=' '
      ich3=' '
      if(msgtype(3).eq.1) then
        ipt1=max(1,min(imsgin(3),npoints))
        ipt2=max(1,min(imsgin(4),npoints))
        ipt3=max(1,min(imsgin(5),npoints))
      else
        ich1=cmsgin(3)
        ich2=cmsgin(4)
        ich3=cmsgin(5)
      endif
c
c     ******************************************************************
c
c     check point limits and translate to valid limits if necessary.
c
      length=npoints
      call mmgetblk('mpary1',isubname,ipmpary1,length,2,icscode)
      if(msgtype(3).eq.1) then
         call pntlimn(ipt1,ipt2,ipt3,ipmpary1,mpno,npoints,isetwd,itp1)
      elseif(msgtype(3).ne.1) then
         call pntlimc(ich1,ich2,ich3,ipmpary1,mpno,npoints,isetwd,itp1)
      endif
c
c     ******************************************************************
c
c     set requested values for mass points.
c
      do 100 j1=1,mpno
         i1=mpary1(j1)
         iword=cmsgin(2)
         len=icharlnf(iword)
c
c        _______________________________________________________________
c
c        write data for points.
c
         if(iword(1:len).eq.'all') then
            iprt1=1
            iprt2=1
            iprt3=1
            goto 9995
         endif
         if(iword(1:len).eq.'imt') then
            len=len+1
            iword(1:len)='imt1'
         endif
         if(iword(1:len).eq.'isn') then
            len=len+1
            iword(1:len)='isn1'
         endif
c
c
c        set type in common.
c
         if(iword(1:len).eq.'itp'.or.iword(1:len).eq.'itp1') then
            ierror=0
            if(nwds.lt.6) then
               iprt1=1
               goto 9995
            endif
            if(msgtype(6).eq.1) then
              itp1(i1)=imsgin(6)
            elseif (msgtype(6).eq.3) then
               call getptyp(cmsgin(6),itp1(i1),ierror)
               if (ierror.ne.0) call x3d_error(isubname,'getptyp')
            else
               write(logmess,'(a30)') 'illegal argument in zq command'
               call writloga ('default',0,logmess,0,ierr)
            endif
            goto 9998
         endif
c
c
c        set constriants index in common.
c
         if(iword(1:len).eq.'icr'.or.iword(1:len).eq.'icr1') then
            ierror=0
            if(nwds.lt.6) then
               iprt1=1
               goto 9995
            endif
            if(msgtype(6).le.2) then
               icr=imsgin(6)
            elseif(msgtype(6).eq.3) then
               iword=cmsgin(6)
               if1=5
               if2=8
               do 120 i=if1,if2
                  if(iword(i:i).eq.' ') then
                     iformat=' '
                     write(iformat,"('(a',i1,',i',i1,')')") if1-1,i-if1
                     goto 121
                  endif
 120           continue
 121           continue
               read(iword,iformat) iword1,icr
            endif
            icr1(i1)=imsgin(6)
            goto 9998
         endif
c
c        _______________________________________________________________
c
c        set xyz-coordinates.
c
         if(iword(1:len).eq.'xyz') then
            ierror=0
            if(nwds.lt.6) then
               iprt2=1
               goto 9995
            endif
            xic(i1)=xmsgin(6)
            yic(i1)=xmsgin(7)
            zic(i1)=xmsgin(8)
            goto 9998
         endif
c
c        _______________________________________________________________
c
c        set x-coordinate.
c
         if(iword(1:len).eq.'x'.or.iword(1:len).eq.'xic') then
            ierror=0
            if(nwds.lt.6) then
               iprt2=1
               goto 9995
            endif
            xic(i1)=xmsgin(6)
            goto 9998
         endif
c
c        _______________________________________________________________
c
c        set y-coordinate.
c
         if(iword(1:len).eq.'y'.or.iword(1:len).eq.'yic') then
            ierror=0
            if(nwds.lt.6) then
               iprt2=1
               goto 9995
            endif
            yic(i1)=xmsgin(6)
            goto 9998
         endif
c
c        _______________________________________________________________
c
c        set z-coordinate.
c
         if(iword(1:len).eq.'z'.or.iword(1:len).eq.'zic') then
            ierror=0
            if(nwds.lt.6) then
               iprt2=1
               goto 9995
            endif
            zic(i1)=xmsgin(6)
            goto 9998
         endif
c        _______________________________________________________________
c      look for added attribute
c
         if(nwds.le.5) then
           iprt1=1
           iprt2=1
           go to 9995
         endif
         if(j1.eq.1) then
            call  cmo_get_attparam(iword,cmo,iatt_index,ctype,crank,
     *      clen,cinter,cpers,cio,ierror_return)
            if(icscode.eq.0) then
C   found existing attribute
               if(ctype(1:4).eq.'VINT') ityp=1
               if(ctype(1:7).eq.'VDOUBLE') ityp=2
               call mmgetpr(iword(1:len),cmo,ipxfield,icscode)
               call cmo_get_info(crank,cmo,irank,lin,itin,ier1)
            else
               write(logmess,9005) iword(1:len),cmo
               call writloga('default',0,logmess,0,ierrw)
               go to 9998
            endif
         endif
         ierror=0
         if((j1.eq.1.and.icscode.eq.0.and.ier.eq.0.and.ier1.eq.0.and.
     *            ier2.eq.0).or.j1.ne.1) then
            if(msgtype(6).eq.2) then
               do l=1,irank
                  xfield((i1-1)*irank+l)=xmsgin(6)
               enddo
            elseif(msgtype(6).eq.1) then
               do l=1,irank
                  ifield((i1-1)*irank+l)=imsgin(6)
               enddo
            elseif(nwds.lt.6) then
               do l=1,irank
                  if(ityp.eq.1) write(logmess,'(i6,i6)') i1,
     *                 ifield((i1-1)*irank+l)
                  if(ityp.eq.2) write(logmess,'(i6,e10.3)') i1,
     *                 xfield((i1-1)*irank+l)
                  call writloga('default',0,logmess,0,ierrw)
               enddo
            endif
         else
            write(logmess,9005) iword(1:len),cmo
 9005       format(' zq cannot find attribute ',a,' in ',a)
            call writloga('default',0,logmess,0,ierrw)
            go to 9997
         endif
c        _______________________________________________________________
c
c        invalid option used if this location reached.
c
         goto 9998
c
c        _______________________________________________________________
c
c        write point data to log files.
c
 9995    continue
c
         ierror=0
         if(ihead1.eq.0) then
            ihead1=1
            if(iprt1.eq.1) then
               write(logmess,9000)
               call writloga('default',0,logmess,0,ierr)
            endif
         endif
         if(ihead2.eq.0) then
            ihead2=1
            if(iprt2.eq.1) then
               write(logmess,9030)
               call writloga('default',0,logmess,0,ierr)
            endif
         endif
         if(ihead3.eq.0) then
            ihead3=1
            if(iprt3.eq.1) then
               write(logmess,9050)
               call writloga('default',0,logmess,0,ierr)
            endif
         endif
         if(ihead4.eq.0) then
            ihead4=1
            if(iprt4.eq.1) then
               write(logmess,9070)
               call writloga('default',0,logmess,0,ierr)
            endif
         endif
         isq=0
         imt=imt1(i1)
         itp=itp1(i1)
         ihb=0
         irb=0
         nn1=0
         nn2=0
Cdcg     int=int1(i1)
         icr=icr1(i1)
Cdcg     icn=icn1(i1)
         isn=isn1(i1)
c
         if(iprt1.eq.1) then
            write(logmess,9010) i1,isq,imt,itp,icr,isn
            call writloga('default',0,logmess,0,ierr)
         endif
         if(iprt2.eq.1) then
            xa=xic(i1)
            ya=yic(i1)
            za=zic(i1)
            if(j1.eq.1)
     *        call cmo_get_attinfo('velname',cmo,iout,rout,cvelnm,
     *         ip,length,itype,ier)
            if(ier.eq.0.and.j1.eq.1)
     *        call cmo_get_info(cvelnm,cmo,ipvels,length,itype,ier)
            if(ier.eq.0) then
               ua=vels(1,i1)
               va=vels(2,i1)
               wa=vels(3,i1)
            else
               ua=0.
               va=0.
               wa=0.
            endif
            write(logmess,9040) i1,xa,ya,za,ua,va,wa
            call writloga('default',0,logmess,0,ierr)
         endif
Cdcg     pmat=pic(i1)
Cdcg     rho=ric(i1)
Cdcg     emat=eic(i1)
         tmat=0.0
         trad=0.0
         prad=0.0
         if(iprt3.eq.1) then
            write(logmess,9040) i1,rho,emat,tmat,trad,pmat,prad
            call writloga('default',0,logmess,0,ierr)
         endif
c
         if(iprt4.eq.1) then
            radius=sqrt(xic(i1)*xic(i1)+yic(i1)*yic(i1)
     *              +zic(i1)*zic(i1) )
            call cmo_get_attinfo('velname',cmo,iout,rout,cvelnm,
     *       ip,length,itype,ier)
            call cmo_get_info(cvelnm,cmo,ipvels,length,itype,ier)
            if(ier.eq.0) then
               uradius=
     *          sqrt(vels(1,i1)*vels(1,i1)+vels(2,i1)*vels(2,i1)
     *              +vels(3,i1)*vels(3,i1) )
            else
               uradius=0.
            endif
            write(logmess,9040) i1,radius,uradius
            call writloga('default',0,logmess,0,ierr)
         endif

        if(iprt5.eq.1) then
           if(nptstmp.le.0) then
              write(logmess,9082)
              call writloga('default',0,logmess,0,ierr)
9082          format('no source points have been assigned.')
           else
             iss=ies1(i1)
              tmdss=tmes1(i1)
              its=its1(i1)
              tmdts=tmts1(i1)
              ips=ips1(i1)
              tmdps=tmps1(i1)
              ivs=ivs1(i1)
              tmdvs=tmvs1(i1)
              write(logmess,9080) i1,iss,tmdss,its,tmdts,ips,tmdps,ivs,
     *                            tmdvs
              call writloga('default',0,logmess,0,ierr)
9080          format(i5,2x,4(i4,2x,1pe12.4))
           endif
        endif
 9000    format('      ipt',' isq',' imt',' itp',' icr',' isn')
 9010    format(i9,4i4,i9)
 9030    format(1x,'ipt',9x,'xa',10x,'ya',10x,'za',
     *                         10x,'ua',10x,'va',10x,'wa')
 9040    format(i5,2x,6(1pe12.4))
 9050    format(1x,'ipt',8x,'rho ',8x,'emat',8x,'tmat',
     *                         8x,'trad',8x,'pmat',8x,'prad')
 9070    format(1x,'ipt',5x,'radius',4x,'velocity')
         goto 9998
 9998    continue
c
c        _______________________________________________________________
c
 100  continue
c
c     ******************************************************************
c
c     release any local memory allocated for this routine.
c
 9997 call mmrelblk('mpary1',isubname,ipmpary1,icscode)
c
c     if (mmrel2 .ne. 0) call mmrelprt(cnamz, ics)
c
c
c     ******************************************************************
c
c     error returns transfer to this statement 9999
c
      goto 9999
 9999 continue
c
c     ******************************************************************
c
      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
!##############################################################################
!# ****************************************************************************
!# <name> AllenCahn </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains the actual nonstationary solver for the Allen Cahn
!# problem, i.e. the timeloop and time stepping. The following routines can
!# be found here:
!#
!# 1.) AC5_initparameters
!#     -> Initialise the parameters for the time stepping.
!#
!# 2.) AC5_timestep
!#     -> Calculate the solution of the next timestep.
!#
!# 3.) AC5_postprocessing
!#     -> Perform postprocessing (write GMV's,...)
!# </purpose>
!##############################################################################

module AllenCahn

  use fsystem
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use sortstrategy
  use coarsegridcorrection
  use ucd
  use timestepping
  use genoutput

  use collection
  use paramlist

  use AllenCahn_callback
  use AllenCahn_basic
  use AllenCahn_nonlinearcore

  use AllenCahn_matvec
  use AllenCahn_boundarycondition
  use AllenCahn_partridiscr
  use AllenCahn_solver

  IMPLICIT NONE

CONTAINS
  ! ***************************************************************************
!<subroutine>

  subroutine AC5_timestep (rACproblem,rACvector,rACrhs,rNSproblem,rNSvector,rNSrhs)

!<description>
  ! Performs one time step: $t^n -> t^n+1$.
  ! Assembles system matrix and RHS vector.
  ! Solves the corresponding time-step equation and returns the solution vector
  ! at the end of the time step.
  ! Solves the given problem by applying a linear solver.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem

  ! The current solution vector at time $t^n$. Is replaced by the
  ! solution vector at time $t^{n+1}.
  type(t_vectorBlock), INTENT(INOUT) :: rACvector

  ! The RHS vector at time $t^n$. Is replaced by the RHS at time $t^{n+1}$.
  type(t_vectorBlock), INTENT(INOUT) :: rACrhs

!Mcai
!~~~~~~~~~~~~~~~~~~~~~NS problem~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  type(t_problem), intent(inout), optional :: rNSproblem
  type(t_vectorBlock), intent(inout), optional :: rNSvector
  type(t_vectorBlock), intent(inout), optional :: rNSrhs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!</inputoutput>
!</subroutine>

  ! local variables
    INTEGER :: NLMIN,NLMAX
    INTEGER :: i

    ! Error indicator during initialisation of the solver
    INTEGER :: ierror

    ! A filter chain to filter the vectors and the matrix during the
    ! solution process.
    type(t_filterChain), DIMENSION(1), TARGET :: RfilterChain
    type(t_filterChain), DIMENSION(:), POINTER :: p_RfilterChain

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_matrixBlock), POINTER :: p_rmatrix
    type(t_vectorBlock), POINTER :: p_rrhs
    type(t_vectorBlock)  :: rtempBlock

    ! A solver node that accepts parameters for the linear solver
    type(t_linsolNode), POINTER :: p_rsolverNode,p_rsmoother
    type(t_linsolNode), POINTER :: p_rcoarseGridSolver,p_rpreconditioner

    ! An array for the system matrix(matrices) during the initialisation of
    ! the linear solver.
    type(t_matrixBlock), DIMENSION(1:rACproblem%NLMAX) :: Rmatrices

    ! An interlevel projection structure for changing levels
    type(t_interlevelProjectionBlock) :: rprojection

    ! One level of multigrid
    type(t_linsolMGLevelInfo), POINTER :: p_rlevelInfo

!~~~~~~~~~~~~~~local variables for AC problem~~~~~~~~~~~~~~~~~~~~~~~~~~
    !local variables for AC problem
    type(t_ACNonlinearIteration) :: rACnonlinearIteration
	! A tmp rnonlinearIterarion

    type(t_ACNonlinearIteration) :: rnonlinearIterationTmp
    ! The nonlinear solver configuration
    type(t_nlsolNode) :: rACnlSol

    ! Debug
    real(DP), dimension(:), pointer ::  p_vectordata

! Mcai
	! A parameter to determine, whether we use implicit scheme or explicit scheme for
	! Allen-Cahn problem, if it is 0, we treat convective term explictly. Otherwise,implicitly
	integer :: Implicit_Convective=0

!~~~~~~~~~~~~~~If we only have Allen-Cahn problem, activate the following one~~~~~~~~
    !MCai, For testing the Allen-Cahn solver only, we set rNSvector = 0 vector
	call lsysbl_clearVector(rNSvector)

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    ! We have an equation of the type
    !
    !   d/dt u(x,t)  +  N(u(x,t))  =  f(x,t)
    !
    ! Which is discretised in time with a Theta scheme, leading to
    !
    !   $$ u_{n+1} + w_1*N(u_n+1)
    !      =   u_n + w_2*N(u_n)  +  w_3*f_{n+1}  +  w_4*f_n $$
    !
    ! with k=time step size, u_{n+1} = u(.,t_{n+1}),etc., c.f. timestepping.f90.
    !
    ! The RHS of that equation therefore contains parts of the solution
    ! u_n, of the old RHS f_n and the new RHS f_{n+1}. At first, we make
    ! a weighted copy of the current RHS f_n to the 'global' RHS vector
    ! according to the time stepping scheme.

    NLMIN = rACproblem%NLMIN
    NLMAX = rACproblem%NLMAX

! MCai~~~~~~~Geneate Poly matrix and Conv matrix on all levels~~~~~~~~~~~~~~~
    ! We first generate the matrix based on nonlinear term
!    call AC_assembleMatPoly(rACproblem, rACvector)

!~~~~~~~~~MCai, implicitly treat convection term?~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~think about how to implement this	*)*)*)*)*)*)*)**)*)*)*)*)*)*)*)*)*)*)
! we may also need to update the system matrix of AC problem by using new obtained
! velocity field in NS problem.

    ! We can also generate convetive matrix here
    call AC_assembleMatConv(rACproblem, rACvector, rNSproblem, rNSvector)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~We need to specify rACnonlinearIteration, rACnlSolver~~~~~~~
    call AC_getNonlinearSolver (rACnlSol)

    ! Initialise the nonlinear loop. This is to prepare everything for
    ! our callback routines that are called from the nonlinear solver.
    ! The preconditioner in that structure is initialised later.
    call AC_initNonlinearLoop (rACproblem,rACproblem%NLMIN,rACproblem%NLMAX,&
        rACvector,rACrhs,rACnonlinearIteration)

!~~~~~~~~~~~~~~~call nonlinear solver here~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Set up both RHS and preconditioner for core equation:

    ! First, p_rrhs is to save RHS for core equation
    p_rrhs => rACproblem%rrhs
    call lsysbl_assignDiscrIndirect(rACproblem%rrhs, p_rrhs)

    ! Create a temporary vector we need for some preparations.
    call lsysbl_createVecBlockIndirect (p_rrhs, rtempBlock, .FALSE.)

    ! Set up w_2*N(u_n) + w_4*f_n.
	! First, set up w_4*f_n
    call lsysbl_vectorLinearComb(rACrhs,p_rrhs,&
         rACproblem%rtimedependence%rtimestepping%dweightOldRHS,0.0_DP)

    ! Set up w_2*N(u_n). Unlike cc2d, we do not need rnonlinearACmatrix...
    ! We pass those coefficient in front of matrices through rACnonlinearIteration
    ! dalpha is for Mass matrix
    rACnonlinearIteration%dalpha=-1.0_DP
	! dtheta is for Laplace term
    rACnonlinearIteration%dtheta=-rACproblem%rtimedependence%rtimestepping%dweightMatrixRHS
	! dconv is for the convetive term
    rACnonlinearIteration%dconv=-rACproblem%rtimedependence%rtimestepping%dweightMatrixRHS
	! dgamma is for nonlinear term
    rACnonlinearIteration%dgamma=-rACproblem%rtimedependence%rtimestepping%dweightMatrixRHS

    call AC_nonlinearMatMul (rACproblem,rACvector,rACnonlinearIteration, &
                                           p_rrhs, -1.0_DP,1.0_DP)

    ! Switch to the next point in time.
    rACproblem%rtimedependence%dtime = rACproblem%rtimedependence%dtime + &
          rACproblem%rtimedependence%rtimestepping%dtstep

!~~~~~~~Mcai, we need to rewrite AC5_calRHS, because RHS may depending on rACvector
    ! Generate f_n+1 into the rACrhs overwriting the previous RHS.
	call AC5_calcRHS (rACproblem, rACvector, rACrhs,&
                          rNSproblem, rNSvector, rNSrhs)
! Debugging, to make sure that the integration of coeff_RHS3 is correct.
!     call lsyssc_getbaseVector_double(p_rrhs%rvectorBlock(1), p_vectordata)
!     print *, p_vectordata(1)
!     print *, p_vectordata(10)
!     print *, p_vectordata(100)
!     stop
    ! Add w_3*f_{n+1} to the p_rrhs
    call lsysbl_vectorLinearComb(rACrhs,p_rrhs,&
         rACproblem%rtimedependence%rtimestepping%dweightOldRHS,1.0_DP)

!     call AC5_calcRHS (rACproblem, rACvector, rACrhs,&
!                           rNSproblem, rNSvector, rNSrhs)
!
!    !print *, rACproblem%rtimedependence%rtimestepping%dweightOldRHS
!
!     call lsysbl_vectorLinearComb(rACrhs,p_rrhs,&
!          rACproblem%rtimedependence%rtimestepping%dweightOldRHS,1.0_DP)

    ! Discretise the boundary conditions at the new time instant
    call AC5_initDiscreteBC (rACproblem)

    ! Implement boundary conditions into the RHS vector, the solution vector
    ! and the current system matrices.
    call AC5_implementBC (rACproblem,rACvector,p_rrhs,1.0_DP)

    ! That's it for the RHS vector. (saved in p_rrhs)
    !
    ! The LHS "u_{n+1} + w_1*N(u_n+1)" results in the system matrix
    ! "M + w_1 N(.)" for the next linear system to solve. Set up that system
    ! on every level of the discretisation.

    ! we can use rnonlinearIterationTmp for assembling the matrix of LHS, or we directly
	! assemble system matrix for all levels and save them in rACproblem%RlevelInfo(i)%rmatrix
!    rnonlinearIterationTmp = rACnonlinearIteration
!	! A = \alpha M + \theta Laplace + \conv Conv + \gamma Nonlinear
    rnonlinearIterationTmp=rACnonlinearIteration
    rnonlinearIterationTmp%dalpha = 1.0_DP
    rnonlinearIterationTmp%dtheta = rACproblem%rtimedependence%rtimestepping%dweightMatrixLHS
    rnonlinearIterationTmp%dconv = rACproblem%rtimedependence%rtimestepping%dweightMatrixLHS
    rnonlinearIterationTmp%dgamma = rACproblem%rtimedependence%rtimestepping%dweightMatrixLHS

! The convection matrix is passed through rACproblem%RlevelInfo(i)%rmatrixConv
	!Preparing the system matrix for all levels, since we are going to use MG
!MCai, do we need it here?
    call AC_updatePreconditioner (rACproblem,rnonlinearIterationTmp,&
                                  rACvector)

! MCai, make sure that rnonlinearIterationTmp%dtheta is implemented correctly

! We have different ways to define p_rsolverNode, one is using the following
!    p_rsolverNode => rACproblem%p_rsolverNode
! The other is to use
!    p_rsolverNode => rACnonlinearIteration%p_rsolverNode
! remember to call
!    call linsol_doneData (p_rsolverNode)

    !print *, 'before core eqn, we have assembled system matrices and p_rsolverNode, saved in rnonlinearIterationTmp'
    call AC_solveCoreEquation (rACproblem,rnonlinearIterationTmp,rACnlSol,&
      rACvector,p_rrhs,rtempBlock)

    ! We need to clear system matrix
	do i= NLMIN, NLMAX
           call lsysbl_clearMatrix(rACproblem%RlevelInfo(i)%rmatrix)
	end do

    ! Release the temporary vector
    call lsysbl_clearVector(p_rrhs)
    call lsysbl_releaseVector (rtempBlock)
    call lsyssc_releaseVector(rnonlinearIterationTmp%p_rtempVectorSc)
    call lsyssc_releaseVector(rnonlinearIterationTmp%p_rtempVectorSc2)
    call linsol_doneData (rnonlinearIterationTmp%p_rsolverNode)
    ! there are memory problem,how to solve this
!     call linsol_doneStructure (rnonlinearIterationTmp%p_rsolverNode)
!     call linsol_releaseSolver (rnonlinearIterationTmp%p_rsolverNode)

    ! Finally tell the time stepping scheme that we completed the time step.
    call timstp_nextSubstep (rACproblem%rtimedependence%rtimestepping)

  end subroutine

  ! ***************************************************************************
!<subroutine>
  subroutine AC5_initparameters (rparams,rACproblem)

!<description>
  ! Reads the DAT file from disc into the parameter list rparams and
  ! initialises basic variables (number of levels, time stepping technique)
  ! in rACproblem according to these settings.
!</description>

!<inputoutput>
  ! A parameter list structure accepting the parameters from the DAT file.
  type(t_parlist), INTENT(INOUT) :: rparams

  ! A problem structure saving problem-dependent information.
  type(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem
!</inputoutput>

!</subroutine>

  ! local variables
  INTEGER :: cscheme,niterations
  real(DP) :: dtheta,dtstep,dtimemin,dtimemax

    ! Read the parameters from disc and put a reference to it
    ! to the collection
!    call parlst_readfromfile(rparams, 'data/AllenCahn.dat')

    ! Get the parameters for the time stepping scheme from the parameter list
    call parlst_readfromfile(rparams, 'data/timediscr.dat')
    call parlst_getvalue_int (rparams, 'TIME-DISCRETISATION', 'itimeStepScheme', cscheme, 0)
    call parlst_getvalue_int (rparams, 'TIME-DISCRETISATION ', 'NITERATIONS', niterations, 1000)
    call parlst_getvalue_double (rparams, 'TIME-DISCRETISATION', 'DTIMESTEPTHETA', dtheta, 1.0_DP)
    call parlst_getvalue_double (rparams, 'TIME-DISCRETISATION', 'dtimeStep', dtstep, 0.1_DP)
    call parlst_getvalue_double (rparams, 'TIME-DISCRETISATION', 'DTIMEINIT', dtimemin, 0.0_DP)
    call parlst_getvalue_double (rparams, 'TIME-DISCRETISATION', 'DTIMEMAX', dtimemax, 1.0_DP)

    call  parlst_readfromfile(rparams, 'data/discretisation.dat')
    ! We want to solve our Laplace problem on level...
    call parlst_getvalue_int (rparams, 'CC-DISCRETISATION', 'NLMIN', rACproblem%NLMIN, 7)
    call parlst_getvalue_int (rparams, 'CC-DISCRETISATION', 'NLMAX', rACproblem%NLMAX, 7)

    ! Initialise the time stepping in the problem structure
    call timstp_init (rACproblem%rtimedependence%rtimestepping, &
                      cscheme, dtimemin, dtstep, dtheta)

    rACproblem%rtimedependence%niterations = niterations

    rACproblem%rtimedependence%dtimemin = dtimemin
    rACproblem%rtimedependence%dtime = dtimemin
    rACproblem%rtimedependence%dtimemax = dtimemax

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine AC5_postprocessing (rACproblem,rvector,iiteration,dtime)

!<description>
  ! Writes the solution into a GMV file.
!</description>

!<input>
  ! Number of current iteration
  INTEGER, INTENT(IN) :: iiteration

  ! Current simulation time
  real(DP), INTENT(IN) :: dtime
!</input>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem

  ! The current solution vector.
  type(t_vectorBlock), INTENT(INOUT) :: rvector
!</inputoutput>

!</subroutine>

  ! local variables
     ! file name for output
     character(SYS_STRLEN) :: sfile,sfilename

    ! We need some more variables for postprocessing
    real(DP), DIMENSION(:), POINTER :: p_Ddata

    ! Output block for UCD output to GMV file
    type(t_ucdExport) :: rexport

    ! A pointer to the solution vector and to the triangulation.
    type(t_triangulation), POINTER :: p_rtriangulation

! Debug
!    write(*,*)'time dependent'
!    print *, dtime
!    print *, iiteration

    ! Basic filename
    sfilename='AC5_phi.gmv'

!    SFILENAMEUCD='gmv/AC5_phi.gmv'
    sfile = trim(adjustl(sfilename))//'.'//sys_si0(iiteration,5)
!    sfile = trim(adjustl(sfilename))//'.'//sys_si0(iiteration,5)
!    sfile = trim(adjustl(sfilename))//'.'//sys_si0(iiteration,5)

    ! From the attached discretisation, get the underlying triangulation
    p_rtriangulation => &
      rvector%RvectorBlock(1)%p_rspatialDiscr%p_rtriangulation

    ! p_rvector now contains our solution. We can now
    ! start the postprocessing.
    ! Start UCD export to GMV file:

!    call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,&
!                       'gmv/ACu5.gmv.'//TRIM(sys_si0L(iiteration,5)))

    call ucd_startGMV (rexport,UCD_FLAG_STANDARD,p_rtriangulation,&
                       'gmv/'//sfile//TRIM(sys_si0L(iiteration,5)))

!    call ucd_startVTK (rexport,UCD_FLAG_STANDARD,p_rtriangulation,sfile)

    call ucd_setSimulationTime (rexport,dtime)

    call lsyssc_getbase_double (rvector%RvectorBlock(1),p_Ddata)
    call ucd_addVariableVertexBased (rexport,'sol',UCD_VAR_STANDARD, p_Ddata)

    ! Write the file to disc, that's it.
    call ucd_write (rexport)
    call ucd_release (rexport)

!    if (present(dtime)) then
      ! Update time stamp of last written out GMV.
!      rpostprocessing%dnextTimeUCD = rpostprocessing%dnextTimeUCD+dtimeDifferenceUCD
!      rpostprocessing%inextFileSuffixUCD = rpostprocessing%inextFileSuffixUCD + 1
!       iiteration=iiteration+1
!    end if

  end subroutine

  ! ***************************************************************************

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
-makelib ies/xil_defaultlib -sv \
  "C:/Xilinx/Vivado/2017.2/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
  "C:/Xilinx/Vivado/2017.2/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib ies/xpm \
  "C:/Xilinx/Vivado/2017.2/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies/axi_infrastructure_v1_1_0 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/7e3a/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib ies/xil_common_vip_v1_0_0 -sv \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/2ad9/hdl/xil_common_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib ies/smartconnect_v1_0 -sv \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/39ca/hdl/sc_util_v1_0_vl_rfs.sv" \
-endlib
-makelib ies/axi_protocol_checker_v1_1_14 -sv \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/a1b2/hdl/axi_protocol_checker_v1_1_vl_rfs.sv" \
-endlib
-makelib ies/axi_vip_v1_0_2 -sv \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/6eb1/hdl/axi_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib ies/axi_vip_v1_0_1 -sv \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/1d61/hdl/processing_system7_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_processing_system7_0_0/sim/zynq_design_1_processing_system7_0_0.v" \
-endlib
-makelib ies/axi_lite_ipif_v3_0_4 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/0ba0/hdl/axi_lite_ipif_v3_0_vh_rfs.vhd" \
-endlib
-makelib ies/lib_cdc_v1_0_2 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/52cb/hdl/lib_cdc_v1_0_rfs.vhd" \
-endlib
-makelib ies/interrupt_control_v3_1_4 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/e956/hdl/interrupt_control_v3_1_vh_rfs.vhd" \
-endlib
-makelib ies/axi_gpio_v2_0_15 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/cb07/hdl/axi_gpio_v2_0_vh_rfs.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_axi_gpio_0_0/sim/zynq_design_1_axi_gpio_0_0.vhd" \
-endlib
-makelib ies/blk_mem_gen_v8_3_6 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/4158/simulation/blk_mem_gen_v8_3.v" \
-endlib
-makelib ies/axi_bram_ctrl_v4_0_11 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/9183/hdl/axi_bram_ctrl_v4_0_rfs.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_axi_bram_ctrl_0_0/sim/zynq_design_1_axi_bram_ctrl_0_0.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_0/sim/bd_350b_ila_lib_0.v" \
-endlib
-makelib ies/gigantic_mux \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/c0d4/hdl/gigantic_mux_v1_0_cntr.v" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_1/bd_350b_g_inst_0_gigantic_mux.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_1/sim/bd_350b_g_inst_0.v" \
-endlib
-makelib ies/xlconcat_v2_1_1 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/73b7/hdl/xlconcat_v2_1_vl_rfs.v" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_2/sim/bd_350b_slot_0_aw_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_3/sim/bd_350b_slot_0_w_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_4/sim/bd_350b_slot_0_b_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_5/sim/bd_350b_slot_0_ar_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_6/sim/bd_350b_slot_0_r_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_7/sim/bd_350b_slot_2_aw_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_8/sim/bd_350b_slot_2_w_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_9/sim/bd_350b_slot_2_b_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_10/sim/bd_350b_slot_2_ar_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/ip/ip_11/sim/bd_350b_slot_2_r_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/bd_0/hdl/bd_350b.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_system_ila_0/sim/zynq_design_1_system_ila_0.v" \
-endlib
-makelib ies/proc_sys_reset_v5_0_11 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/5db7/hdl/proc_sys_reset_v5_0_vh_rfs.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_rst_ps7_0_100M_0/sim/zynq_design_1_rst_ps7_0_100M_0.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_axi_bram_ctrl_0_bram_0/sim/zynq_design_1_axi_bram_ctrl_0_bram_0.v" \
-endlib
-makelib ies/generic_baseblocks_v2_1_0 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/f9c1/hdl/generic_baseblocks_v2_1_vl_rfs.v" \
-endlib
-makelib ies/axi_register_slice_v2_1_13 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/55c0/hdl/axi_register_slice_v2_1_vl_rfs.v" \
-endlib
-makelib ies/fifo_generator_v13_1_4 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/ebc2/simulation/fifo_generator_vlog_beh.v" \
-endlib
-makelib ies/fifo_generator_v13_1_4 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/ebc2/hdl/fifo_generator_v13_1_rfs.vhd" \
-endlib
-makelib ies/fifo_generator_v13_1_4 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/ebc2/hdl/fifo_generator_v13_1_rfs.v" \
-endlib
-makelib ies/axi_data_fifo_v2_1_12 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/95b9/hdl/axi_data_fifo_v2_1_vl_rfs.v" \
-endlib
-makelib ies/axi_crossbar_v2_1_14 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/f582/hdl/axi_crossbar_v2_1_vl_rfs.v" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_xbar_0/sim/zynq_design_1_xbar_0.v" \
-endlib
-makelib ies/axi_protocol_converter_v2_1_13 \
  "../../../../embedded_lab_2.srcs/sources_1/bd/zynq_design_1/ipshared/f0ae/hdl/axi_protocol_converter_v2_1_vl_rfs.v" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../bd/zynq_design_1/ip/zynq_design_1_auto_pc_0/sim/zynq_design_1_auto_pc_0.v" \
  "../../../bd/zynq_design_1/ip/zynq_design_1_auto_pc_1/sim/zynq_design_1_auto_pc_1.v" \
  "../../../bd/zynq_design_1/hdl/zynq_design_1.v" \
-endlib
-makelib ies/xil_defaultlib \
  glbl.v
-endlib


text:
Real Function prbbar(srt)
  Parameter (pimass=0.140, arho=0.77, aomega=0.782)
  Common /ppb1/ene, factr2(6), fsum, ppinnb, s, wtot
  Common /ppmm/pprr, ppee, pppe, rpre, xopoe, rree
  Save
  sppb3p = xppbar(srt)*factr2(3)*ene/fsum
  pi2 = (s-(pimass+arho)**2)*(s-(pimass-arho)**2)/4/s
  prbbar = 4./27.*sppb3p/pi2*wtot
  Return
End Function prbbar

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


!=====================================================================*
!                                                                     *
!   Software Name : ppohBEM                                           *
!         Version : 0.1                                               *
!                                                                     *
!   License                                                           *
!     This file is part of ppohBEM.                                   *
!     ppohBEM is a free software, you can use it under the terms      *
!     of The MIT License (MIT). See LICENSE file and User's guide     *
!     for more details.                                               *
!                                                                     *
!   ppOpen-HPC project:                                               *
!     Open Source Infrastructure for Development and Execution of     *
!     Large-Scale Scientific Applications on Post-Peta-Scale          *
!     Supercomputers with Automatic Tuning (AT).                      *
!                                                                     *
!   Organizations:                                                    *
!     The University of Tokyo                                         *
!       - Information Technology Center                               *
!       - Atmosphere and Ocean Research Institute (AORI)              *
!       - Interfaculty Initiative in Information Studies              *
!         /Earthquake Research Institute (ERI)                        *
!       - Graduate School of Frontier Science                         *
!     Kyoto University                                                *
!       - Academic Center for Computing and Media Studies             *
!     Hokkaido University                                             *
!       - Information Initiative Center                               *
!     Japan Agency for Marine-Earth Science and Technology (JAMSTEC)  *
!                                                                     *
!   Sponsorship:                                                      *
!     Japan Science and Technology Agency (JST), Basic Research       *
!     Programs: CREST, Development of System Software Technologies    *
!     for post-Peta Scale High Performance Computing.                 *
!                                                                     *
!   Copyright (c) 2014 <Takeshi Iwashita, Takeshi Mifune, Yuki Noseda,*
!                    Yasuhito Takahashi, Masatoshi Kawai, Akihiro Ida>*
!                                                                     *
!=====================================================================*
module m_ppohBEM_bembb2hacapk
  use m_HACApK_use
  use m_HACApK_hdump
  implicit real*8(a-h,o-z)
contains

!*** bembb2hacapk
! np                ; coordinates of vertexes: x,y,z
! face2node         ; nodes which compose face
! dble_para_fc      ; real*8 parameters sat on each face
! rhs               ; right hand side vector
! sol               ; solution vector
! param             ; parameter for hacapk
! lpmd              ; data for mpi
! int_para_fc       ; integer parameters sat on each face
! nond              ; number of vertexes
! nofc              ; number of faces
! nond_on_face      ; number of vertexes which compose a face
! number_element_dof; degree of freedom on a face
! ndble_para_fc     ; number of real*8 parameters
! nint_para_fc      ; number of integer parameters
  integer function bembb2hacapk(st_bemv, st_ctl, np, face2node,dble_para_fc,rhs,sol,ztol, &
       int_para_fc,nond,nofc,nond_on_face,number_element_dof,ndble_para_fc,nint_para_fc)
    include 'mpif.h'
    type(coordinate) np(nond)
    real*8 :: rhs(nofc*number_element_dof),sol(nofc*number_element_dof),dble_para_fc(ndble_para_fc,nofc),ztol
    integer*4 :: face2node(nond_on_face,nofc),int_para_fc(nint_para_fc,nofc)
    integer*4, dimension(:), allocatable :: lwww
    integer*4,pointer :: lpmd(:)
    type(st_HACApK_leafmtxp) :: st_leafmtxp
    type(st_HACApK_lcontrol) :: st_ctl
    type(st_HACApK_calc_entry) :: st_bemv
    real*8,dimension(:,:),allocatable :: zgmin,zgmax,zgmid
1000 format(5(a,i10)/)
2000 format(5(a,f10.4)/)

    lpmd => st_ctl%lpmd(:)
    mpinr=lpmd(3); mpilog=lpmd(4); nrank=lpmd(2); icomm=lpmd(1)
    if(st_ctl%param(1)>1 .and. mpinr==0)  print*,'func bembb2hacapk start'
    allocate(st_bemv%int_para_fc(nint_para_fc,nofc), st_bemv%face2node(nond_on_face,nofc))
    allocate(st_bemv%dble_para_fc(ndble_para_fc,nofc))
    allocate(st_bemv%np(nond))
    st_bemv%nond=nond; st_bemv%nofc=nofc;
    st_bemv%np=np;
    if(nint_para_fc  /= 0) then
       st_bemv%int_para_fc=int_para_fc
    endif
    st_bemv%nint_para_fc=nint_para_fc
    if(ndble_para_fc /= 0) then
       st_bemv%dble_para_fc=dble_para_fc
    endif
    st_bemv%ndble_para_fc=ndble_para_fc; st_bemv%nond_on_face=nond_on_face;
    st_bemv%number_element_dof=number_element_dof; st_bemv%face2node=face2node

    ndim=3
    allocate(st_bemv%zx(nond),st_bemv%zy(nond),st_bemv%zz(nond),stat=ierr)
    if(ierr/=0)then
       goto 9999
    endif
    do il=1,nond
       st_bemv%zx(il)=np(il)%x; st_bemv%zy(il)=np(il)%y; st_bemv%zz(il)=np(il)%z;
    enddo

    allocate(zgmid(nofc,ndim))
    do il=1,nofc
       n1 = face2node(1,il)+1; n2 = face2node(2,il)+1; n3 = face2node(3,il)+1;
       zgmid(il,1) = (np(n1)%x+np(n2)%x+np(n3)%x)/3.0d0; ! center of balance
       zgmid(il,2) = (np(n1)%y+np(n2)%y+np(n3)%y)/3.0d0
       zgmid(il,3) = (np(n1)%z+np(n2)%z+np(n3)%z)/3.0d0
!***************************************************
!    rhs(il)=zgmid(il,3)
!***************************************************
    enddo
    ! default
    ! bembb2hacapk= hacapk_gensolv(st_leafmtxp,st_bemv,st_ctl,zgmid,rhs,sol,ztol)
    ! dump H matrix
    call hacapk_hdump(st_leafmtxp,st_bemv,st_ctl,zgmid,rhs,sol,ztol)
    lrtrn=HACApK_free_leafmtxp(st_leafmtxp)
    if(st_ctl%param(1)>1 .and. mpinr==0)  print*,'func bembb2hacapk end'
9999 continue
    bembb2hacapk=lrtrn
  endfunction bembb2hacapk

endmodule m_ppohBEM_bembb2hacapk


text:
      program bridge

      parameter(maxmodes=5000)
      implicit real*8 (a-h,o-z)
      character*3 labs(maxmodes)
      dimension freq(maxmodes),aint(maxmodes),anew(maxmodes)
      character*132 card
      read(5,*)natoms
      nmode=(natoms-2)*3
      lp=(nmode-1)/3+1
      now=0
      do 1 i=1,lp
      read(5,'(a132)')card
      read(card,*)(labs(now+j),j=1,3)
      read(5,'(a132)')card
      read(card(17:132),*)(freq(now+j),j=1,3)
      read(5,'(a132)')card
      read(card(17:132),*)(aint(now+j),j=1,3)
      now=now+3
1     continue
      do 2 i=1,nmode
      write(6,3)i,labs(i),freq(i),aint(i)
3      format(i4,a3,2f10.3)
2      continue
      energy=6.0d00
      ec=energy*1.60217653d-12
      call fortran_c_wrapper(freq,aint,ec,anew,nmode,T)
      write(6,10)T
10     format(' temperature is ',f10.3)
      do 12 i=1,nmode
      write(6,15)i,anew(i)
15     format(i5,e15.5)
12     continue
      stop
      end

max_stars_repo_name: PAHdb/CPP-Backend

max_issues_repo_name: PAHdb/CPP-Backend

max_forks_repo_name: PAHdb/CPP-Backend


text:
************************************************************************
*
*     Subroutine TACIT_H                              Called by: DHUMSL
*
*     evaluate a particular parameter set with respect to linear
*     regression and Tobit models
*
C
C      NP     NUMBER OF PARAMETERS IN MODEL
C      Y      VECTOR OF DEPENDENT OBSERVATIONS
C      YD     VECTOR OF CENSORING THRESHOLDS; Y IS OBSERVED IFF Y>YD
C      X      MATRIX OF X.  NOTE THAT X(@,1) IS A VECTOR OF 1'S
C      B_MLE   VECTOR OF PARAMETER ESTIMATES
C      HF      VALUE OF THE L.F. GIVEN B_MLE,X,Y AND YD. (OUTPUT)
*
************************************************************************
      SUBROUTINE TACIT_H(B_MLE,HF,NOBSC,NP,X,Y,YD,CENSFLAG,NOBSCI)
*
*     dimensional parameters
*
      INCLUDE 'fmodules.inc'
*
*     subroutine arguments
*
      INTEGER*4 NOBSC,NP,NOBSCI
      LOGICAL CENSFLAG(*)
      DOUBLE PRECISION B_MLE(*),Y(*),YD(*),X(NOBSCI,*),
     &                 HF(MAXPARMS+1,*)
*
*     local vars
*
      INTEGER*4 I,K,K2
      DOUBLE PRECISION S,S_2,XSI,A,B,GT2,Z
      DOUBLE PRECISION HT(MAXPARMS,MAXPARMS)
*
*     function declaration
*
      DOUBLE PRECISION PRED
*
*
*
      S = EXP(B_MLE(NP+1))
      S_2 = (1.D0/(S*S))

      DO 20 K=1,NP+1
         DO 10 K2=1,NP+1
            HF(K,K2) = 0.D0
 10      CONTINUE
 20   CONTINUE

      DO 50 I=1,NOBSC
         IF (CENSFLAG(I)) THEN
            XSI = (YD(I)-PRED(NP,NOBSCI,I,X,B_MLE))/S
            CALL TACIT_CALC(XSI,A)
            B = -A*(XSI+A)
            GT2 = -XSI*A/S
            HT(1,1) = S_2*B
            HT(1,2) = S_2*(A+XSI*B)
            HT(2,2) = XSI*S_2*(2.D0*A + XSI*B)
         ELSE
            Z = (Y(I)-PRED(NP,NOBSCI,I,X,B_MLE))/S
            GT2 = (Z**2-1.D0)/S
            HT(1,1) = -S_2
            HT(1,2) = -S_2*(2.D0*Z)
            HT(2,2) = -S_2*(3.D0*Z**2-1.D0)
         ENDIF
*
*        construct negatives for IMSL routine DUMIAH (a minimizer)
*
         DO 40 K2=1,NP
            HF(NP+1,K2) = HF(NP+1,K2) - S*HT(1,2)*X(I,K2)
            HF(K2,NP+1) = HF(NP+1,K2)
            DO 30 K=1,NP
               HF(K,K2) = HF(K,K2) - HT(1,1)*X(I,K)*X(I,K2)
 30         CONTINUE
 40      CONTINUE
         HF(NP+1,NP+1) =  HF(NP+1,NP+1) - S*(GT2 + S*HT(2,2))
 50   CONTINUE

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!        DECLARATION OF SHARED VARIABLES (PARAMETERS, VARIABLES, TEMP)       !!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	MODULE PARAMETERS
	IMPLICIT NONE
	INTEGER			:: ne,ni1,ni2,tot_n,total_time,ne_pw,ni_pw1,ni_pw2,n_pw,one_sec,pref_pw
	INTEGER			:: n_trials, n_rows, row_count
	REAL*8			:: pr(300),dt,dpi,seed0,max_step
	CHARACTER*6		:: folder
	END MODULE

	MODULE VARIABLES
	IMPLICIT NONE
	REAL*8,ALLOCATABLE	:: msynw(:),spkt(:),inp_hz(:) !inp_hz is the step rate
	LOGICAL,ALLOCATABLE	:: spk(:)
	LOGICAL			:: spk_post
	REAL*8			:: x(30),tr,spkt_post,inp_fr
	REAL*8,ALLOCATABLE	:: patt_time(:),inp_patt(:),spkr(:), step(:), ex_spk(:), in1_spk(:), in2_spk(:)
	REAL*8,ALLOCATABLE	:: msynw0(:),msynw1(:),msynw2(:),msynw5(:),msynw10(:),msynw20(:)
	INTEGER			:: sims_type
	INTEGER,ALLOCATABLE  		:: pulse_sigs(:), weight_times(:),cond(:), weight_tm(:)
	REAL,ALLOCATABLE	:: phasic(:,:),tonic(:,:),phasic_f(:,:),tonic_f(:,:)
	INTEGER,ALLOCATABLE	:: post_phas(:), post_tonic(:), counter(:), pathwy(:)
	END MODULE

	MODULE TEMP
	IMPLICIT NONE
	REAL*8,ALLOCATABLE	:: tmp_inp(:),tmpi1(:),tmpi2(:)
	END MODULE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


module module_fr_fire_atm

use module_model_constants, only: cp,xlv
use module_fr_fire_util

contains

SUBROUTINE fire_tendency( &
    ids,ide, kds,kde, jds,jde, &
    ims,ime, kms,kme, jms,jme, &
    its,ite, kts,kte, jts,jte, &
    grnhfx,grnqfx,canhfx,canqfx, &
    alfg,alfc,z1can, &
    zs,z_at_w,dz8w,mu,rho, &
    rthfrten,rqvfrten)
   IMPLICIT NONE
   INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
                           ims,ime, kms,kme, jms,jme, &
                           its,ite, kts,kte, jts,jte
   REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx
   REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx
   REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs
   REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu
   REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w
   REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w
   REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho
   REAL, INTENT(in) :: alfg
   REAL, INTENT(in) :: alfc
   REAL, INTENT(in) :: z1can
   REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) :: &
       rthfrten, &
       rqvfrten
   INTEGER :: i,j,k
   INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
   REAL :: cp_i
   REAL :: rho_i
   REAL :: xlv_i
   REAL :: z_w
   REAL :: fact_g, fact_c
   REAL :: alfg_i, alfc_i
   REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
        do j=jts,jte
            do k=kts,min(kte+1,kde)
               do i=its,ite
                   rthfrten(i,k,j)=0.
                   rqvfrten(i,k,j)=0.
               enddo
            enddo
        enddo
   cp_i = 1./cp
   xlv_i = 1./xlv
   alfg_i = 1./alfg
   alfc_i = 1./alfc
   call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
   call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
   i_st = MAX(its,ids+1)
   i_en = MIN(ite,ide-1)
   k_st = kts
   k_en = MIN(kte,kde-1)
   j_st = MAX(jts,jds+1)
   j_en = MIN(jte,jde-1)
   DO j = j_st,j_en
      DO k = k_st,k_en
         DO i = i_st,i_en
            z_w = z_at_w(i,k,j) - zs(i,j)
            fact_g = cp_i * EXP( - alfg_i * z_w )
            IF ( z_w < z1can ) THEN
               fact_c = cp_i
            ELSE
               fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
            END IF
            hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
            fact_g = xlv_i * EXP( - alfg_i * z_w )
            IF (z_w < z1can) THEN
               fact_c = xlv_i
            ELSE
               fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
            END IF
            qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
         END DO
      END DO
   END DO
   DO j = j_st,j_en
      DO k = k_st,k_en-1
         DO i = i_st,i_en
            rho_i = 1./rho(i,k,j)
            rthfrten(i,k,j) = - mu(i,j) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
            rqvfrten(i,k,j) = - mu(i,j) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
         END DO
      END DO
   END DO
   call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
   call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
   RETURN
END SUBROUTINE fire_tendency
end module module_fr_fire_atm


text:
module prec
!
! Define 8-byte size
! Use 6 for real*4 arithmetic.  Reason is to get bfb agreement with
! Karl Taylor's original code.
!
!  integer, parameter :: r8 = selected_real_kind(6)
  integer, parameter :: r8 = selected_real_kind(12)
end module prec

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD