!
! CDDL HEADER START
!
! The contents of this file are subject to the terms of the Common Development
! and Distribution License Version 1.0 (the "License").
!
! You can obtain a copy of the license at
! http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
! specific language governing permissions and limitations under the License.
!
! When distributing Covered Code, include this CDDL HEADER in each file and
! include the License file in a prominent location with the name LICENSE.CDDL.
! If applicable, add the following below this CDDL HEADER, with the fields
! enclosed by brackets "[]" replaced with your own identifying information:
!
! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
!
! CDDL HEADER END
!

!
! Copyright (c) 2013--2018, Regents of the University of Minnesota.
! All rights reserved.
!
! Contributors:
!    Ryan S. Elliott
!    Ellad B. Tadmor
!    Valeriu Smirichinski
!


!*******************************************************************************
!**
!**  MODULE ex_model_Ar_P_MLJ_NEIGH_PURE_H_F
!**
!**  Lennard-Jones pair potential model for argon
!**  (modified by adding a quadratic function to have a smooth cutoff)
!**
!**  Release: This file is part of the kim-api-v1.9.3 package.
!**
!*******************************************************************************

#include "KIM_API_status.h"
#define THIS_FILE_NAME __FILE__

module ex_model_Ar_P_MLJ_NEIGH_PURE_H_F
  use, intrinsic :: iso_c_binding
  use KIM_API_F03
  implicit none

  save
  private
  integer(c_int), parameter :: cd = c_double ! used for literal constants

  public Compute_Energy_Forces
  public ReInit
  public Destroy

contains

!-------------------------------------------------------------------------------
!
! Compute energy and forces on particles from the positions.
!
!-------------------------------------------------------------------------------
  integer(c_int) function Compute_Energy_Forces(pkim) bind(c)
    implicit none

    !-- Transferred variables
    type(c_ptr), intent(in)  :: pkim

    !-- Local variables
    integer(c_int), parameter :: DIM=3
    real(c_double) r,Rsqij,phi,dphi,d2phi,dEidr
    integer(c_int) idum
    integer(c_int) i,j,jj,numnei,part,part_ret
    character(len=KIM_KEY_STRING_LENGTH) :: NBC_Method
    integer(c_int) HalfOrFull ! 1 - half, 2 - full
    real(c_double), dimension(DIM) :: Rij
    type(c_ptr) :: pRij
    
    !-- KIM variables
    integer(c_int), pointer :: numberOfParticles;  type(c_ptr) :: pnParts
    integer(c_int), pointer :: nOfSpecies;         type(c_ptr) :: pnOfSpecies
    integer(c_int), pointer :: particleSpecies(:); type(c_ptr) :: pparticleSpecies
    real(c_double), pointer :: model_cutoff;       type(c_ptr) :: pcutoff
    real(c_double), pointer :: model_epsilon;      type(c_ptr) :: pepsilon
    real(c_double), pointer :: model_sigma;        type(c_ptr) :: psigma
    real(c_double), pointer :: model_cutnorm;      type(c_ptr) :: pcutnorm
    real(c_double), pointer :: model_A;            type(c_ptr) :: pA
    real(c_double), pointer :: model_B;            type(c_ptr) :: pB
    real(c_double), pointer :: model_C;            type(c_ptr) :: pC
    real(c_double), pointer :: model_sigmasq;      type(c_ptr) :: psigmasq
    real(c_double), pointer :: model_cutsq;        type(c_ptr) :: pcutsq
    real(c_double), pointer :: energy;             type(c_ptr) :: penergy
    real(c_double), pointer :: coor(:,:);          type(c_ptr) :: pcoor
    real(c_double), pointer :: force(:,:);         type(c_ptr) :: pforce
    real(c_double), pointer :: enepot(:);          type(c_ptr) :: penepot
    real(c_double), pointer :: virial(:);          type(c_ptr) :: pvirial
    integer(c_int) comp_energy, comp_force, comp_enepot, comp_virial
    integer(c_int), pointer :: numContrib;              type(c_ptr) :: pnumContrib
    integer(c_int), pointer :: nei1part(:);             type(c_ptr) :: pnei1part
    
    ! Check to see if we have been asked to compute the energy, forces,
    ! energyperpart, and virial
    !
    call kim_api_getm_compute(pkim, Compute_Energy_Forces, &
         "energy",         comp_energy, 1, &
         "forces",         comp_force,  1, &
         "particleEnergy", comp_enepot, 1, &
         "virial",         comp_virial, 1)
    if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_compute",   &
                                   Compute_Energy_Forces)
       return
    endif

    ! Unpack data from KIM object
    !
    call kim_api_getm_data(pkim,  Compute_Energy_Forces,          &
         "numberOfParticles",     pnParts,           1,           &
         "numberOfSpecies",       pnOfSpecies,       1,           &
         "particleSpecies",       pparticleSpecies,  1,           &
         "cutoff",                pcutoff,           1,           &
         "coordinates",           pcoor,             1,           &
         "energy",                penergy,           comp_energy, &
         "forces",                pforce,            comp_force,  &
         "particleEnergy",        penepot,           comp_enepot, &
         "virial",                pvirial,           comp_virial)
    if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_data", Compute_Energy_Forces)
       return
    endif

    call c_f_pointer(pnParts,          numberOfParticles)
    call c_f_pointer(pnOfSpecies,      nOfSpecies)
    call c_f_pointer(pparticleSpecies, particleSpecies, [numberOfParticles])
    call c_f_pointer(pcutoff,          model_cutoff)
    call c_f_pointer(pcoor,            coor,            [DIM,numberOfParticles])
    if (comp_energy.eq.1) call c_f_pointer(penergy, energy)
    if (comp_force.eq.1)  call c_f_pointer(pforce,  force,  &
                                                        [DIM,numberOfParticles])
    if (comp_enepot.eq.1) call c_f_pointer(penepot, enepot, [numberOfParticles])
    if (comp_virial.eq.1) call c_f_pointer(pvirial, virial, [6])

    call kim_api_getm_data(pkim, Compute_Energy_Forces, &
         "PARAM_FREE_epsilon",  pepsilon,      1,         &
         "PARAM_FREE_sigma",    psigma,        1,         &
         "PARAM_FIXED_cutnorm", pcutnorm,      1,         &
         "PARAM_FIXED_A",       pA,            1,         &
         "PARAM_FIXED_B",       pB,            1,         &
         "PARAM_FIXED_C",       pC,            1,         &
         "PARAM_FIXED_sigmasq", psigmasq,      1,         &
         "PARAM_FIXED_cutsq",   pcutsq,        1)
    if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_data", Compute_Energy_Forces)
       return
    endif

    call c_f_pointer(pepsilon, model_epsilon)
    call c_f_pointer(psigma,   model_sigma)
    call c_f_pointer(pcutnorm, model_cutnorm)
    call c_f_pointer(pA,       model_A)
    call c_f_pointer(pB,       model_B)
    call c_f_pointer(pC,       model_C)
    call c_f_pointer(psigmasq, model_sigmasq)
    call c_f_pointer(pcutsq,   model_cutsq)

    ! Check to be sure that the species are correct by comparing
    ! the provided species codes to the value given here (which should
    ! be the same as that given in the .kim file).
    !
    Compute_Energy_Forces = KIM_STATUS_FAIL ! assume an error
    do i = 1,numberOfParticles
       if (.not. (particleSpecies(i) .eq. 1)) then
          idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                      "Wrong Species", Compute_Energy_Forces)
          return
       endif
    enddo
    Compute_Energy_Forces = KIM_STATUS_OK ! everything is ok


    ! Initialize potential energies, forces, virial term
    !
    if (comp_enepot.eq.1) enepot = 0.0_cd
    if (comp_energy.eq.1) energy = 0.0_cd
    if (comp_force.eq.1)  force  = 0.0_cd
    if (comp_virial.eq.1) virial = 0.0_cd
    
    
    ! determine whether half or full neighbor lists are being used
    Compute_Energy_Forces = kim_api_get_nbc_method(pkim, NBC_Method)
    if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME,   &
                                   "kim_api_get_nbc_method", &
                                   Compute_Energy_Forces)
       return
    endif
    if (index(NBC_Method,"NEIGH_PURE_H").eq.1) then
       HalfOrFull = 1
       ! get numberContributingParticles
       pnumContrib = kim_api_get_data(pkim,"numberContributingParticles", &
                                      Compute_Energy_Forces)
       if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
          idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                      "kim_api_get_data", Compute_Energy_Forces)
          return
       endif
       call c_f_pointer(pnumContrib, numContrib)
    elseif (index(NBC_Method,"NEIGH_PURE_F").eq.1) then
       HalfOrFull = 2
       allocate( numContrib )
    else
       Compute_Energy_Forces = KIM_STATUS_FAIL
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "Unsupported NBC type",   &
                                   Compute_Energy_Forces)
       return
    endif

    !  Compute energy and forces
    !
    do i = 1,numberOfParticles

       ! Get neighbors for particle i
       !
       part = i ! request neighbors for particle i

       Compute_Energy_Forces = kim_api_get_neigh(pkim,1,part,part_ret,numnei, &
                                                 pnei1part,pRij)
       if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
          idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                      "kim_api_get_neigh",      &
                                      Compute_Energy_Forces)
          return
       endif
       call c_f_pointer(pnei1part, nei1part, [numnei])

       ! Loop over the neighbors of particle i
       !
       do jj = 1, numnei
          j = nei1part(jj)
          Rij(:) = coor(:,j) - coor(:,i)    ! distance vector between i j
          Rsqij = dot_product(Rij,Rij)      ! compute square distance
          if ( Rsqij < model_cutsq ) then   ! particles are interacting?
             r = sqrt(Rsqij)                ! compute distance
             call pair(model_epsilon,model_sigma,model_A,model_B, model_C, &
                  r,phi,dphi,d2phi)         ! compute pair potential
             if ((HalfOrFull.eq.1) .and. &
                 (j .le. numContrib)) then  ! HALF mode
                dEidr = dphi                !      double contribution
             else                           ! FULL mode
                dEidr = 0.5_cd*dphi         !      regular contribution
             endif
             if (comp_enepot.eq.1) then
                enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
                if ((HalfOrFull.eq.1) .and. &
                     (j .le. numContrib)) then
                   enepot(j) = enepot(j) + 0.5_cd*phi  ! (i and j share it)
                endif
             endif
             if (comp_energy.eq.1) then
               if ((HalfOrFull.eq.1) .and. &
                   (j .le. numContrib)) then           ! HALF mode
                   energy = energy + phi               ! half neigh case
                else
                   energy = energy + 0.5_cd*phi        ! full neigh case
                endif
             endif
             if (comp_virial.eq.1) then                ! accumul. virial
                virial(1) = virial(1) + Rij(1)*Rij(1)*dEidr/r
                virial(2) = virial(2) + Rij(2)*Rij(2)*dEidr/r
                virial(3) = virial(3) + Rij(3)*Rij(3)*dEidr/r
                virial(4) = virial(4) + Rij(2)*Rij(3)*dEidr/r
                virial(5) = virial(5) + Rij(1)*Rij(3)*dEidr/r
                virial(6) = virial(6) + Rij(1)*Rij(2)*dEidr/r
             endif
             if (comp_force.eq.1) then
                force(:,i) = force(:,i) + dEidr*Rij/r  ! accumulate forces on i
                force(:,j) = force(:,j) - dEidr*Rij/r  ! accumulate forces on j
             endif
          endif
       enddo
    enddo

    if (HalfOrFull.eq.2) deallocate( numContrib )
  end function Compute_Energy_Forces

!-------------------------------------------------------------------------------
!
! Pair potential: Lennard-Jones with smooth cutoff imposed by Ar^2 + Br + C
!
!-------------------------------------------------------------------------------
  subroutine pair(epsilon,sigma,A,B,C,r,phi,dphi,d2phi)
    implicit none
    integer(c_int), parameter :: cd = c_double ! used for literal constants

    !-- Transferred variables
    real(c_double), intent(in)  :: epsilon, sigma, A, B, C
    real(c_double), intent(in)  :: r
    real(c_double), intent(out) :: phi, dphi, d2phi

    !-- Local variables
    real(c_double) rsq,sor,sor6,sor12

    rsq  = r*r             !  r^2
    sor  = sigma/r         !  (sig/r)
    sor6 = sor*sor*sor     !
    sor6 = sor6*sor6       !  (sig/r)^6
    sor12= sor6*sor6       !  (sig/r)^12

    phi   =  4.0_cd*epsilon*(sor12-sor6) + A*rsq + B*r + C
    dphi  = 24.0_cd*epsilon*(-2.0_cd*sor12+sor6)/r  + 2.0_cd*A*r + B
    d2phi = 24.0_cd*epsilon*(26.0_cd*sor12-7.0_cd*sor6)/rsq + 2.0_cd*A

  end subroutine pair

!-------------------------------------------------------------------------------
!
! Model reinitialization routine
!
!-------------------------------------------------------------------------------
  integer(c_int) function ReInit(pkim) bind(c)
    implicit none
    integer(c_int), parameter :: cd = c_double ! used for literal constants

    !-- Transferred variables
    type(c_ptr), intent(in) :: pkim

    !-- Local variables
    real(c_double), pointer :: model_cutoff;  type(c_ptr) :: pcutoff
    real(c_double), pointer :: model_epsilon; type(c_ptr) :: pepsilon
    real(c_double), pointer :: model_sigma;   type(c_ptr) :: psigma
    real(c_double), pointer :: model_Pcutoff; type(c_ptr) :: pparamcut
    real(c_double), pointer :: model_cutnorm; type(c_ptr) :: pcutnorm
    real(c_double), pointer :: model_A;       type(c_ptr) :: pA
    real(c_double), pointer :: model_B;       type(c_ptr) :: pB
    real(c_double), pointer :: model_C;       type(c_ptr) :: pC
    real(c_double), pointer :: model_sigmasq; type(c_ptr) :: psigmasq
    real(c_double), pointer :: model_cutsq;   type(c_ptr) :: pcutsq
    integer idum

    ! Get (changed) parameters from KIM object ---------------------------------

    ! get stuff from KIM object
    call kim_api_getm_data(pkim, ReInit, &
         "PARAM_FREE_sigma",   psigma,    1, &
         "PARAM_FREE_epsilon", pepsilon,  1, &
         "PARAM_FREE_cutoff",  pparamcut, 1)
    if (ReInit.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_data", ReInit)
       return
    endif

    call c_f_pointer(psigma,    model_sigma)
    call c_f_pointer(pepsilon,  model_epsilon)
    call c_f_pointer(pparamcut, model_Pcutoff)

    ! Set new values in KIM object ---------------------------------------------

    ! Set stuff in KIM object
    call kim_api_getm_data(pkim, ReInit, &
         "cutoff",              pcutoff,  1, &
         "PARAM_FIXED_cutnorm", pcutnorm, 1, &
         "PARAM_FIXED_A",       pA,       1, &
         "PARAM_FIXED_B",       pB,       1, &
         "PARAM_FIXED_C",       pC,       1, &
         "PARAM_FIXED_sigmasq", psigmasq, 1, &
         "PARAM_FIXED_cutsq",   pcutsq,   1)
    if (ReInit.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_data", ReInit)
       return
    endif

    call c_f_pointer(pcutoff,  model_cutoff)
    call c_f_pointer(pcutnorm, model_cutnorm)
    call c_f_pointer(pA,       model_A)
    call c_f_pointer(pB,       model_B)
    call c_f_pointer(pC,       model_C)
    call c_f_pointer(psigmasq, model_sigmasq)
    call c_f_pointer(pcutsq,   model_cutsq)

    model_cutoff = model_Pcutoff
    model_cutnorm = model_cutoff/model_sigma
    model_A = 12.0_cd*model_epsilon*(-26.0_cd + 7.0_cd*model_cutnorm**6)/ &
         (model_cutnorm**14*model_sigma**2)
    model_B = 96.0_cd*model_epsilon*(7.0_cd-2.0_cd*model_cutnorm**6)/     &
         (model_cutnorm**13*model_sigma)
    model_C = 28.0_cd*model_epsilon*(-13.0_cd+4.0_cd*model_cutnorm**6)/   &
         (model_cutnorm**12)
    model_sigmasq = model_sigma**2
    model_cutsq = model_cutoff**2

    ReInit = KIM_STATUS_OK
    return

  end function ReInit

!-------------------------------------------------------------------------------
!
! Model destroy routine
!
!-------------------------------------------------------------------------------
  integer(c_int) function Destroy(pkim) bind(c)
    implicit none

    !-- Transferred variables
    type(c_ptr), intent(in) :: pkim

    !-- Local variables
    real(c_double), pointer :: model_epsilon; type(c_ptr) :: pepsilon
    real(c_double), pointer :: model_sigma;   type(c_ptr) :: psigma
    real(c_double), pointer :: model_Pcutoff; type(c_ptr) :: pparamcut
    real(c_double), pointer :: model_cutnorm; type(c_ptr) :: pcutnorm
    real(c_double), pointer :: model_A;       type(c_ptr) :: pA
    real(c_double), pointer :: model_B;       type(c_ptr) :: pB
    real(c_double), pointer :: model_C;       type(c_ptr) :: pC
    real(c_double), pointer :: model_sigmasq; type(c_ptr) :: psigmasq
    real(c_double), pointer :: model_cutsq;   type(c_ptr) :: pcutsq
    integer(c_int) idum


    call kim_api_getm_data(pkim, Destroy, &
         "PARAM_FREE_sigma",    psigma,    1, &
         "PARAM_FREE_epsilon",  pepsilon,  1, &
         "PARAM_FREE_cutoff",   pparamcut, 1, &
         "PARAM_FIXED_cutnorm", pcutnorm,  1, &
         "PARAM_FIXED_A",       pA,        1, &
         "PARAM_FIXED_B",       pB,        1, &
         "PARAM_FIXED_C",       pC,        1, &
         "PARAM_FIXED_sigmasq", psigmasq,  1, &
         "PARAM_FIXED_cutsq",   pcutsq,    1)
    if (Destroy.lt.KIM_STATUS_OK) then
       idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                   "kim_api_getm_data", Destroy)
       return
    endif

    call c_f_pointer(psigma,    model_sigma)
    call c_f_pointer(pepsilon,  model_epsilon)
    call c_f_pointer(pparamcut, model_Pcutoff)
    call c_f_pointer(pcutnorm, model_cutnorm)
    call c_f_pointer(pA,       model_A)
    call c_f_pointer(pB,       model_B)
    call c_f_pointer(pC,       model_C)
    call c_f_pointer(psigmasq, model_sigmasq)
    call c_f_pointer(pcutsq,   model_cutsq)

    deallocate( model_sigma   )
    deallocate( model_epsilon )
    deallocate( model_Pcutoff )
    deallocate( model_cutnorm )
    deallocate( model_A       )
    deallocate( model_B       )
    deallocate( model_C       )
    deallocate( model_sigmasq )
    deallocate( model_cutsq   )

    Destroy = KIM_STATUS_OK
    return

  end function Destroy

end module ex_model_Ar_P_MLJ_NEIGH_PURE_H_F


!-------------------------------------------------------------------------------
!
! Model initialization routine (REQUIRED)
!
!-------------------------------------------------------------------------------
integer(c_int) function model_init(pkim) bind(c)
  use, intrinsic :: iso_c_binding
  use ex_model_Ar_P_MLJ_NEIGH_PURE_H_F
  use KIM_API_F03
  implicit none
  integer(c_int), parameter :: cd = c_double ! used for literal constants

  !-- Transferred variables
  type(c_ptr), intent(in) :: pkim

  !-- Local variables
  integer(c_int), parameter :: one=1
  real(c_double), pointer :: model_cutoff;  type(c_ptr) :: pcutoff
  real(c_double), pointer :: model_epsilon
  real(c_double), pointer :: model_sigma
  real(c_double), pointer :: model_Pcutoff
  real(c_double), pointer :: model_cutnorm
  real(c_double), pointer :: model_A
  real(c_double), pointer :: model_B
  real(c_double), pointer :: model_C
  real(c_double), pointer :: model_sigmasq
  real(c_double), pointer :: model_cutsq
  integer(c_int) ier, idum

  ! store pointers in KIM object
  call kim_api_setm_method(pkim, ier, &
       "compute", one, c_funloc(Compute_Energy_Forces), 1, &
       "reinit",  one, c_funloc(ReInit),                1, &
       "destroy", one, c_funloc(Destroy),               1)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_setm_method", ier)
     goto 42
  endif

  ! store model cutoff in KIM object
  pcutoff =  kim_api_get_data(pkim,"cutoff",ier)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_get_data", ier)
     goto 42
  endif
  call c_f_pointer(pcutoff, model_cutoff)
  model_cutoff = 8.15_cd ! cutoff distance in angstroms

  ! Allocate memory for parameters
  allocate( model_epsilon )
  allocate( model_sigma   )
  allocate( model_Pcutoff )
  allocate( model_cutnorm )
  allocate( model_A       )
  allocate( model_B       )
  allocate( model_C       )
  allocate( model_sigmasq )
  allocate( model_cutsq   )

  model_sigma = 3.40_cd ! LJ Argon sigma in angstroms

  model_epsilon = 0.0104_cd ! LJ Argon epsilon in eV

  model_Pcutoff = model_cutoff

  model_cutnorm = model_cutoff/model_sigma

  model_A = 12.0_cd*model_epsilon*(-26.0_cd + 7.0_cd*model_cutnorm**6)/ &
       (model_cutnorm**14*model_sigma**2)

  model_B = 96.0_cd*model_epsilon*(7.0_cd-2.0_cd*model_cutnorm**6)/     &
       (model_cutnorm**13*model_sigma)

  model_C = 28.0_cd*model_epsilon*(-13.0_cd+4.0_cd*model_cutnorm**6)/   &
       (model_cutnorm**12)

  model_sigmasq = model_sigma**2

  model_cutsq = model_cutoff**2


  call kim_api_setm_data(pkim, ier, &
       "PARAM_FREE_sigma",    one, c_loc(model_sigma),    1, &
       "PARAM_FREE_epsilon",  one, c_loc(model_epsilon),  1, &
       "PARAM_FREE_cutoff",   one, c_loc(model_Pcutoff),  1, &
       "PARAM_FIXED_cutnorm", one, c_loc(model_cutnorm),  1, &
       "PARAM_FIXED_A",       one, c_loc(model_A),        1, &
       "PARAM_FIXED_B",       one, c_loc(model_B),        1, &
       "PARAM_FIXED_C",       one, c_loc(model_C),        1, &
       "PARAM_FIXED_sigmasq", one, c_loc(model_sigmasq),  1, &
       "PARAM_FIXED_cutsq",   one, c_loc(model_cutsq),    1)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_setm_data", ier)
     goto 42
  endif

  ier = KIM_STATUS_OK
42 continue
  model_init = ier
  return

end function model_init
