text:
!-------------------------------------------------------------------------------
!>
!! Prognostic variable module
!!
!! @par Description
!!         This module contains the prognostic variables
!!
!! @author H.Tomita
!!
!! @par History
!! @li      2004-02-17 (H.Tomita)   Imported from igdc-4.34
!! @li      2006-05-17 (H.Tomita)   Add optiion 'no_input'.
!! @li      2006-06-02 (H.Tomita)   Avoid the blocking of communication in prgvar_set_in.
!! @li      2007-02-09 (M.Satoh)    Change default: input_direct_access, output_direct_access
!! @li      2007-06-05 (M.Satoh)    add option: TRC_VMAX_INPUT
!! @li      2008-04-12 (T.Mitsui)   add prgvar_get_noq
!! @li      2008-04-12 (T.Mitsui)   add Nc,Nr,Ni,Ns,Ng restart and diagnosis option
!! @li      2008-05-24 (T.Mitsui)   trivial fix (only for log file)
!! @li      2009-04-14 (T.Mitsui)   add option diag_qi (Re-analysis data implicitly contains ice)
!! @li      2009-08-18 (T.Mitusi)   add option opt_qcqi_to_qv (qc and qi derived from Re-analysis may not be reliable)
!! @li      2011-07-22 (T.Ohno)     add subroutines for plane hgrid systems
!! @li      2011-09-03 (H.Yashiro)  New I/O
!! @li      2011-11-30 (S.Iga)      Avoid irregal isend/irecv buffer array (thanks to T.Inoue)
!! @li      2012-07-23 (H.Yashiro)  [add] Water Isotope by K.Yoshimura
!!
!<
module mod_prgvar
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use mod_precision
  use mod_debug
  use mod_adm, only: &
     ADM_LOG_FID,  &
     ADM_MAXFNAME, &
     ADM_NSYS
  use mod_runconf, only: &
     PRG_vmax0,  &
     I_RHOG,     &
     I_RHOGVX,   &
     I_RHOGVY,   &
     I_RHOGVZ,   &
     I_RHOGW,    &
     I_RHOGE,    &
     I_RHOGQstr, &
     I_RHOGQend, &
     DIAG_vmax0, &
     I_pre,      &
     I_tem,      &
     I_vx,       &
     I_vy,       &
     I_vz,       &
     I_w,        &
     I_qstr,     &
     I_qend
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: prgvar_setup
  public :: prgvar_get
  public :: prgvar_get_noq
  public :: prgvar_get_withdiag
  public :: prgvar_set
  public :: prgvar_get_in
  public :: prgvar_get_in_withdiag
  public :: prgvar_set_in

  public :: restart_input
  public :: restart_output

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, allocatable :: PRG_var (:,:,:,:)
  real(RP), public, allocatable :: PRG_var1(:,:,:,:)
  real(RP), public, allocatable :: DIAG_var(:,:,:,:)

  character(len=ADM_MAXFNAME), public :: restart_input_basename  = ''
  character(len=ADM_MAXFNAME), public :: restart_output_basename = ''

  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private, allocatable :: PRG_var_pl (:,:,:,:)
  real(RP), private, allocatable :: PRG_var1_pl(:,:,:,:)
  real(RP), private, allocatable :: DIAG_var_pl(:,:,:,:)

  integer, private :: TRC_vmax_input ! number of input tracer variables

  character(len=ADM_MAXFNAME), private :: layername      = ''
  character(len=ADM_MAXFNAME), private :: input_io_mode  = 'ADVANCED'
  character(len=ADM_MAXFNAME), private :: output_io_mode = 'ADVANCED'
  logical,                     private :: allow_missingq = .false.

  !-----------------------------------------------------------------------------
contains

  !-----------------------------------------------------------------------------
  !>
  !> Description of the subroutine prgvar_setup
  !>
  subroutine prgvar_setup
    use mod_adm, only: &
       ADM_CTL_FID,   &
       ADM_proc_stop, &
       ADM_gall,      &
       ADM_gall_pl,   &
       ADM_kall,      &
       ADM_lall,      &
       ADM_lall_pl
    use mod_runconf, only: &
       PRG_vmax,  &
       DIAG_vmax, &
       TRC_vmax
    implicit none

    character(len=ADM_MAXFNAME) :: input_basename    = ''
    character(len=ADM_MAXFNAME) :: output_basename   = 'restart'
    character(len=ADM_MAXFNAME) :: restart_layername = ''

    namelist / RESTARTPARAM / &
       TRC_vmax_input,    &
       input_basename,    &
       output_basename,   &
       restart_layername, &
       input_io_mode,     &
       output_io_mode,    &
       allow_missingq

    integer :: ierr
    !---------------------------------------------------------------------------

    TRC_vmax_input = TRC_vmax

    !--- read parameters
    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '+++ Module[prgvar]/Category[nhm share]'
    rewind(ADM_CTL_FID)
    read(ADM_CTL_FID,nml=RESTARTPARAM,iostat=ierr)
    if ( ierr < 0 ) then
       write(ADM_LOG_FID,*) '*** RESTARTPARAM is not specified. use default.'
    elseif( ierr > 0 ) then
       write(*,          *) 'xxx Not appropriate names in namelist RESTARTPARAM. STOP.'
       write(ADM_LOG_FID,*) 'xxx Not appropriate names in namelist RESTARTPARAM. STOP.'
       call ADM_proc_stop
    endif
    write(ADM_LOG_FID,nml=RESTARTPARAM)

    restart_input_basename  = input_basename
    restart_output_basename = output_basename
    layername               = restart_layername

    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '*** io_mode for restart, input : ', trim(input_io_mode)
    if    ( input_io_mode == 'ADVANCED'     ) then
    elseif( input_io_mode == 'LEGACY'       ) then
    elseif( input_io_mode == 'IDEAL'        ) then
    elseif( input_io_mode == 'IDEAL_TRACER' ) then
    else
       write(ADM_LOG_FID,*) 'xxx Invalid input_io_mode. STOP.'
       call ADM_proc_stop
    endif

    write(ADM_LOG_FID,*) '*** io_mode for restart, output: ', trim(output_io_mode)
    if    ( output_io_mode == 'ADVANCED' ) then
    elseif( output_io_mode == 'LEGACY'   ) then
    else
       write(ADM_LOG_FID,*) 'xxx Invalid output_io_mode. STOP'
       call ADM_proc_stop
    endif

    if ( allow_missingq ) then
       write(ADM_LOG_FID,*) '*** Allow missing tracer in restart file.'
       write(ADM_LOG_FID,*) '*** Value will set to zero for missing tracer.'
    endif

    allocate( PRG_var    (ADM_gall,   ADM_kall,ADM_lall,   PRG_vmax) )
    allocate( PRG_var_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,PRG_vmax) )

    allocate( PRG_var1   (ADM_gall,   ADM_kall,ADM_lall,   PRG_vmax) )
    allocate( PRG_var1_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl,PRG_vmax) )

    allocate( DIAG_var   (ADM_gall,   ADM_kall,ADM_lall,   DIAG_vmax) )
    allocate( DIAG_var_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl,DIAG_vmax) )

    return
  end subroutine prgvar_setup

  !-----------------------------------------------------------------------------
  !>
  !> get prognostic variables from prg[num].
  !>
  subroutine prgvar_get( &
       rhog,   rhog_pl,   &
       rhogvx, rhogvx_pl, &
       rhogvy, rhogvy_pl, &
       rhogvz, rhogvz_pl, &
       rhogw,  rhogw_pl,  &
       rhoge,  rhoge_pl,  &
       rhogq,  rhogq_pl,  &
       num                )
    use mod_adm, only: &
       ADM_have_pl, &
       ADM_gall,    &
       ADM_gall_pl, &
       ADM_kall,    &
       ADM_lall,    &
       ADM_lall_pl
    use mod_runconf, only: &
       TRC_vmax
    implicit none

    real(RP), intent(out) :: rhog     (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhog_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvx   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvy   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvz   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogw    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhoge    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhoge_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogq    (ADM_gall,   ADM_kall,ADM_lall   ,TRC_vmax)
    real(RP), intent(out) :: rhogq_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,TRC_vmax)
    integer, intent(in)  :: num

    integer :: n, k, l, nq
    !---------------------------------------------------------------------------

    if ( num == 0 ) then

       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do n = 1, ADM_gall
          rhog  (n,k,l) = PRG_var(n,k,l,I_RHOG  )
          rhogvx(n,k,l) = PRG_var(n,k,l,I_RHOGVX)
          rhogvy(n,k,l) = PRG_var(n,k,l,I_RHOGVY)
          rhogvz(n,k,l) = PRG_var(n,k,l,I_RHOGVZ)
          rhogw (n,k,l) = PRG_var(n,k,l,I_RHOGW )
          rhoge (n,k,l) = PRG_var(n,k,l,I_RHOGE )
       enddo
       enddo
       enddo

       do nq = 1, TRC_vmax
       do l  = 1, ADM_lall
       do k  = 1, ADM_kall
       do n  = 1, ADM_gall
          rhogq(n,k,l,nq) = PRG_var(n,k,l,PRG_vmax0+nq)
       enddo
       enddo
       enddo
       enddo

       if ( ADM_have_pl ) then

          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do n = 1, ADM_gall_pl
             rhog_pl  (n,k,l) = PRG_var_pl(n,k,l,I_RHOG  )
             rhogvx_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVX)
             rhogvy_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVY)
             rhogvz_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVZ)
             rhogw_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGW )
             rhoge_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGE )
          enddo
          enddo
          enddo

          do nq = 1, TRC_vmax
          do l  = 1, ADM_lall_pl
          do k  = 1, ADM_kall
          do n  = 1, ADM_gall_pl
             rhogq_pl(n,k,l,nq) = PRG_var_pl(n,k,l,PRG_vmax0+nq)
          enddo
          enddo
          enddo
          enddo

       endif

    elseif( num == 1 ) then

       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do n = 1, ADM_gall
          rhog  (n,k,l) = PRG_var1(n,k,l,I_RHOG  )
          rhogvx(n,k,l) = PRG_var1(n,k,l,I_RHOGVX)
          rhogvy(n,k,l) = PRG_var1(n,k,l,I_RHOGVY)
          rhogvz(n,k,l) = PRG_var1(n,k,l,I_RHOGVZ)
          rhogw (n,k,l) = PRG_var1(n,k,l,I_RHOGW )
          rhoge (n,k,l) = PRG_var1(n,k,l,I_RHOGE )
       enddo
       enddo
       enddo

       do nq = 1, TRC_vmax
       do l  = 1, ADM_lall
       do k  = 1, ADM_kall
       do n  = 1, ADM_gall
          rhogq(n,k,l,nq) = PRG_var1(n,k,l,PRG_vmax0+nq)
       enddo
       enddo
       enddo
       enddo

       if ( ADM_have_pl ) then

          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do n = 1, ADM_gall_pl
             rhog_pl  (n,k,l) = PRG_var1_pl(n,k,l,I_RHOG  )
             rhogvx_pl(n,k,l) = PRG_var1_pl(n,k,l,I_RHOGVX)
             rhogvy_pl(n,k,l) = PRG_var1_pl(n,k,l,I_RHOGVY)
             rhogvz_pl(n,k,l) = PRG_var1_pl(n,k,l,I_RHOGVZ)
             rhogw_pl (n,k,l) = PRG_var1_pl(n,k,l,I_RHOGW )
             rhoge_pl (n,k,l) = PRG_var1_pl(n,k,l,I_RHOGE )
          enddo
          enddo
          enddo

          do nq = 1, TRC_vmax
          do l  = 1, ADM_lall_pl
          do k  = 1, ADM_kall
          do n  = 1, ADM_gall_pl
             rhogq_pl(n,k,l,nq) = PRG_var1_pl(n,k,l,PRG_vmax0+nq)
          enddo
          enddo
          enddo
          enddo

       endif

    endif

    return
  end subroutine prgvar_get

  !----------------------------------------------------------------------------
  subroutine prgvar_get_noq( &
       rhog,   rhog_pl,   &
       rhogvx, rhogvx_pl, &
       rhogvy, rhogvy_pl, &
       rhogvz, rhogvz_pl, &
       rhogw,  rhogw_pl,  &
       rhoge,  rhoge_pl   )
    use mod_adm, only: &
       ADM_have_pl, &
       ADM_gall,    &
       ADM_gall_pl, &
       ADM_kall,    &
       ADM_lall,    &
       ADM_lall_pl
    implicit none

    real(RP), intent(out) :: rhog     (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhog_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvx   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvy   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvz   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogw    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhoge    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhoge_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)

    integer :: n, k, l
    !---------------------------------------------------------------------------

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do n = 1, ADM_gall
       rhog  (n,k,l) = PRG_var(n,k,l,I_RHOG  )
       rhogvx(n,k,l) = PRG_var(n,k,l,I_RHOGVX)
       rhogvy(n,k,l) = PRG_var(n,k,l,I_RHOGVY)
       rhogvz(n,k,l) = PRG_var(n,k,l,I_RHOGVZ)
       rhogw (n,k,l) = PRG_var(n,k,l,I_RHOGW )
       rhoge (n,k,l) = PRG_var(n,k,l,I_RHOGE )
    enddo
    enddo
    enddo

    if ( ADM_have_pl ) then

       do l = 1, ADM_lall_pl
       do k = 1, ADM_kall
       do n = 1, ADM_gall_pl
          rhog_pl  (n,k,l) = PRG_var_pl(n,k,l,I_RHOG  )
          rhogvx_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVX)
          rhogvy_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVY)
          rhogvz_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVZ)
          rhogw_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGW )
          rhoge_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGE )
       enddo
       enddo
       enddo

    endif

    return
  end subroutine prgvar_get_noq

  !-----------------------------------------------------------------------------
  !>
  !> get prognostic variables from prg[num].
  !>
  subroutine prgvar_get_withdiag( &
       rhog,   rhog_pl,   &
       rhogvx, rhogvx_pl, &
       rhogvy, rhogvy_pl, &
       rhogvz, rhogvz_pl, &
       rhogw,  rhogw_pl,  &
       rhoge,  rhoge_pl,  &
       rhogq,  rhogq_pl,  &
       rho,    rho_pl,    &
       pre,    pre_pl,    &
       tem,    tem_pl,    &
       vx,     vx_pl,     &
       vy,     vy_pl,     &
       vz,     vz_pl,     &
       w,      w_pl,      &
       q,      q_pl       )
    use mod_adm, only: &
       ADM_have_pl, &
       ADM_gall,    &
       ADM_gall_pl, &
       ADM_kall,    &
       ADM_lall,    &
       ADM_lall_pl
    use mod_vmtr, only: &
       VMTR_RGSGAM2,    &
       VMTR_RGSGAM2_pl
    use mod_runconf, only: &
       TRC_vmax
    use mod_cnvvar, only: &
       cnvvar_prg2diag
    implicit none

    real(RP), intent(out) :: rhog     (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhog_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvx   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvy   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogvz   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogvz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogw    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhoge    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rhoge_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: rhogq    (ADM_gall,   ADM_kall,ADM_lall   ,TRC_vmax)
    real(RP), intent(out) :: rhogq_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,TRC_vmax)
    real(RP), intent(out) :: rho      (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: rho_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: pre      (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: pre_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: tem      (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: tem_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: vx       (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: vx_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: vy       (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: vy_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: vz       (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: vz_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: w        (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(out) :: w_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: q        (ADM_gall,   ADM_kall,ADM_lall,   TRC_vmax)
    real(RP), intent(out) :: q_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl,TRC_vmax)

    integer :: n, k, l, nq
    !---------------------------------------------------------------------------

    call cnvvar_prg2diag( PRG_var (:,:,:,:), PRG_var_pl (:,:,:,:), & !--- [IN]
                          DIAG_var(:,:,:,:), DIAG_var_pl(:,:,:,:)  ) !--- [OUT]

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do n = 1, ADM_gall
       rhog  (n,k,l) = PRG_var(n,k,l,I_RHOG  )
       rhogvx(n,k,l) = PRG_var(n,k,l,I_RHOGVX)
       rhogvy(n,k,l) = PRG_var(n,k,l,I_RHOGVY)
       rhogvz(n,k,l) = PRG_var(n,k,l,I_RHOGVZ)
       rhogw (n,k,l) = PRG_var(n,k,l,I_RHOGW )
       rhoge (n,k,l) = PRG_var(n,k,l,I_RHOGE )

       rho   (n,k,l) = PRG_var(n,k,l,I_RHOG) * VMTR_RGSGAM2(n,k,l)

       pre   (n,k,l) = DIAG_var(n,k,l,I_pre)
       tem   (n,k,l) = DIAG_var(n,k,l,I_tem)
       vx    (n,k,l) = DIAG_var(n,k,l,I_vx )
       vy    (n,k,l) = DIAG_var(n,k,l,I_vy )
       vz    (n,k,l) = DIAG_var(n,k,l,I_vz )
       w     (n,k,l) = DIAG_var(n,k,l,I_w  )
    enddo
    enddo
    enddo

    do nq = 1, TRC_vmax
    do l  = 1, ADM_lall
    do k  = 1, ADM_kall
    do n  = 1, ADM_gall
       rhogq(n,k,l,nq) = PRG_var(n,k,l,PRG_vmax0+nq)

       q    (n,k,l,nq) = DIAG_var(n,k,l,DIAG_vmax0+nq)
    enddo
    enddo
    enddo
    enddo

    if ( ADM_have_pl ) then

       do l = 1, ADM_lall_pl
       do k = 1, ADM_kall
       do n = 1, ADM_gall_pl
          rhog_pl  (n,k,l) = PRG_var_pl(n,k,l,I_RHOG  )
          rhogvx_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVX)
          rhogvy_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVY)
          rhogvz_pl(n,k,l) = PRG_var_pl(n,k,l,I_RHOGVZ)
          rhogw_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGW )
          rhoge_pl (n,k,l) = PRG_var_pl(n,k,l,I_RHOGE )

          rho_pl   (n,k,l) = PRG_var_pl(n,k,l,I_RHOG) * VMTR_RGSGAM2_pl(n,k,l)

          pre_pl   (n,k,l) = DIAG_var_pl(n,k,l,I_pre)
          tem_pl   (n,k,l) = DIAG_var_pl(n,k,l,I_tem)
          vx_pl    (n,k,l) = DIAG_var_pl(n,k,l,I_vx )
          vy_pl    (n,k,l) = DIAG_var_pl(n,k,l,I_vy )
          vz_pl    (n,k,l) = DIAG_var_pl(n,k,l,I_vz )
          w_pl     (n,k,l) = DIAG_var_pl(n,k,l,I_w  )
       enddo
       enddo
       enddo

       do nq = 1, TRC_vmax
       do l  = 1, ADM_lall_pl
       do k  = 1, ADM_kall
       do n  = 1, ADM_gall_pl
          rhogq_pl(n,k,l,nq) = PRG_var_pl(n,k,l,PRG_vmax0+nq)

          q_pl    (n,k,l,nq) = DIAG_var_pl(n,k,l,DIAG_vmax0+nq)
       enddo
       enddo
       enddo
       enddo

    endif

    return
  end subroutine prgvar_get_withdiag

  !-----------------------------------------------------------------------------
  !>
  !> set prognostic variables to prg[num]- and COMMUNICATION.
  !>
  subroutine prgvar_set( &
       rhog,   rhog_pl,   &
       rhogvx, rhogvx_pl, &
       rhogvy, rhogvy_pl, &
       rhogvz, rhogvz_pl, &
       rhogw,  rhogw_pl,  &
       rhoge,  rhoge_pl,  &
       rhogq,  rhogq_pl,  &
       num                )
    use mod_adm, only: &
       ADM_have_pl, &
       ADM_gall,    &
       ADM_gall_pl, &
       ADM_kall,    &
       ADM_lall,    &
       ADM_lall_pl, &
       ADM_gall_1d, &
       ADM_gmin,    &
       ADM_gmax
    use mod_comm, only: &
       COMM_data_transfer
    use mod_runconf, only: &
       TRC_vmax
    implicit none

    real(RP), intent(in)  :: rhog     (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhog_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvx   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhogvx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvy   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhogvy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvz   (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhogvz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogw    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhoge    (ADM_gall,   ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: rhoge_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogq    (ADM_gall,   ADM_kall,ADM_lall   ,TRC_vmax)
    real(RP), intent(in)  :: rhogq_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,TRC_vmax)
    integer, intent(in)  :: num

    integer :: i, j, suf
    suf(i,j) = ADM_gall_1d * ((j)-1) + (i)

    integer :: n, k, l, nq
    !---------------------------------------------------------------------------

    if ( num == 0 ) then

       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do n = 1, ADM_gall
          PRG_var(n,k,l,I_RHOG  ) = rhog  (n,k,l)
          PRG_var(n,k,l,I_RHOGVX) = rhogvx(n,k,l)
          PRG_var(n,k,l,I_RHOGVY) = rhogvy(n,k,l)
          PRG_var(n,k,l,I_RHOGVZ) = rhogvz(n,k,l)
          PRG_var(n,k,l,I_RHOGW ) = rhogw (n,k,l)
          PRG_var(n,k,l,I_RHOGE ) = rhoge (n,k,l)
       enddo
       enddo
       enddo

       do nq = 1, TRC_vmax
       do l  = 1, ADM_lall
       do k  = 1, ADM_kall
       do n  = 1, ADM_gall
          PRG_var(n,k,l,PRG_vmax0+nq) = rhogq(n,k,l,nq)
       enddo
       enddo
       enddo
       enddo

       if ( ADM_have_pl ) then

          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do n = 1, ADM_gall_pl
             PRG_var_pl(n,k,l,I_RHOG  ) = rhog_pl  (n,k,l)
             PRG_var_pl(n,k,l,I_RHOGVX) = rhogvx_pl(n,k,l)
             PRG_var_pl(n,k,l,I_RHOGVY) = rhogvy_pl(n,k,l)
             PRG_var_pl(n,k,l,I_RHOGVZ) = rhogvz_pl(n,k,l)
             PRG_var_pl(n,k,l,I_RHOGW ) = rhogw_pl (n,k,l)
             PRG_var_pl(n,k,l,I_RHOGE ) = rhoge_pl (n,k,l)
          enddo
          enddo
          enddo

          do nq = 1, TRC_vmax
          do l  = 1, ADM_lall_pl
          do k  = 1, ADM_kall
          do n  = 1, ADM_gall_pl
             PRG_var_pl(n,k,l,PRG_vmax0+nq) = rhogq_pl(n,k,l,nq)
          enddo
          enddo
          enddo
          enddo

       endif

       ! communication
       call COMM_data_transfer( PRG_var, PRG_var_pl )

       PRG_var(suf(ADM_gmax+1,ADM_gmin-1),:,:,:) = PRG_var(suf(ADM_gmax+1,ADM_gmin),:,:,:)
       PRG_var(suf(ADM_gmin-1,ADM_gmax+1),:,:,:) = PRG_var(suf(ADM_gmin,ADM_gmax+1),:,:,:)

    elseif( num == 1 ) then

       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do n = 1, ADM_gall
          PRG_var1(n,k,l,I_RHOG  ) = rhog  (n,k,l)
          PRG_var1(n,k,l,I_RHOGVX) = rhogvx(n,k,l)
          PRG_var1(n,k,l,I_RHOGVY) = rhogvy(n,k,l)
          PRG_var1(n,k,l,I_RHOGVZ) = rhogvz(n,k,l)
          PRG_var1(n,k,l,I_RHOGW ) = rhogw (n,k,l)
          PRG_var1(n,k,l,I_RHOGE ) = rhoge (n,k,l)
       enddo
       enddo
       enddo

       do nq = 1, TRC_vmax
       do l  = 1, ADM_lall
       do k  = 1, ADM_kall
       do n  = 1, ADM_gall
          PRG_var1(n,k,l,PRG_vmax0+nq) = rhogq(n,k,l,nq)
       enddo
       enddo
       enddo
       enddo

       if ( ADM_have_pl ) then

          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do n = 1, ADM_gall_pl
             PRG_var1_pl(n,k,l,I_RHOG  ) = rhog_pl  (n,k,l)
             PRG_var1_pl(n,k,l,I_RHOGVX) = rhogvx_pl(n,k,l)
             PRG_var1_pl(n,k,l,I_RHOGVY) = rhogvy_pl(n,k,l)
             PRG_var1_pl(n,k,l,I_RHOGVZ) = rhogvz_pl(n,k,l)
             PRG_var1_pl(n,k,l,I_RHOGW ) = rhogw_pl (n,k,l)
             PRG_var1_pl(n,k,l,I_RHOGE ) = rhoge_pl (n,k,l)
          enddo
          enddo
          enddo

          do nq = 1, TRC_vmax
          do l  = 1, ADM_lall_pl
          do k  = 1, ADM_kall
          do n  = 1, ADM_gall_pl
             PRG_var1_pl(n,k,l,PRG_vmax0+nq) = rhogq_pl(n,k,l,nq)
          enddo
          enddo
          enddo
          enddo

       endif

       ! communication
       call COMM_data_transfer( PRG_var1, PRG_var1_pl )

       PRG_var1(suf(ADM_gmax+1,ADM_gmin-1),:,:,:) = PRG_var1(suf(ADM_gmax+1,ADM_gmin),:,:,:)
       PRG_var1(suf(ADM_gmin-1,ADM_gmax+1),:,:,:) = PRG_var1(suf(ADM_gmin,ADM_gmax+1),:,:,:)

    endif

    return
  end subroutine prgvar_set

  !-----------------------------------------------------------------------------
  !>
  !> Description of the subroutine prgvar_get_in
  !>
  subroutine prgvar_get_in( &
       rhog,   &
       rhogvx, &
       rhogvy, &
       rhogvz, &
       rhogw,  &
       rhoge,  &
       rhogq   )
    use mod_adm, only: &
       ADM_kall,        &
       ADM_lall,        &
       ADM_IopJop_nmax, &
       ADM_IopJop,      &
       ADM_GIoJo
    use mod_runconf, only: &
       TRC_vmax
    implicit none

    real(RP), intent(out) :: rhog  (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvx(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvy(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvz(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogw (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhoge (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogq (ADM_IopJop_nmax,ADM_kall,ADM_lall,TRC_vmax)

    integer :: n, k, l, nq, ij
    !---------------------------------------------------------------------------

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do n = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       rhog  (n,k,l) = PRG_var(ij,k,l,I_RHOG  )
       rhogvx(n,k,l) = PRG_var(ij,k,l,I_RHOGVX)
       rhogvy(n,k,l) = PRG_var(ij,k,l,I_RHOGVY)
       rhogvz(n,k,l) = PRG_var(ij,k,l,I_RHOGVZ)
       rhogw (n,k,l) = PRG_var(ij,k,l,I_RHOGW )
       rhoge (n,k,l) = PRG_var(ij,k,l,I_RHOGE )
    enddo
    enddo
    enddo

    do nq = 1, TRC_vmax
    do l  = 1, ADM_lall
    do k  = 1, ADM_kall
    do n  = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       rhogq(n,k,l,nq) = PRG_var(ij,k,l,PRG_vmax0+nq)
    enddo
    enddo
    enddo
    enddo

    return
  end subroutine prgvar_get_in

  !-----------------------------------------------------------------------------
  !>
  !> Description of the subroutine prgvar_get_in
  !>
  subroutine prgvar_get_in_withdiag( &
       rhog,   &
       rhogvx, &
       rhogvy, &
       rhogvz, &
       rhogw,  &
       rhoge,  &
       rhogq,  &
       rho,    &
       pre,    &
       tem,    &
       vx,     &
       vy,     &
       vz,     &
       w,      &
       q       )
    use mod_adm, only: &
       ADM_kall,        &
       ADM_lall,        &
       ADM_IopJop_nmax, &
       ADM_IopJop,      &
       ADM_GIoJo
    use mod_vmtr, only: &
       VMTR_RGSGAM2
    use mod_runconf, only: &
       TRC_vmax
    use mod_cnvvar, only: &
       cnvvar_prg2diag
    implicit none

    real(RP), intent(out) :: rhog  (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvx(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvy(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogvz(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogw (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhoge (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: rhogq (ADM_IopJop_nmax,ADM_kall,ADM_lall,TRC_vmax)
    real(RP), intent(out) :: rho   (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: pre   (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: tem   (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: vx    (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: vy    (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: vz    (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: w     (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(out) :: q     (ADM_IopJop_nmax,ADM_kall,ADM_lall,TRC_vmax)

    integer :: n, k, l, nq, ij
    !---------------------------------------------------------------------------

    call cnvvar_prg2diag( PRG_var (:,:,:,:), PRG_var_pl (:,:,:,:), & !--- [IN]
                          DIAG_var(:,:,:,:), DIAG_var_pl(:,:,:,:)  ) !--- [OUT]

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do n = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       rhog  (n,k,l) = PRG_var(ij,k,l,I_RHOG  )
       rhogvx(n,k,l) = PRG_var(ij,k,l,I_RHOGVX)
       rhogvy(n,k,l) = PRG_var(ij,k,l,I_RHOGVY)
       rhogvz(n,k,l) = PRG_var(ij,k,l,I_RHOGVZ)
       rhogw (n,k,l) = PRG_var(ij,k,l,I_RHOGW )
       rhoge (n,k,l) = PRG_var(ij,k,l,I_RHOGE )

       rho   (n,k,l) = PRG_var(ij,k,l,I_RHOG) * VMTR_RGSGAM2(ij,k,l)

       pre   (n,k,l) = DIAG_var(ij,k,l,I_pre)
       tem   (n,k,l) = DIAG_var(ij,k,l,I_tem)
       vx    (n,k,l) = DIAG_var(ij,k,l,I_vx )
       vy    (n,k,l) = DIAG_var(ij,k,l,I_vy )
       vz    (n,k,l) = DIAG_var(ij,k,l,I_vz )
       w     (n,k,l) = DIAG_var(ij,k,l,I_w  )
    enddo
    enddo
    enddo

    do nq = 1, TRC_vmax
    do l  = 1, ADM_lall
    do k  = 1, ADM_kall
    do n  = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       rhogq(n,k,l,nq) = PRG_var(ij,k,l,PRG_vmax0+nq)

       q    (n,k,l,nq) = DIAG_var(ij,k,l,DIAG_vmax0+nq)
    enddo
    enddo
    enddo
    enddo

    return
  end subroutine prgvar_get_in_withdiag

  !-----------------------------------------------------------------------------
  !>
  !> set prognostic variables to prg[num] and COMMUNICATION.
  !>
  subroutine prgvar_set_in( &
       rhog,   &
       rhogvx, &
       rhogvy, &
       rhogvz, &
       rhogw,  &
       rhoge,  &
       rhogq   )
    use mod_adm, only : &
       ADM_kall,        &
       ADM_lall,        &
       ADM_IopJop_nmax, &
       ADM_IopJop,      &
       ADM_GIoJo,       &
       ADM_gall_1d,     &
       ADM_gmax,        &
       ADM_gmin
    use mod_comm, only: &
       COMM_var
    use mod_runconf, only: &
       PRG_vmax, &
       TRC_vmax
    implicit none

    real(RP), intent(in) :: rhog  (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhogvx(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhogvy(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhogvz(ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhogw (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhoge (ADM_IopJop_nmax,ADM_kall,ADM_lall)
    real(RP), intent(in) :: rhogq (ADM_IopJop_nmax,ADM_kall,ADM_lall,TRC_vmax)

    integer :: n, k, l, nq, ij

    integer :: i,j,suf
    suf(i,j) = ADM_gall_1d * ((j)-1) + (i)
    !---------------------------------------------------------------------------

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do n = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       PRG_var(ij,k,l,I_RHOG  ) = rhog  (n,k,l)
       PRG_var(ij,k,l,I_RHOGVX) = rhogvx(n,k,l)
       PRG_var(ij,k,l,I_RHOGVY) = rhogvy(n,k,l)
       PRG_var(ij,k,l,I_RHOGVZ) = rhogvz(n,k,l)
       PRG_var(ij,k,l,I_RHOGW ) = rhogw (n,k,l)
       PRG_var(ij,k,l,I_RHOGE ) = rhoge (n,k,l)
    enddo
    enddo
    enddo

    do nq = 1, TRC_vmax
    do l  = 1, ADM_lall
    do k  = 1, ADM_kall
    do n = 1, ADM_IopJop_nmax
       ij = ADM_IopJop(n,ADM_GIoJo)

       PRG_var(ij,k,l,PRG_vmax0+nq) = rhogq(n,k,l,nq)
    enddo
    enddo
    enddo
    enddo

    ! communication
    call COMM_var( PRG_var, PRG_var_pl, ADM_kall, PRG_vmax )

    PRG_var(suf(ADM_gmax+1,ADM_gmin-1),:,:,:) = PRG_var(suf(ADM_gmax+1,ADM_gmin),:,:,:)
    PRG_var(suf(ADM_gmin-1,ADM_gmax+1),:,:,:) = PRG_var(suf(ADM_gmin,ADM_gmax+1),:,:,:)

    return
  end subroutine prgvar_set_in

  !-----------------------------------------------------------------------------
  subroutine restart_input( basename )
    use mod_misc, only: &
       MISC_get_available_fid, &
       MISC_make_idstr
    use mod_adm, only: &
       ADM_prc_tab, &
       ADM_prc_me,  &
       ADM_gall,    &
       ADM_lall,    &
       ADM_kall,    &
       ADM_kmax,    &
       ADM_kmin
    use mod_fio, only: &
       FIO_input
    use mod_comm, only: &
       COMM_var
    use mod_gtl, only: &
       GTL_max, &
       GTL_min
    use mod_runconf, only: &
       PRG_name,  &
       DIAG_vmax, &
       DIAG_name, &
       TRC_vmax,  &
       TRC_name
    use mod_cnvvar, only: &
       cnvvar_diag2prg
    use mod_ideal_init, only :  & ! [add] R.Yoshida 20121019
       dycore_input, &
       tracer_input
    implicit none

    character(len=ADM_MAXFNAME), intent(in) :: basename

    character(len=ADM_MAXFNAME) :: fname

    real(RP) :: val_max, val_min
    logical :: nonzero

    integer :: fid
    integer :: l, rgnid, nq
    !---------------------------------------------------------------------------

    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '*** read restart/initial data'

    if ( input_io_mode == 'ADVANCED' ) then

       do nq = 1, DIAG_vmax0
          call FIO_input( DIAG_var(:,:,:,nq),basename,DIAG_name(nq), &
                          layername,1,ADM_kall,1                     )
       enddo

       do nq = 1, TRC_vmax_input
          call FIO_input( DIAG_var(:,:,:,DIAG_vmax0+nq),basename,TRC_name(nq), &
                          layername,1,ADM_kall,1,                              &
                          allow_missingq=allow_missingq                        )
       enddo

    elseif( input_io_mode == 'LEGACY' ) then

       do l = 1, ADM_lall
          rgnid = ADM_prc_tab(l,ADM_prc_me)
          call MISC_make_idstr(fname,trim(basename),'rgn',rgnid)
          fid = MISC_get_available_fid()
          open( unit   = fid,                 &
                file   = trim(fname),         &
                form   = 'unformatted',       &
                access = 'direct',            &
                recl   = ADM_gall*ADM_kall*8, &
                status = 'old'                )

          do nq = 1, DIAG_vmax0+TRC_vmax_input
             read(fid,rec=nq) DIAG_var(:,:,l,nq)
          enddo

          close(fid)
       enddo

    elseif( input_io_mode == 'IDEAL' ) then

       write(ADM_LOG_FID,*) '*** make ideal initials'

       call dycore_input( DIAG_var(:,:,:,:) ) !--- [OUT]

    elseif( input_io_mode == 'IDEAL_TRACER' ) then

       do nq = 1, DIAG_vmax0
          call FIO_input( DIAG_var(:,:,:,nq),basename,DIAG_name(nq), &
                          layername,1,ADM_kall,1                     )
       enddo

       write(ADM_LOG_FID,*) '*** make ideal initials for tracer'

       call tracer_input( DIAG_var(:,:,:,DIAG_vmax0+1:DIAG_vmax0+TRC_vmax) ) !--- [OUT]

    endif !--- io_mode

    call COMM_var( DIAG_var, DIAG_var_pl, ADM_kall, DIAG_vmax )

    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '====== data range check : diagnostic variables ======'
    do nq = 1, DIAG_vmax0
       val_max = GTL_max( DIAG_var(:,:,:,nq), DIAG_var_pl(:,:,:,nq), ADM_kall, ADM_kmin, ADM_kmax )
       val_min = GTL_min( DIAG_var(:,:,:,nq), DIAG_var_pl(:,:,:,nq), ADM_kall, ADM_kmin, ADM_kmax )
       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- ', DIAG_name(nq), ': max=', val_max, ', min=', val_min
    enddo

    do nq = 1, TRC_vmax
       val_max = GTL_max( DIAG_var   (:,:,:,DIAG_vmax0+nq), &
                          DIAG_var_pl(:,:,:,DIAG_vmax0+nq), &
                          ADM_kall, ADM_kmin, ADM_kmax      )

       if ( val_max <= 0.0_RP ) then
          nonzero = .false.
       else
          nonzero = .true.
       endif

       val_min = GTL_min( DIAG_var   (:,:,:,DIAG_vmax0+nq), &
                          DIAG_var_pl(:,:,:,DIAG_vmax0+nq), &
                          ADM_kall, ADM_kmin, ADM_kmax, nonzero)

       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- ', TRC_name(nq),  ': max=', val_max, ', min=', val_min
    enddo

    call cnvvar_diag2prg( PRG_var (:,:,:,:), PRG_var_pl (:,:,:,:), & !--- [OUT]
                          DIAG_var(:,:,:,:), DIAG_var_pl(:,:,:,:)  ) !--- [IN]

    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '====== data range check : prognostic variables ======'
    do nq = 1, PRG_vmax0
       val_max = GTL_max( PRG_var(:,:,:,nq), PRG_var_pl(:,:,:,nq), ADM_kall, ADM_kmin, ADM_kmax )
       val_min = GTL_min( PRG_var(:,:,:,nq), PRG_var_pl(:,:,:,nq), ADM_kall, ADM_kmin, ADM_kmax )
       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- ', PRG_name(nq), ': max=', val_max, ', min=', val_min
    enddo

    do nq = 1, TRC_vmax
       val_max = GTL_max( PRG_var   (:,:,:,PRG_vmax0+nq), &
                          PRG_var_pl(:,:,:,PRG_vmax0+nq), &
                          ADM_kall, ADM_kmin, ADM_kmax    )

       if ( val_max <= 0.0_RP ) then
          nonzero = .false.
       else
          nonzero = .true.
       endif

       val_min = GTL_min( PRG_var   (:,:,:,PRG_vmax0+nq),      &
                          PRG_var_pl(:,:,:,PRG_vmax0+nq),      &
                          ADM_kall, ADM_kmin, ADM_kmax, nonzero)

       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- rhog * ', TRC_name(nq),  ': max=', val_max, ', min=', val_min
    enddo

    return
  end subroutine restart_input

  !-----------------------------------------------------------------------------
  subroutine restart_output( basename )
    use mod_misc, only: &
       MISC_get_available_fid, &
       MISC_make_idstr
    use mod_adm, only: &
       ADM_prc_tab, &
       ADM_prc_me,  &
       ADM_gall,    &
       ADM_lall,    &
       ADM_kall,    &
       ADM_kmax,    &
       ADM_kmin
    use mod_fio, only : & ! [add] H.Yashiro 20110819
       FIO_output, &
       FIO_HSHORT, &
       FIO_HMID,   &
       FIO_REAL8
    use mod_time, only : &
       TIME_CTIME
    use mod_gtl, only: &
       GTL_max, &
       GTL_min
    use mod_runconf, only: &
       DIAG_vmax, &
       DIAG_name, &
       TRC_vmax,  &
       TRC_name,  &
       WLABEL
    use mod_cnvvar, only: &
       cnvvar_prg2diag
    implicit none

    character(len=ADM_MAXFNAME), intent(in) :: basename

    character(len=FIO_HMID)   :: desc = 'INITIAL/RESTART_data_of_prognostic_variables' ! [add] H.Yashiro 20110819

    character(len=FIO_HSHORT) :: DLABEL(DIAG_vmax0)
    data DLABEL / 'Pressure ',        &
                  'Temperature ',     &
                  'H-Velocity(XDIR)', &
                  'H-Velocity(YDIR)', &
                  'H-Velocity(ZDIR)', &
                  'V-Velocity '       /

    character(len=FIO_HSHORT) :: DUNIT(DIAG_vmax0)
    data DUNIT /  'Pa',  &
                  'K',   &
                  'm/s', &
                  'm/s', &
                  'm/s', &
                  'm/s'  /

    character(len=FIO_HSHORT) :: WUNIT = 'kg/kg'

    character(len=ADM_MAXFNAME) :: fname

    real(RP) :: val_max, val_min
    logical :: nonzero

    integer :: fid
    integer :: l, rgnid, nq
    !---------------------------------------------------------------------------

    call cnvvar_prg2diag( PRG_var (:,:,:,:), PRG_var_pl (:,:,:,:), & !--- [IN]
                          DIAG_var(:,:,:,:), DIAG_var_pl(:,:,:,:)  ) !--- [OUT]

    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '====== data range check : prognostic variables ======'
    do nq = 1, DIAG_vmax0
       val_max = GTL_max( DIAG_var   (:,:,:,nq),       &
                          DIAG_var_pl(:,:,:,nq),       &
                          ADM_kall, ADM_kmin, ADM_kmax )
       val_min = GTL_min( DIAG_var   (:,:,:,nq),       &
                          DIAG_var_pl(:,:,:,nq),       &
                          ADM_kall, ADM_kmin, ADM_kmax )

       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- ', DIAG_name(nq), ': max=', val_max, ', min=', val_min
    enddo

    do nq = 1, TRC_vmax
       val_max = GTL_max( DIAG_var   (:,:,:,DIAG_vmax0+nq), &
                          DIAG_var_pl(:,:,:,DIAG_vmax0+nq), &
                          ADM_kall, ADM_kmin, ADM_kmax      )

       if ( val_max <= 0.0_RP ) then
          nonzero = .false.
       else
          nonzero = .true.
       endif

       val_min = GTL_min( DIAG_var   (:,:,:,DIAG_vmax0+nq), &
                          DIAG_var_pl(:,:,:,DIAG_vmax0+nq), &
                          ADM_kall, ADM_kmin, ADM_kmax, nonzero)

       write(ADM_LOG_FID,'(1x,A,A16,2(A,1PE24.17))') '--- ', TRC_name(nq),  ': max=', val_max, ', min=', val_min
    enddo

    if ( output_io_mode == 'ADVANCED' ) then

       do nq = 1, DIAG_vmax0
          call FIO_output( DIAG_var(:,:,:,nq), basename, desc, '', DIAG_name(nq), DLABEL(nq), '', DUNIT(nq), &
                           FIO_REAL8, layername, 1, ADM_kall, 1, TIME_CTIME, TIME_CTIME                      )
       enddo

       do nq = 1, TRC_vmax
          call FIO_output( DIAG_var(:,:,:,DIAG_vmax0+nq), basename, desc, '', TRC_name(nq), WLABEL(nq), '', WUNIT, &
                           FIO_REAL8, layername, 1, ADM_kall, 1, TIME_CTIME, TIME_CTIME                            )
       enddo

    elseif( output_io_mode == 'LEGACY' ) then

       do l = 1, ADM_lall
          rgnid = ADM_prc_tab(l,ADM_prc_me)
          call MISC_make_idstr(fname,trim(basename),'rgn',rgnid)
          fid = MISC_get_available_fid()
          open( unit   = fid,                 &
                file   = trim(fname),         &
                form   = 'unformatted',       &
                access = 'direct',            &
                recl   = ADM_gall*ADM_kall*8, &
                status = 'unknown'            )

          do nq = 1, DIAG_vmax
             write(fid,rec=nq) DIAG_var(:,:,l,nq)
          enddo

          close(fid)
       enddo

    endif  !--- io_mode

    return
  end subroutine restart_output

end module mod_prgvar
!-------------------------------------------------------------------------------

max_stars_repo_name: hisashiyashiro/nicamdc

max_issues_repo_name: nicamdc-dev/nicamdc

max_forks_repo_name: nicamdc-dev/nicamdc


text: #define included_print_statement print*, "I was included."

max_stars_repo_name: ROCmSoftwarePlatform/gpufort

max_issues_repo_name: odellus/gpufort

max_forks_repo_name: ROCmSoftwarePlatform/gpufort


text:
subroutine sldsup_output(a,itask)
    !-----------------------------------------------------------------------
    !****f* SOLIDS/sldsup_output
    ! NAME
    !    sldsup_output
    ! DESCRIPTION
    !    End of a SOLIDS time step
    !    itask = 0  When timemarching is true. There is output or post-process
    !               of results if required.
    !    itask = 1  When timemarching is false. Output and/or post-process of
    !               results is forced if they have not been written previously.
    ! USES
    !    output
    !    a%FilePostpr%postpr
    ! USED BY
    !    sldsup_endste (itask=0)
    !    sldsup_turnof (itask=1)
    !***
    !-----------------------------------------------------------------------
    use typre
    use Mod_Postpr
    use Mod_SUPSolids
    implicit none
    class(SUPSolidsProblem):: a
    integer(ip)            :: itask
    integer(ip), save      :: dopost(49)

    select case(itask)

    !End of a time step.
   case(0)
      a%pos_alrea=0

   case(1)

   end select

  !Decide which postprocesses need to be done
  call SetupDoPost(itask,a%istep,size(a%npp_stepi),a%npp_inits,a%npp_stepi,a%pos_alrea,dopost)

  !Deviatoric strain
  if (dopost(14) == 1) then
      call a%FilePostpr%postpr(a%devstrain(:,:),'Dev strain',a%istep,a%ctime,a%Mesh,'voigtT')
  endif

  !J2 stresses
  if (dopost(15) == 1) then
      call a%FilePostpr%postpr(a%J2(:),'J2 Stresses',a%istep,a%ctime,a%Mesh)
  endif

  !Residual
  if (dopost(16) == 1) then
      call a%FilePostpr%postgp(a%residualU,'ResidualU',a%istep,a%ctime,a%Mesh)
      call a%FilePostpr%postgp(a%residualP,'ResidualP',a%istep,a%ctime,a%Mesh)
      call a%FilePostpr%postgp(a%residualS,'ResidualS',a%istep,a%ctime,a%Mesh,'voigtT')
  end if

end subroutine sldsup_output

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
Real Function rrtoee(srt)
  rrtoee = 5.
  Return
End Function rrtoee

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X IPModel_TTM_nF
!X
!% TTM_nF module
!%
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"

module IPModel_TTM_nF_module

use error_module
use system_module, only : dp, inoutput, print, verbosity_push_decrement, verbosity_pop, operator(//), lower_case
use dictionary_module
use units_module, only : KCAL_MOL
use linearalgebra_module, only: is_diagonal
use topology_module, only: find_water_monomer
use paramreader_module
use linearalgebra_module
use atoms_types_module
use atoms_module

use mpi_context_module
use QUIP_Common_module

use iso_c_binding, only: C_INT, C_SIZE_T, C_DOUBLE

implicit none

#ifdef HAVE_TTM_NF

integer(C_INT), protected, bind(C,name="POT_QTIP4PF") :: POT_QTIP4PF
integer(C_INT), protected, bind(C,name="POT_TTM2F") :: POT_TTM2F
integer(C_INT), protected, bind(C,name="POT_TTM3F") :: POT_TTM3F
integer(C_INT), protected, bind(C,name="POT_TTM4F") :: POT_TTM4F

interface
   subroutine ttm_nf_energy(pot_type,n,pos,energy) bind(C)
      use iso_c_binding, only: C_INT, C_SIZE_T, C_DOUBLE
      integer(C_INT), value :: pot_type
      integer(C_SIZE_T), value :: n
      real(C_DOUBLE), dimension(*) :: pos
      real(C_DOUBLE) :: energy
   endsubroutine ttm_nf_energy
endinterface

interface
   subroutine ttm_nf_energy_gradient(pot_type,n,pos,energy,gradient) bind(C)
      use iso_c_binding, only: C_INT, C_SIZE_T, C_DOUBLE
      integer(C_INT), value :: pot_type
      integer(C_SIZE_T), value :: n
      real(C_DOUBLE), dimension(*) :: pos
      real(C_DOUBLE) :: energy
      real(C_DOUBLE), dimension(*) :: gradient
   endsubroutine ttm_nf_energy_gradient
endinterface

#endif

private

include 'IPModel_interface.h'

public :: IPModel_TTM_nF
type IPModel_TTM_nF
  !integer :: n_types = 0
  !integer, allocatable :: atomic_num(:), type_of_atomic_num(:)

  real(dp) :: cutoff = 0.0_dp
  integer :: potential_type

  character(len=STRING_LENGTH) :: label

end type IPModel_TTM_nF

logical, private :: parse_in_ip, parse_matched_label
type(IPModel_TTM_nF), private, pointer :: parse_ip

interface Initialise
  module procedure IPModel_TTM_nF_Initialise_str
end interface Initialise

interface Finalise
  module procedure IPModel_TTM_nF_Finalise
end interface Finalise

interface Print
  module procedure IPModel_TTM_nF_Print
end interface Print

interface Calc
  module procedure IPModel_TTM_nF_Calc
end interface Calc

contains

subroutine IPModel_TTM_nF_Initialise_str(this, args_str, param_str)
  type(IPModel_TTM_nF), intent(inout) :: this
  character(len=*), intent(in) :: args_str, param_str

  character(len=STRING_LENGTH) :: potential_type

  type(Dictionary) :: params

  call Finalise(this)

  call initialise(params)
  this%label=''
  call param_register(params, 'label', '', this%label, help_string="No help yet.  This source file was $LastChangedBy$")
  call param_register(params, 'potential_type', 'TTM4F', potential_type, help_string="Type of potential. Allowed values: QTIP4PF, TTM2F, TTM3F, TTM4F$")

  if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='IPModel_TTM_nF_Initialise_str args_str')) then
    call system_abort("IPModel_TTM_nF_Initialise_str failed to parse label from args_str="//trim(args_str))
  endif
  call finalise(params)

  this%cutoff = 2.0_dp
#ifdef HAVE_TTM_NF
  select case(lower_case(trim(potential_type)))
  case("qtip4pf")
     this%potential_type = POT_QTIP4PF
  case("ttm2f")
     this%potential_type = POT_TTM2F
  case("ttm3f")
     this%potential_type = POT_TTM3F
  case("ttm4f")
     this%potential_type = POT_TTM4F
  case default
     call system_abort("IPModel_TTM_nF_Initialise_str: unknown potential_type "//potential_type)
  endselect
#endif

  !  Add initialisation code here

end subroutine IPModel_TTM_nF_Initialise_str

subroutine IPModel_TTM_nF_Finalise(this)
  type(IPModel_TTM_nF), intent(inout) :: this

  ! Add finalisation code here
  this%cutoff = 0.0_dp

  this%label = ''
end subroutine IPModel_TTM_nF_Finalise

subroutine IPModel_TTM_nF_Calc(this, at, e, local_e, f, virial, local_virial, args_str, mpi, error)
   type(IPModel_TTM_nF), intent(inout):: this
   type(Atoms), intent(inout)      :: at
   real(dp), intent(out), optional :: e, local_e(:)
   real(dp), intent(out), optional :: f(:,:), local_virial(:,:)   !% Forces, dimensioned as \texttt{f(3,at%N)}, local virials, dimensioned as \texttt{local_virial(9,at%N)}
   real(dp), intent(out), optional :: virial(3,3)
   character(len=*), optional      :: args_str
   type(MPI_Context), intent(in), optional :: mpi
   integer, intent(out), optional :: error

   integer(C_SIZE_T) :: nWater
   integer :: i, a
   integer, dimension(:,:), allocatable :: water_monomer_index
   real(dp) :: energy
   real(dp), dimension(3) :: lattice
   real(dp), dimension(:), allocatable :: pos, gradient

   INIT_ERROR(error)
#ifndef HAVE_TTM_NF
   RAISE_ERROR('IPModel_TTM_nF_Calc - not linked to the TTM_nF library',error)
#endif

   if (present(e)) e = 0.0_dp

   if (present(local_e)) then
      RAISE_ERROR('IPModel_TTM_nF_Calc - local energies not implemented',error)
      call check_size('Local_E',local_e,(/at%N/),'IPModel_TTM_nF_Calc', error)
      local_e = 0.0_dp
   endif

   if (present(f)) then
      call check_size('Force',f,(/3,at%Nbuffer/),'IPModel_TTM_nF_Calc', error)
      f = 0.0_dp
   end if

   if (present(virial)) then
      RAISE_ERROR('IPModel_TTM_nF_Calc - virials not implemented',error)
      virial = 0.0_dp
   endif

   if (present(local_virial)) then
      RAISE_ERROR('IPModel_TTM_nF_Calc - local virials not implemented',error)
      call check_size('Local_virial',local_virial,(/9,at%Nbuffer/),'IPModel_TTM_nF_Calc', error)
      local_virial = 0.0_dp
   endif

   nWater = count(at%Z==8)
   allocate(water_monomer_index(3,nWater))
   call find_water_monomer(at,water_monomer_index,error=error)

   allocate(pos(3*at%N),gradient(3*at%N))

   do i = 1, nWater
      do a = 1, 3
         pos(a + (i-1)*9)  = at%pos(a,water_monomer_index(1,i))
         pos(a + (i-1)*9 + 3)  = at%pos(a,water_monomer_index(2,i))
         pos(a + (i-1)*9 + 6)  = at%pos(a,water_monomer_index(3,i))
      enddo
   enddo

#ifdef HAVE_TTM_NF
   if(present(f)) then
      call ttm_nf_energy_gradient(this%potential_type,nWater,pos,energy,gradient)
   endif

   if(present(e)) then
      call ttm_nf_energy(this%potential_type,nWater,pos,energy)
   endif
#endif

   if (present(e)) e = energy * KCAL_MOL

   if(present(f)) then
      do i = 1, nWater
         do a = 1, 3
            f(a,water_monomer_index(1,i)) = - gradient(a + (i-1)*9)
            f(a,water_monomer_index(2,i)) = - gradient(a + (i-1)*9 + 3)
            f(a,water_monomer_index(3,i)) = - gradient(a + (i-1)*9 + 6)
         enddo
      enddo
       f = f * KCAL_MOL
   endif

   if(allocated(water_monomer_index)) deallocate(water_monomer_index)
   if(allocated(pos)) deallocate(pos)
   if(allocated(gradient)) deallocate(gradient)

end subroutine IPModel_TTM_nF_Calc

subroutine IPModel_TTM_nF_Print(this, file)
  type(IPModel_TTM_nF), intent(in) :: this
  type(Inoutput), intent(inout),optional :: file

  integer :: ti

  call Print("IPModel_TTM_nF : TTM_nF Potential", file=file)

#ifdef HAVE_TTM_NF
  if (this%potential_type == POT_QTIP4PF) then
     call print("IPModel_TTM_nF : type is QTIP4PF")
  elseif( this%potential_type == POT_TTM2F ) then
     call print("IPModel_TTM_nF : type is TTM2F")
  elseif( this%potential_type == POT_TTM3F ) then
     call print("IPModel_TTM_nF : type is TTM3F")
  elseif( this%potential_type == POT_TTM4F ) then
     call print("IPModel_TTM_nF : type is TTM4F")
  else
     call system_abort("IPModel_TTM_nF_Print: unknown potential_type "//this%potential_type)
  endif
#endif

end subroutine IPModel_TTM_nF_Print

subroutine IPModel_TTM_nF_read_params_xml(this, param_str)
  type(IPModel_TTM_nF), intent(inout), target :: this
  character(len=*), intent(in) :: param_str

  type(xml_t) :: fxml

  if (len(trim(param_str)) <= 0) return

  parse_in_ip = .false.
  parse_matched_label = .false.
  parse_ip => this

  call open_xml_string(fxml, param_str)
  call parse(fxml,  &
    startElement_handler = IPModel_startElement_handler, &
    endElement_handler = IPModel_endElement_handler)
  call close_xml_t(fxml)

end subroutine IPModel_TTM_nF_read_params_xml

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!% XML param reader functions
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
subroutine IPModel_startElement_handler(URI, localname, name, attributes)
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name
  type(dictionary_t), intent(in) :: attributes

  integer :: status
  character(len=STRING_LENGTH) :: value

  if (name == 'TTM_nF_params') then ! new TTM_nF stanza

    if (parse_matched_label) return ! we already found an exact match for this label

    call QUIP_FoX_get_value(attributes, 'label', value, status)
    if (status /= 0) value = ''

    if (len(trim(parse_ip%label)) > 0) then ! we were passed in a label
      if (value == parse_ip%label) then ! exact match
        parse_matched_label = .true.
        parse_in_ip = .true.
      else ! no match
        parse_in_ip = .false.
      endif
    else ! no label passed in
      parse_in_ip = .true.
    endif

    if(parse_in_ip) then
       call finalise(parse_ip)
    endif
  endif

end subroutine IPModel_startElement_handler

subroutine IPModel_endElement_handler(URI, localname, name)
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name

  if (parse_in_ip) then
    if (name == 'TTM_nF_params') then
      parse_in_ip = .false.
    end if
  endif

end subroutine IPModel_endElement_handler

end module IPModel_TTM_nF_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


! NOTES: Variables have subscript suffixes

module md
    implicit none
    private
    public :: Particle, Particles, Simulation

    !--------- Start of type declarations -----------
    !------------------------------------------------
    type Particle
        real :: mass_
        real, dimension(3) :: force_
        real, dimension(3) :: prev_force_
        real, dimension(3) :: position_
        real, dimension(3) :: velocity_

    contains
        procedure, public :: zero_force
        procedure, public :: update_prev_force
        procedure, public :: update_velocity
        procedure, public :: update_position

    end type Particle

    !------------------------------------------------
    type, public :: Particles
        type(Particle), dimension(12), public :: list_

    contains
        procedure, public :: set_velocities => pset_velocities
        procedure, public :: print_positions => pprint_positions
        procedure, public :: calculate_forces
        procedure, public :: update_velocities
        procedure, public :: update_positions

    end type Particles

    !------------------------------------------------
    type, public :: LJPotential
        real, dimension(3), public :: f_

    contains
        procedure, public :: add_force => padd_force

    end type LJPotential

    !------------------------------------------------
    type, public :: Simulation
        type(Particles), public:: particles_
        type(LJPotential), public :: potential_
        integer, private :: num_steps_
        real, private :: time_step_

    contains
        procedure, public :: run => prun

    end type Simulation
    !---------- End of type declarations -----------

    !-------- Start of constructor defns -----------
    !-----------------------------------------------
    interface Particles
        procedure :: new_particles
    end interface Particles

    interface Simulation
        procedure :: new_simulation
    end interface Simulation

    interface LJPotential
        procedure :: new_potential
    end interface LJPotential
    !------------ End of constructors --------------

    !------------ Start of methods -----------------
    !-----------------------------------------------
    contains

        subroutine zero_force(this)
            ! Zero the force on a particle

            class(Particle), intent(inout) :: this
            integer :: k

            do k = 1, 3
                this%force_(k) = 0.0
            enddo

        end subroutine zero_force

        subroutine update_prev_force(this)
            ! Set the previous force on a particle
            ! from the current value of the force

            class(Particle), intent(inout) :: this
            integer :: k

            do k = 1, 3
                this%prev_force_(k) = this%force_(k)
            enddo

        end subroutine update_prev_force

        subroutine update_velocity(this, dt)
            ! Update the current velocity based on the
            ! current and previous force

            class(Particle), intent(inout) :: this
            real, intent(in) :: dt
            integer :: k

            do k = 1, 3
                ! v_k += (f_k(i) + f_k(i-1)) dt / (2M)
                this%velocity_(k) = (this%velocity_(k) &
                                    + ((this%force_(k) + this%prev_force_(k)) / 2.0 &
                                       * dt / this%mass_)                           &
                                    )
            enddo

        end subroutine update_velocity

        subroutine update_position(this, dt)
            ! Update the current position based on the
            ! velocity and acceleration, for a timstep dt

            class(Particle), intent(inout) :: this
            real, intent(in) :: dt
            integer :: k

            ! self.velocity * dt + self.acceleration * dt**2

            do k = 1, 3
                ! x_k += v_k dt  + (f_k / M) dt^2
                this%position_(k) = (this%position_(k)                        &
                                     + (this%velocity_(k) * dt)               &
                                     + (this%force_(k) / this%mass_) * dt*dt  &
                                     )
            enddo

        end subroutine update_position

        type(Particles) function new_particles(position_filename)
            ! Implementation of the constructor
            ! Generate a set of particles given a file with
            ! their positions defined by:
            ! x0  y0  z0
            ! x1  y1  z1

            character(len = 100), intent(in) :: position_filename
            integer :: i, k, io_status
            integer :: read_unit = 9

            open(unit=read_unit,         &
                 action="read",          &
                 file=position_filename, &
                 status="old",           &
                 iostat=io_status)

            if (io_status /= 0) stop "Error opening file"

            do i = 1,12
                read(read_unit, *) &
                        new_particles%list_(i)%position_(1), &
                             new_particles%list_(i)%position_(2), &
                               new_particles%list_(i)%position_(3)

                do k = 1, 3
                    new_particles%list_(i)%prev_force_(k) = 0.0
                enddo
                new_particles%list_(i)%mass_ = 1.0

             enddo

            close(read_unit) ! Close position_filename

            ! call new_particles%set_pair_list()

        end function new_particles

        subroutine pset_velocities(this, filename_)
            ! Set the velocities for each particle from a
            ! file containing velocities for each particle

            class(Particles), intent(inout) :: this
            character(len = 100), intent(in) :: filename_
            integer :: i
            integer :: read_unit = 9

            open(unit=read_unit, action="read", file=filename_, status="old")

            do i = 1,size(this%list_)
                read(read_unit, *) &
                        this%list_(i)%velocity_(1), &
                             this%list_(i)%velocity_(2), &
                               this%list_(i)%velocity_(3)
            enddo

            close(read_unit)

        end subroutine pset_velocities

        subroutine pprint_positions(this)
            ! Print the positions into a .txt file

            class(Particles), intent(in) :: this

            integer :: print_unit = 10
            integer :: i

            open(unit=print_unit,            &
                 action="write",             &
                 file="final_positions.txt", &
                 status="new")

            do i = 1, size(this%list_)
                write(print_unit, *)   &
                    this%list_(i)%position_(1),            &
                        this%list_(i)%position_(2),        &
                            this%list_(i)%position_(3)

            enddo

            close(print_unit, status="keep")

        end subroutine pprint_positions

        subroutine calculate_forces(this, potential_)
            ! Calculate the forces on all the particles
            class(Particles), intent(inout) :: this
            class(LJPotential), intent(in) :: potential_
            integer :: i, j

            do i = 1, size(this%list_)
                call this%list_(i)%update_prev_force()
                call this%list_(i)%zero_force()

                do j = 1, size(this%list_)
                    if (i /= j) then

                        call potential_%add_force(this%list_(i), &
                                                this%list_(j))
                    end if
                enddo

            enddo

        end subroutine calculate_forces

        subroutine update_velocities(this, dt)
            ! Update the velocities of all particles

            class(Particles), intent(inout) :: this
            real, intent(in) :: dt
            integer :: i

            do i = 1, size(this%list_)
                call this%list_(i)%update_velocity(dt)
            enddo

        end subroutine update_velocities

        subroutine update_positions(this, dt)
            ! Update the positions of all the particles

            class(Particles), intent(inout) :: this
            real, intent(in) :: dt
            integer :: i

            do i = 1, size(this%list_)
                call this%list_(i)%update_position(dt)
            enddo

        end subroutine update_positions

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

        type(LJPotential) function new_potential(epsilon_, &
                                                 sigma_)
            ! Construct a new LJ potential from the
            ! well depth and particle interaction range

            real, intent(in) :: epsilon_, sigma_

            new_potential%f_(1) = epsilon_/2.0
            new_potential%f_(2) = 12.0 * sigma_**12
            new_potential%f_(3) = -6.0 * sigma_**6

        end function new_potential

        subroutine padd_force(this, particle_i_, particle_j_)
            ! Add the force on particle i due to the presence
            ! of particle j
            class(LJPotential), intent(in) :: this
            class(Particle), intent(inout) :: particle_i_
            class(Particle), intent(in) :: particle_j_

            real, dimension(5) :: v

            v(1) = particle_i_%position_(1) - particle_j_%position_(1)  ! dx
            v(2) = particle_i_%position_(2) - particle_j_%position_(2)  ! dy
            v(3) = particle_i_%position_(3) - particle_j_%position_(3)  ! dz

            v(4) = sqrt(v(1) * v(1) + v(2) * v(2) + v(3) * v(3))        ! r
            v(5) = this%f_(1) * (this%f_(2) * v(4)**(-14) + this%f_(3) * v(4)**(-8))

            particle_i_%force_(1) = particle_i_%force_(1) + v(5) * v(1)   ! f_x + c * dx
            particle_i_%force_(2) = particle_i_%force_(2) + v(5) * v(2)   ! f_y + c * dy
            particle_i_%force_(3) = particle_i_%force_(3) + v(5) * v(3)   ! f_z + c * dz

        end subroutine padd_force

        type(Simulation) function new_simulation(particles_, &
                                                 potential_, &
                                                 num_steps_, &
                                                 time_step_)
            ! Create a new simulation from a set of particles
            ! along with a potential and

            class(Particles), intent(in) :: particles_
            class(LJPotential), intent(in) :: potential_
            integer, intent(in) :: num_steps_
            real, intent(in) :: time_step_

            new_simulation%particles_ = particles_
            new_simulation%potential_ = potential_

            if (num_steps_ <= 0) stop "Number of simulation steps must be positive"
            new_simulation%num_steps_ = num_steps_

            if (time_step_ <= 0.0) stop "Timestep must be positive"
            new_simulation%time_step_ = time_step_

        end function new_simulation

        subroutine prun(this)
            ! Run the simulation using the defined number
            ! of steps

            class(Simulation), intent(inout) :: this
            integer :: step

            call this%particles_%calculate_forces(this%potential_)

            do step = 0, this%num_steps_

                call this%particles_%update_positions(this%time_step_)
                call this%particles_%calculate_forces(this%potential_)
                call this%particles_%update_velocities(this%time_step_)

            enddo

        end subroutine prun

end module md

program run_md
    use md
    implicit none

    !----------------------------------------- Type definitions
    character(len=100) :: filename_
    type(Particles) :: particles_
    type(Simulation) :: simulation_
    type(LJPotential) :: potential_
    !-----------------------------------------

    ! --------- Initalise particles ----------
    filename_ = "positions.txt"
    particles_ = Particles(position_filename=filename_)

    filename_ = "velocities.txt"
    call particles_%set_velocities(filename_)

    ! --------- Initalise potential ----------
    potential_ = LJPotential(epsilon_=100.0, &
                             sigma_=1.7)

    ! --------- Initalise simulation ---------
    simulation_ = Simulation(particles_=particles_, &
                             potential_=potential_, &
                             num_steps_=10000,      &
                             time_step_=0.01)

    call simulation_%run()
    call simulation_%particles_%print_positions()

end program


text:
	SUBROUTINE FFDRIV  ( nlvl, clvl, lincol, lintyp, linwid, linlbl,
     +			     linflg, zgrid, iret )
C************************************************************************
C* FFDRIV								*
C*									*
C* This subroutine draws filled contours.				*
C*									*
C* FFDRIV  ( NLVL, CLVL, LINCOL, LINTYP, LINWID, LINLBL, LINFLG, ZGRID, *
C*           IRET )							*
C*									*
C* Input parameters:							*
C*	NLVL		INTEGER		Number of contour levels	*
C*	CLVL   (NLVL)	REAL		Contour levels			*
C*	LINCOL (NLVL)	INTEGER		Contour colors			*
C*	LINTYP (NLVL)	INTEGER		Contour line types		*
C*	LINWID (NLVL)	INTEGER		Contour line widths		*
C*	LINLBL (NLVL)	INTEGER		Contour label types		*
C*	LINFLG		LOGICAL		Flag to draw lines 		*
C*									*
C* Work parameters:                                                     *
C*      ZGRID (KX*KY)   REAL            Array for grid points           *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/NMC	11/91						*
C* D.W.Plummer/NCEP	 1/97	Remove call to GEPLOT prior to RETURN	*
C* T. Lee/SAIC		10/01	Query and restore fill attributes	*
C* S. Gilbert/NCEP      03/07   Added work arrays to call sequence      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'CONTUR.CMN'
C*
	REAL		clvl (*), zgrid (*)
	INTEGER		lincol (*), lintyp (*), linwid (*), linlbl (*)
	LOGICAL		linflg
C*
	REAL		box (4), xpt (250), ypt (250)
	LOGICAL		onelev
C------------------------------------------------------------------------
	iret = 0
C
C*	Query current fill type.
C
	CALL GQFILL ( szfil, iftyp, iret )
C
C*	Loop through all the grid boxes.
C
	DO  jll = 1, jsize - 1
	  DO  ill = 1, isize - 1
C
C*	    Get points on this box and check to see if the entire box
C*	    is one color.
C
	    CALL FGTBOX  ( ill, jll, box, xx1, yy1, xx2, yy2, zgrid,
     +                     iret )
C
C*	    If the return code is not 0, there is missing data.
C
	    IF  ( iret .eq. 0 )  THEN
		CALL FGTCOR  ( box, nlvl, clvl, icornr, levl, onelev,
     +			       ier )
	    END IF
C
C*	    Do nothing if there is missing data.
C
	    IF  ( iret .ne. 0 )  THEN
C
C*	      Fill entire box if it is all one color.
C
	     ELSE IF  ( onelev .and. ( .not. linflg ) )  THEN
C
C*	      Plot the entire box.
C
	      ipt = 1
	      xpt (ipt) = xx1
	      ypt (ipt) = yy1
	      DO  jj = 1, numsub - 2
		ipt = ipt + 1
		xpt (ipt) = xx1
		ypt (ipt) = FLOAT ( jj ) * fincxy + yy1
	      END DO
	      ipt = ipt + 1
	      xpt (ipt) = xx1
	      ypt (ipt) = yy2
	      DO  ii = 1, numsub - 2
		ipt = ipt + 1
		xpt (ipt) = FLOAT ( ii ) * fincxy + xx1
		ypt (ipt) = yy2
	      END DO
	      ipt = ipt + 1
	      xpt (ipt) = xx2
	      ypt (ipt) = yy2
	      DO  jj = numsub - 2, 1, -1
		ipt = ipt + 1
		xpt (ipt) = xx2
		ypt (ipt) = FLOAT ( jj ) * fincxy + yy1
	      END DO
	      ipt = ipt + 1
	      xpt (ipt) = xx2
	      ypt (ipt) = yy1
	      DO  ii = numsub - 2, 1, -1
		ipt = ipt + 1
		xpt (ipt) = FLOAT ( ii ) * fincxy + xx1
		ypt (ipt) = yy1
	      END DO
	      CALL FFPLOT  ( ipt, xpt, ypt, linflg, lincol (levl),
     +			     lintyp (levl), linwid (levl), ier )
C
C*	      If box is not subdivided, fill in the whole box.
C
	     ELSE IF  ( numsub .le. 2 )  THEN
	      CALL FILBOX  ( box, xx1, yy1, xx2, yy2, nlvl, clvl,
     +			     lincol, lintyp, linwid, icornr, levl,
     +			     linflg, ier )
C
C*	      Subdivide the box.
C
	     ELSE
C
C*	      Set up the coefficients for this box.
C
	      CALL CCOEFF  ( ill, jll, zgrid, ier )
C
C*	      Loop through all the subboxes.
C
	      DO  jsub = 1, numsub - 1
		DO  isub = 1, numsub - 1
C
C*		  Get the points on this subbox and check for a single
C*		  color.
C
		  CALL FGTSUB  ( ill, jll, isub, jsub, box, xx1, yy1,
     +				 xx2, yy2, zgrid, ier )
		  CALL FGTCOR  ( box, nlvl, clvl, icornr, levl, onelev,
     +				 ier )
C
C*		  Fill the entire box if it is all one color.
C
		  IF  ( onelev )  THEN
C
C*		    Plot the entire box.
C
		    xpt (1) = xx1
		    ypt (1) = yy1
		    xpt (2) = xx1
		    ypt (2) = yy2
		    xpt (3) = xx2
		    ypt (3) = yy2
	      	    xpt (4) = xx2
		    ypt (4) = yy1
		    ipt = 4
		    CALL FFPLOT  ( ipt, xpt, ypt, linflg, lincol (levl),
     +				   lintyp (levl), linwid (levl), ier )
C
C*		    Fill in this box.
C
		   ELSE
		    CALL FILBOX (box, xx1, yy1, xx2, yy2, nlvl, clvl,
     +				 lincol, lintyp, linwid, icornr,
     +				 levl, linflg, ier )
		  END IF
		END DO
	      END DO
	    END IF
	  END DO
	END DO
999	continue
C
C*	Restore fill type setup.
C
	CALL GSFILL ( szfil, iftyp, iret )
C
C*	End plotting.
C
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak