Settings of condensation mass source terms

Introduction

Source terms modelling condensation inside the fluid domain on internal metal structures and at the boundaries can be set respectively through the subroutines cs_user_metal_structures_source_terms and cs_user_boundary_mass_source_terms.

Source terms for condensation on internal metal structures

This model can be enabled in the subroutine usppmo in the file cs_user_parameters.f90 as follows:

if ( ippmod(igmix).ge.0 ) then
  ! Specific condensation modelling
  !      if = -1 module not activated
  !      if =  0 condensation source terms with metal
  !                               structures activate
  icondv = -1
endif

The setting of the condensation source terms is then done in the subroutine cs_user_metal_structures_source_terms as follows below.

The following variables need to be declared:

integer          icmst
integer          ifac, iel, iscal
integer          ivarh
integer          izone
integer          f_id

double precision hvap, tk

type(gas_mix_species_prop) s_h2o_g

double precision, dimension(:), pointer :: cpro_cp
double precision, dimension(:), pointer :: cvar_h

Necessary species physical properties can be retrieved as follows:

call field_get_id_try("y_h2o_g", f_id)
if (f_id.ne.-1) &
  call field_get_key_struct_gas_mix_species_prop(f_id, s_h2o_g)

The zones on which the condensation mass source term will be imposed can be defined as follows:

!===============================================================================
! Select the cells which are associated to the metal structures volume
! with the subroutine getcel
!===============================================================================

izone = 0
met_znb = 1

do izmet = 1 , met_znb

  if (izmet.eq.1) then
    ! Cells associated to the geometric zone
    call getcel('z > -7.0d0 and z < 53.d0', ncmast, ltmast)

    izone = izone + 1

    do icmst = 1, ncmast
      iel = ltmast(icmst)
      izmast(iel) = izone
    enddo
  endif

enddo

Modelling of the metal side (wall thermal model and metal properties can then be specified as follows:

!===============================================================================
! Parameters padding of the wall thermal model and condensation model
! ------------------------------------------------------------------
! The condensation model can be used alone with a constant temperature
! specified by the user at the cold wall (tmet=tmet0 in this case)
! or together with a 0-D thermal model. In the latter case, the two models are
! coupled.
!===============================================================================

if (icondv.eq.0) then

  ! Choice the way to impose the wall temperature (tmet)
  ! at the solid/fluid interface:
  !
  ! with the parameter itagms defined below:
  ! ----------------------------------------
  !  0 : A constant wall temperature imposed is given by the user
  !     ( tmet = tmet0 used as the wall temperature of the metal structures
  !      by the condensation model )
  !  1 : A variable wall temperature is imposed with a 0-D thermal model
  !     ( tmet = tmet(iel,1) computed by cs_metal_structures_tag.f90
  !      and used as the wall temperature of the metal structures
  !      by the condensation model )

  itagms = 1

  ! Wall temperature computed by a 0-D thermal model
  ! with an explicit scheme and variable over time.
  ! ------------------------------------------------
  ! Remark : the wall temperature unit is [°C].
  if(itagms.eq.1) then
    !----------------------------------------
    ! (xem) thickness of the metal structures
    !----------------------------------------
    xem = 0.024d0
    !-------------------------------------------
    ! Initial condition of the 0-D thermal model
    !-------------------------------------------
    tmet0 = 92.d0
    !--------------------------------------------
    ! Physical properties of the metal structures
    !--------------------------------------------
    ! (xro) density (kg.m-3)
    xro_m   = 8000.d0
    ! (xcond) thermal conductivity (W.m-1.C-1)
    xcond_m = 12.8d0
    ! (xcp)   Specific heat (J.kg-1.C-1)
    xcp_m   = 500.0d0
  else
    ! Wall temperature imposed as constant
    ! with a value specified by the user
    tmet = 92.d0
  endif

endif

Finally the source term type and values have to be set as follows:

!===============================================================================
! The user can specify here the values of the following arrays used by the
! modelling of the metal structures condensation:
!         - itypst(:,ivar) to specify the condensation source term type,
!         - svcond(:,ivar) the scalar value to multiply by the sink term array
!                          of the metal structures condensation model.
! These two arrays can be filled for each transported scalar.
!===============================================================================

!-- pointer to the specific heat
if (icp.ge.0) call field_get_val_s(icp, cpro_cp)

!-- pointer to the enthalpy value
ivarh = isca(iscalt)
call field_get_val_s(ivarfl(ivarh), cvar_h)

! loop over the cells associated to the metal structure
! source terms zone
do icmst = 1, ncmast
  iel = ltmast(icmst)

  ! Compute the enthalpy value of vapor gas
  if (ntcabs.le.1) then
    tk = t0
  else
    tk = cvar_h(iel)/cpro_cp(iel)
  endif
  hvap = s_h2o_g%cp*tk

  ! Condensation source terms associated
  ! to the metal structures imposed
  ! for each scalar.
  !----------------------------------------
  if (nscal.gt.0) then
    do iscal = 1, nscal
      if (iscal.eq.iscalt) then

        ! enthalpy value used for
        ! the explicit condensation term
        itypst(iel,isca(iscalt)) = 1
        svcond(iel,isca(iscalt)) = hvap
      else

        ! scalar values used for
        ! the explicit condensation term
        itypst(iel,isca(iscal)) = 1
        svcond(iel,isca(iscal)) = 0.d0
      endif
    enddo
  endif

enddo

Boundary source terms for condensation

The following variables need to be declared:

integer          ieltcd, ii, iz
integer          ifac, iel, iscal
integer          ivarh
integer          ilelt, nlelt
integer          izone
integer          f_id

double precision hvap, tk

type(gas_mix_species_prop) s_h2o_g

integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: cpro_cp
double precision, dimension(:), pointer :: cvar_h

Necessary species physical properties can be retrieved as follows:

! Allocate a temporary array for cells selection
allocate(lstelt(nfabor))

call field_get_id_try("y_h2o_g", f_id)
if (f_id.ne.-1) &
  call field_get_key_struct_gas_mix_species_prop(f_id, s_h2o_g)

The subroutine cs_user_boundary_mass_source_terms is called three times.

At the first call the number of boundary faces and number of zones on which a boundary mass source term is imposed is computed according to the selection criteria prescribed by the user.

if (iappel.eq.1.or.iappel.eq.2) then

!===============================================================================
! 1. One or two calls
! -------------------
!  - iappel = 1: nfbpcd: calculation of the number of faces with
!                             condensation source term
!  - iappel = 2: ifbpcd: index number of faces with condensation source terms
!
! Remarks
! =======
!  - Do not use spcond in this section (it is set on the third call, iappel=3)
!  - Do not use ifbpcd in this section on the first call (iappel=1)
!  - This section (iappel=1 or 2) is only accessed at the beginning of a
!     calculation. Should the localization of the condensation source terms evolve
!     in time, the user must identify at the beginning all cells that can
!     potentially become condensation  source term.
!===============================================================================


  !--To Select the cells with condensation source term
  !---------------------------------------------------

  izone = 0
  ieltcd = 0

  ! Cells with a boundary face of color 60
  call getfbr('60 and box[-0.03,0.0,0.5,0.03,0.6,1.22]',nlelt,lstelt)

  izone = izone + 1

  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    iel  = ifabor(ifac)
    ieltcd = ieltcd + 1
    if (iappel.eq.2) then
      ifbpcd(ieltcd) = ifac
      izzftcd(ieltcd) = izone
    endif
  enddo

  if (irangp.le.0.and.iappel.eq.2.and.ieltcd.gt.0) then
    write(*,*) "izzftcd(",izone,")= ", izzftcd(ieltcd)
  endif

  !========================================================
  ! Faces selection associated to the Concrete wall surface
  !========================================================
  call getfbr('60 and box[-0.03,0.0,1.22,0.03,0.6,2.55]', nlelt, lstelt)

  izone = izone + 1

  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    iel  = ifabor(ifac)
    ieltcd = ieltcd + 1
    if (iappel.eq.2) then
      ifbpcd(ieltcd) = ifac
      izzftcd(ieltcd) = izone
    endif
  enddo

  if (irangp.le.0.and.iappel.eq.2.and.ieltcd.gt.0) then
    write(*,*) "izzftcd(",izone,")= ", izzftcd(ieltcd)
  endif

endif

! For iappel = 1,
! Specification of nfbpcd.
if (iappel.eq.1) then
  nfbpcd = ieltcd
  nzones = izone
endif

The above part of the subroutine is also executed at the second call. In addition, at the second call, condensation models are chosen.

!===============================================================================
! Parameters padding of the wall thermal model and condensation model
! ------------------------------------------------------------------
! The condensation model can be used alone with a constant temperature
! specified by the user at the cold wall (at iappel=3 tpar=tpar0 in this case)
! or together with a 0-D thermal model. In the latter case, the two models are
! coupled.
!===============================================================================

if (iappel.eq.2) then

  do ii = 1, nfbpcd
    iz = izzftcd(ii)
    if (iz.eq.1.and.icondb.eq.0) then
      ! Turbulent law and empiric correlations used to
      ! define the exchange coefficients of the sink
      ! source term and heat transfer to the cooling
      ! wall associated to the condensation phenomenon
      ! given by the COPAIN model.
      !-----------------------------------------------

      ! Choice the way to compute the exchange coefficient (hcond)
      ! associated to the condensation sink source term.

      ! With the parameter icophc defined below:
      ! ----------------------------------------
      !  1 : this one provided by the turbulent flow
      !        at the cooling wall (hcond = hcdt)
      !  2 : this one is given by the copain
      !        correlation (hcond = hcdcop)
      !  3 : this one is obtained by the estimation of
      !        the maximal value between those two previous
      !        coefficient as below : hcond= max(hcdt, hcdcop)
      izcophc(iz) = 3

      ! Choice the way to compute the thermal exchange coefficient
      ! associated to the heat transfer at the cooling wall,
      ! due to the energy loss by condensation phenomenon.

      ! With the parameter icophg defined below:
      ! ----------------------------------------
      !  2 : this one is given by the copain
      !      correlation (hpcond = hw_cop)
      !  3 : this one is obtained by the estimation of
      !      the maximal value between the current
      !      and previous value of hcdcop given by
      !      the copain correlation as below:
      !         hpcond= max(hw_cop^n, hw_cop^n+1)
      izcophg(iz) = 3

      ! Choice the way to impose the wall temperature (tpar)
      ! at the solid/fluid interface:
      !
      ! with the parameter itag1d defined below:
      ! ----------------------------------------
      !  0 : A constant wall temperature imposed is given by the user
      !     ( tpar = tpar0 used as the wall temperature by the condensation model )
      !  1 : A variable wall temperature is imposed with a 1-D thermal model
      !     ( tpar = tmur(ii,1) computed by tagmro.f90 and used as the
      !       wall temperature by the condensation model )
      iztag1d(iz) = 1

      ! Wall temperature computed by a 1-D thermal model
      ! with a implicit scheme and variable over time.
      ! ------------------------------------------------
      ! Remark : the wall temperature is in unit [°C].
      if (iztag1d(iz).eq.1) then
        !---------------------------------------------------
        ! Numerical parameters used by the 1-D thermal model
        !---------------------------------------------------
        ! (theta) parameter of the implicit scheme
        ztheta(iz) = 1.d0
        ! (dxmin) First cell size of the 1D mesh
        ! -> with (dxmin.eq.0) the Dx is constant
        ! -> with (dxmin.ne.0) the Dx is variable
        zdxmin(iz) = 0.d0
        ! (nmur) space steps number of the 1D mesh
        znmur(iz) = 10
        ! (epais) thickness of the 1D wall
        zepais(iz) = 0.024d0

        !-------------------------------------------
        !Initial condition of the 1-D thermal model
        !-------------------------------------------
        ztpar0(iz) = 26.57d0
      endif
    endif

  enddo

  if (irangp.ge.0) then
    call parimx(nzones, izcophc)
    call parimx(nzones, izcophg)
    call parimx(nzones, iztag1d)
    call parrmx(nzones, ztheta)
    call parrmx(nzones, zdxmin)
    call parimx(nzones, znmur )
    call parrmx(nzones, zepais)
    call parrmx(nzones, ztpar0)
  endif

Finally at the third call, the source term types and values have to be set.

elseif (iappel.eq.3) then

!===============================================================================
! 2. For nfbpcd > 0 , third call
!    iappel = 3 : itypcd: type of condensation source term
!                  spcond: condensation source term
! Remark
! ======
! If itypcd(ieltcd,ivar) is set to 1, spcond(ieltcd,ivar) must be set.
!===============================================================================

  !-- pointer to the specific heat
  if (icp.ge.0) call field_get_val_s(icp, cpro_cp)

  !-- pointer to the enthalpy value
  ivarh = isca(iscalt)
  call field_get_val_s(ivarfl(ivarh), cvar_h)

  do ii = 1, nfbpcd

    ifac = ifbpcd(ii)
    iel  = ifabor(ifac)

    iz = izzftcd(ii)

    if (icondb.eq.0) then
      if (iztag1d(iz).eq.1) then
        !-------------------------------------------
        !Boundary conditions of the 1-D thermal model
        !-------------------------------------------
        zhext(iz) =  1.d+8 ; ztext(iz) = 26.57d0
        ! --------------------------------------------
        ! Physical properties of the concrete material
        ! --------------------------------------------
        ! (rob) density (kg.m-3)
        zrob(iz)   = 8000.d0
        ! (condb) thermal conductivity (W.m-1.C-1)
        zcondb(iz) = 12.8d0
        ! (cpb)   Specific heat (J.kg-1.C-1)
        zcpb(iz) = 500.0d0
      else
        ! Wall temperature imposed as constant
        ! with a value specified by the user
        ztpar(iz) = 26.57d0
      endif
    elseif (iz.eq.2.and.icondb.eq.0) then
      if (iztag1d(iz).eq.1) then
        !-------------------------------------------
        !Boundary conditions of the 1-D thermal model
        !-------------------------------------------
        zhext(iz) =  1.d+8 ; ztext(iz) = 26.57d0
        ! --------------------------------------------
        ! Physical properties of the concrete material
        ! --------------------------------------------
        ! (rob) density (kg.m-3)
        zrob(iz)   = 8000.d0
        ! (condb) thermal conductivity (W.m-1.C-1)
        zcondb(iz) = 12.8d0
        ! (cpb)   Specific heat (J.kg-1.C-1)
        zcpb(iz) = 500.0d0
      else
        ! Wall temperature imposed as constant
        ! with a value specified by the user
        ztpar(iz) = 26.57d0
      endif
    endif

    ! To fill the spcond(nfbpcd,ivar) array
    ! if we want to specify a variable value
    !---------------------------------------

    ! Compute the enthalpy value of vapor gas
    if (ntcabs.le.1) then
      tk = t0
    else
      tk = cvar_h(iel)/cpro_cp(iel)
    endif
    hvap = s_h2o_g%cp*tk

    ! any condensation source term
    ! associated to each velocity component
    ! momentum equation in this case.
    !----------------------------------------
    itypcd(ii,iu) = 0
    spcond(ii,iu) = 0.d0
    itypcd(ii,iv) = 0
    spcond(ii,iv) = 0.d0
    itypcd(ii,iw) = 0
    spcond(ii,iw) = 0.d0

    ! any condensation source term
    ! associated to each turbulent variables
    ! for (k -eps) standrad turbulence model
    !----------------------------------------
    if (itytur.eq.2) then
      itypcd(ii,ik ) = 0
      spcond(ii,ik ) = 0.d0
      itypcd(ii,iep) = 0
      spcond(ii,iep) = 0.d0
    endif
    if (nscal.gt.0) then
      do iscal = 1, nscal
        if (iscal.eq.iscalt) then

          ! enthalpy value used for
          ! the explicit condensation term
          itypcd(ii,isca(iscalt)) = 1
          spcond(ii,isca(iscalt)) = hvap
        else

          ! scalar values used for
          ! the explicit condensation term
          itypcd(ii,isca(iscal)) = 1
          spcond(ii,isca(iscal)) = 0.d0
        endif
      enddo
    endif

  enddo

  if (irangp.ge.0) then
    call parrmx(nzones, zhext)
    call parrmx(nzones, ztext)
    call parrmx(nzones, zrob )
    call parrmx(nzones, zcondb)
    call parrmx(nzones, zcpb )
    call parrmx(nzones, ztpar)
  endif

endif

Boundary source terms for condensation

The following variables need to be declared:

integer          ieltcd
integer          ifac, iel, iscal
integer          ivarh
integer          ilelt, nlelt
integer          izone
integer          f_id

double precision hvap, tk

type(gas_mix_species_prop) s_h2o_g

integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: cpro_cp
double precision, dimension(:), pointer :: cvar_h

Necessary species physical properties can be retrieved as follows:

! Allocate a temporary array for cells selection
allocate(lstelt(nfabor))

call field_get_id_try("y_h2o_g", f_id)
if (f_id.ne.-1) &
  call field_get_key_struct_gas_mix_species_prop(f_id, s_h2o_g)

The subroutine cs_user_boundary_mass_source_terms is called three times.

At the first call the number of boundary faces and number of zones on which a boundary mass source term is imposed is computed according to the selection criteria prescribed by the user.

if (iappel.eq.1.or.iappel.eq.2) then

!===============================================================================
! 1. One or two calls
! -------------------
!  - iappel = 1: nfbpcd: calculation of the number of faces with
!                             condensation source term
!  - iappel = 2: ifbpcd: index number of faces with condensation source terms
!
! Remarks
! =======
!  - Do not use spcond in this section (it is set on the third call, iappel=3)
!  - Do not use ifbpcd in this section on the first call (iappel=1)
!  - This section (iappel=1 or 2) is only accessed at the beginning of a
!     calculation. Should the localization of the condensation source terms evolve
!     in time, the user must identify at the beginning all cells that can
!     potentially becomea condensation  source term.
!===============================================================================


  !--To Select the cells with condensation source term
  !---------------------------------------------------

  izone = 0
  ieltcd = 0

  ! Cells with a boundary face of color 60
  call getfbr('60',nlelt,lstelt)

  izone = izone + 1

  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    iel  = ifabor(ifac)
    izftcd(iel) = izone
    ieltcd = ieltcd + 1
    if (iappel.eq.2) ifbpcd(ieltcd) = ifac
  enddo

endif

! For iappel = 1,
! Specification of nfbpcd.
if (iappel.eq.1) then
  nfbpcd = ieltcd
endif

The above part of the subroutine is also executed at the second call. In addition, at the second call, condensation models are chosen.

!===============================================================================
! Parameters padding of the 1-D thermal model and condensation model
! ------------------------------------------------------------------
! Both models can be enabled and coupled together or
! the condensation model can be used without enabling the 1-D thermal model.
! In this latter case a constant wall temperature must be specified by the
! user at the cold wall (at the third call (iappel=3) tpar=tpar0).
!===============================================================================
if (iappel.eq.2) then

  if (icondb.eq.0) then

    ! Turbulent law and empiric correlations used to
    ! define the exchange coefficients of the sink
    ! source term and heat transfer to the cooling
    ! wall associated to the condensation phenomenon
    ! given by the COPAIN model.
    !-----------------------------------------------

    ! Choice the way to compute the exchange coefficient (hcond)
    ! associated to the condensation sink source term.

    ! With the parameter icophc defined below:
    ! ----------------------------------------
    !  1 : this one provided by the turbulent flow
    !        at the cooling wall (hcond = hcdt)
    !  2 : this one is given by the copain
    !        correlation (hcond = hcdcop)
    !  3 : this one is obtained by the estimation of
    !        the maximal value between those two previous
    !        coefficient as below : hcond= max(hcdt, hcdcop)

    icophc = 3

    ! Choice the way to compute the thermal exchange coefficient
    ! associated to the heat transfer at the cooling wall,
    ! due to the energy loss by condensation phenomenon.

    ! With the parameter icophg defined below:
    ! ----------------------------------------
    !  2 : this one is given by the copain
    !      correlation (hpcond = hw_cop)
    !  3 : this one is obtained by the estimation of
    !      the maximal value between the current
    !      and previous value of hcdcop given by
    !      the copain correlation as below:
    !         hpcond= max(hw_cop^n, hw_cop^n+1)

    icophg = 3

    ! Choice the way to impose the wall temperature (tpar)
    ! at the solid/fluid interface:
    !
    ! with the parameter itag1d defined below:
    ! ----------------------------------------
    !  0 : A constant wall temperature imposed is given by the user
    !     ( tpar = tpar0 used as the wall temperature by the condensation model )
    !  1 : A variable wall temperature is imposed with a 1-D thermal model
    !     ( tpar = tmur(ii,1) computed by tagmro.f90 and used as the
    !       wall temperature by the condensation model )

    itag1d = 1

    ! Wall temperature computed by a 1-D thermal model
    ! with a implicit scheme and variable over time.
    ! ------------------------------------------------
    ! Remark : the wall temperature is in unit [°C].
    if(itag1d.eq.1) then

      !---------------------------------------------------
      ! Numerical parameters used by the 1-D thermal model
      !---------------------------------------------------

      ! (theta) parameter of the implicit scheme
      theta = 1.d0
      ! (dxmin) First cell size of the 1D mesh
      ! -> with (dxmin.eq.0) the Dx is constant
      ! -> with (dxmin.ne.0) the Dx is variable
      dxmin = 0.d0
      ! (nmur) space steps number of the 1D mesh
      nmur = 10
      ! (epais) thickness of the 1D wall
      epais = 0.024d0

      !-------------------------------------------
      !Initial condition of the 1-D thermal model
      !-------------------------------------------
      tpar0 = 26.57d0

   endif

  endif

Finally at the third call, the source term types and values have to be set.

elseif (iappel.eq.3) then

!===============================================================================
! 2. For nfbpcd > 0 , third call
!    iappel = 3 : itypcd: type of condensation source term
!                  spcond: condensation source term
! Remark
! ======
! If itypcd(ieltcd,ivar) is set to 1, spcond(ieltcd,ivar) must be set.
!===============================================================================

  !-- pointer to the specific heat
  if (icp.ge.0) call field_get_val_s(icp, cpro_cp)

  !-- pointer to the enthalpy value
  ivarh = isca(iscalt)
  call field_get_val_s(ivarfl(ivarh), cvar_h)

  if (icondb.eq.0) then
    if(itag1d.eq.1) then
      !-------------------------------------------
      !Boundary conditions of the 1-D thermal model
      !-------------------------------------------
      hext =  1.d+8 ; text = 26.57d0
      ! --------------------------------------------
      ! Physical properties of the concrete material
      ! --------------------------------------------
      ! (rob) density (kg.m-3)
      rob   = 8000.d0
      ! (condb) thermal conductivity (W.m-1.C-1)
      condb = 12.8d0
      ! (cpb)   Specific heat (J.kg-1.C-1)
      cpb = 500.0d0

    else
      ! Wall temperature imposed as constant
      ! with a value specified by the user
      tpar = 26.57d0
    endif
  endif

  ! To fill the spcond(nfbpcd,ivar) array
  ! if we want to specify a variable value
  !---------------------------------------
  do ieltcd = 1, nfbpcd

    ifac = ifbpcd(ieltcd)
    iel  = ifabor(ifac)

    ! Compute the enthalpy value of vapor gas
    if (ntcabs.le.1) then
      tk = t0
    else
      tk = cvar_h(iel)/cpro_cp(iel)
    endif
    hvap = s_h2o_g%cp*tk

    ! any condensation source term
    ! associated to each velocity component
    ! momentum equation in this case.
    !----------------------------------------
    itypcd(ieltcd,iu) = 0
    spcond(ieltcd,iu) = 0.d0
    itypcd(ieltcd,iv) = 0
    spcond(ieltcd,iv) = 0.d0
    itypcd(ieltcd,iw) = 0
    spcond(ieltcd,iw) = 0.d0

    ! any condensation source term
    ! associated to each turbulent variables
    ! for (k -eps) standrad turbulence model
    !----------------------------------------
    if (itytur.eq.2) then
      itypcd(ieltcd,ik ) = 0
      spcond(ieltcd,ik ) = 0.d0
      itypcd(ieltcd,iep) = 0
      spcond(ieltcd,iep) = 0.d0
    endif
    if (nscal.gt.0) then
      do iscal = 1, nscal
        if (iscal.eq.iscalt) then

          ! enthalpy value used for
          ! the explicit condensation term
          itypcd(ieltcd,isca(iscalt)) = 1
          spcond(ieltcd,isca(iscalt)) = hvap
        else

          ! scalar values used for
          ! the explicit condensation term
          itypcd(ieltcd,isca(iscal)) = 1
          spcond(ieltcd,isca(iscal)) = 0.d0
        endif
      enddo
    endif


  enddo

endif