set_ic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SET_IC                                                  !
!  Author: M. Syamlal                                 Date: 21-JAN-92  !
!                                                                      !
!  Purpose: This module sets all the initial conditions.               !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE set_ic

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE constant
      USE physprop
      USE ic
      USE fldvar
      USE visc_g
      USE indices
      USE scales
      USE energy
      USE scalars
      USE compar
      USE run
      USE sendrecv
      USE solids_pressure
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! indices
      INTEGER :: I, J, K, IJK, M
! Local index for initial condition
      INTEGER :: L
! Temporary variable for storing IC_EP_g
      DOUBLE PRECISION :: EPGX
! Temporary variable for storing IC_P_g
      DOUBLE PRECISION :: PGX
! Temporary variable for storing P_star
      DOUBLE PRECISION :: PSX
! Temporary variable for storing IC_T_g
      DOUBLE PRECISION :: TGX
! Temporary variable for storing IC_U_g
      DOUBLE PRECISION :: UGX
! Temporary variable for storing IC_V_g
      DOUBLE PRECISION :: VGX
! Temporary variable for storing IC_W_g
      DOUBLE PRECISION :: WGX
! Temporary variable for storing IC_ROP_s
      DOUBLE PRECISION :: ROPSX (dimension_m)
! Temporary variable for storing IC_T_s
      DOUBLE PRECISION :: TSX (dimension_m)
! Temporary variable for storing IC_U_s
      DOUBLE PRECISION :: USX (dimension_m)
! Temporary variable for storing IC_V_s
      DOUBLE PRECISION :: VSX (dimension_m)
! Temporary variable for storing IC_W_s
      DOUBLE PRECISION :: WSX (dimension_m)
! number density for GHD theory
      DOUBLE PRECISION :: nM, nTOT
!-----------------------------------------------

!  Set the initial conditions.
      DO l = 1, dimension_ic
         IF (ic_defined(l)) THEN

            epgx = ic_ep_g(l)
            pgx = ic_p_g(l)
            psx = ic_p_star(l)
            IF (psx==undefined .AND. ic_type(l)/='PATCH') psx = zero
            tgx = ic_t_g(l)
            ugx = ic_u_g(l)
            vgx = ic_v_g(l)
            wgx = ic_w_g(l)

            m = 1
            IF (mmax > 0) THEN
              ropsx(:mmax) = ic_rop_s(l,:mmax)
              tsx(:mmax) = ic_t_s(l,:mmax)
              usx(:mmax) = ic_u_s(l,:mmax)
              vsx(:mmax) = ic_v_s(l,:mmax)
              wsx(:mmax) = ic_w_s(l,:mmax)
              m = mmax + 1
            ENDIF

            DO k = ic_k_b(l), ic_k_t(l)
            DO j = ic_j_s(l), ic_j_n(l)
            DO i = ic_i_w(l), ic_i_e(l)
               IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
               IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
               ijk = funijk(i,j,k)

               IF (.NOT.wall_at(ijk)) THEN
                  IF (epgx /= undefined) ep_g(ijk) = epgx

                  IF (ic_type(l) == 'PATCH') THEN
                      IF (pgx /= undefined) p_g(ijk) = scale_pressure(pgx)
                  ELSE
                     p_g(ijk) = merge(scale_pressure(pgx), undefined,           &
                        pgx /= undefined)
                  ENDIF

                  IF (psx /= undefined) p_star(ijk) = psx
                  IF (tgx /= undefined) t_g(ijk) = tgx
                  IF (ic_l_scale(l) /= undefined) l_scale(ijk) =       &
                      ic_l_scale(l)

                  IF (nmax(0) > 0) THEN
                     WHERE (ic_x_g(l,:nmax(0)) /= undefined)           &
                        x_g(ijk,:nmax(0)) = ic_x_g(l,:nmax(0))
                  ENDIF

                  IF (nscalar > 0) THEN
                     WHERE (ic_scalar(l,:nscalar) /= undefined)        &
                        scalar(ijk,:nscalar) = ic_scalar(l,:nscalar)
                  ENDIF

                  IF (k_epsilon) THEN
                      IF (ic_k_turb_g(l) /= undefined)                 &
                          k_turb_g(ijk) = ic_k_turb_g(l)
                      IF (ic_e_turb_g(l) /= undefined)                 &
                          e_turb_g(ijk) = ic_e_turb_g(l)
                  ENDIF

                  IF (ugx /= undefined) u_g(ijk) = ugx
                  IF (vgx /= undefined) v_g(ijk) = vgx
                  IF (wgx /= undefined) w_g(ijk) = wgx

                  gama_rg(ijk) = ic_gama_rg(l)
                  t_rg(ijk) = merge(ic_t_rg(l), zero,                  &
                     ic_t_rg(l) /= undefined)

                  DO m = 1, mmax
                    IF (ropsx(m) /= undefined) rop_s(ijk,m) = ropsx(m)
                    IF (tsx(m) /= undefined) t_s(ijk,m) = tsx(m)
                    IF (ic_theta_m(l,m) /= undefined)                  &
                       theta_m(ijk,m) = ic_theta_m(l,m)
                    IF (usx(m) /= undefined) u_s(ijk,m) = usx(m)
                    IF (vsx(m) /= undefined) v_s(ijk,m) = vsx(m)
                    IF (wsx(m) /= undefined) w_s(ijk,m) = wsx(m)

                    gama_rs(ijk,m) = ic_gama_rs(l,m)
                    t_rs(ijk,m) = merge(ic_t_rs(l,m),zero,             &
                       ic_t_rs(l,m) /= undefined)

                    IF (nmax(m) > 0) THEN
                        WHERE (ic_x_s(l,m,:nmax(m)) /= undefined)      &
                           x_s(ijk,m,:nmax(m)) = ic_x_s(l,m,:nmax(m))
                    ENDIF
                  ENDDO

! for GHD theory to compute mixture IC of velocity and density
                  IF(kt_type_enum == ghd_2007) THEN
                     rop_s(ijk,mmax) = zero
                     u_s(ijk,mmax) = zero
                     v_s(ijk,mmax) = zero
                     w_s(ijk,mmax) = zero
                     theta_m(ijk,mmax) = zero
                     ntot = zero
                     nm = zero
                     DO m = 1, smax
                        IF (ropsx(m) /= undefined) THEN
                           rop_s(ijk,mmax) = rop_s(ijk,mmax) + ropsx(m)
                           nm = ropsx(m)*6d0 /                         &
                              (pi*d_p(ijk,m)**3*ro_s(ijk,m))
                           ntot = ntot + nm
                        ENDIF
                        IF (ic_theta_m(l,m) /= undefined)              &
                           theta_m(ijk,mmax) = theta_m(ijk,mmax) +     &

                           nm*ic_theta_m(l,m)
                        IF(usx(m) /= undefined .AND.                   &
                           ropsx(m) /= undefined)  u_s(ijk,mmax) =     &
                           u_s(ijk,mmax) + ropsx(m)*usx(m)

                        IF(vsx(m) /= undefined .AND.                   &
                           ropsx(m) /= undefined) v_s(ijk,mmax) =      &
                           v_s(ijk,mmax) +  ropsx(m)*vsx(m)

                        IF(wsx(m) /= undefined .AND.                   &
                           ropsx(m) /= undefined) w_s(ijk,mmax) =      &
                           w_s(ijk,mmax) +  ropsx(m)*wsx(m)
                     ENDDO

! If ropsTotal > 0 then RoN_T > 0
                     IF(rop_s(ijk,mmax) > zero) THEN
                        u_s(ijk,mmax) = u_s(ijk,mmax) / rop_s(ijk,mmax)
                        v_s(ijk,mmax) = v_s(ijk,mmax) / rop_s(ijk,mmax)
                        w_s(ijk,mmax) = w_s(ijk,mmax) / rop_s(ijk,mmax)
                        theta_m(ijk,mmax) = theta_m(ijk,mmax) / ntot
! For initially empty bed:
                     ELSE
                        u_s(ijk,mmax) = u_s(ijk,mmax)
                        v_s(ijk,mmax) = v_s(ijk,mmax)
                        w_s(ijk,mmax) = w_s(ijk,mmax)
! Set T > 0 in case Ti > 0
                        DO m = 1, smax
                           theta_m(ijk,mmax) = theta_m(ijk,m)
                        ENDDO
                        IF(theta_m(ijk,mmax)==zero)                     &
                           theta_m(ijk,mmax) = small_number
                    ENDIF
                  ENDIF
! end of modifications for GHD theory
               ENDIF     ! Fluid at
            ENDDO   ! over i
            ENDDO   ! over j
            ENDDO   ! over k
         ENDIF   ! if (ic_defined)
      ENDDO   ! over dimension_ic

      CALL send_recv(l_scale,2)

      RETURN
      END SUBROUTINE set_ic