#include "cppdefs.h"
      MODULE ice_frazil_mod
#if (defined ICE_MODEL && defined ICE_THERMO && !defined TS_FIXED) \
      || defined CICE_MODEL
!
!=======================================================================
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!================================================== Hernan G. Arango ===
!                                                                      !
!  This routine computes the frazil ice growth in the water when the
!  water temperature gets below freezing. It adjusts the water
!  temperature and salinity accordingly.
!
!  Reference: Steele et al. (1989). JPO, 19, 139-147.
!                                                                      !
!=======================================================================
!
      implicit none

      PRIVATE
      PUBLIC ice_frazil

      CONTAINS

      SUBROUTINE ice_frazil (ng, tile)

      USE mod_param
      USE mod_grid
      USE mod_ocean
      USE mod_ice
      USE mod_stepping
# ifdef CICE_MODEL
      USE mod_forces
# endif

      integer, intent(in) :: ng, tile
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, iNLM, 80, __LINE__, __FILE__)
# endif
!
      CALL ice_frazil_tile (ng, tile,                                   &
     &                      LBi, UBi, LBj, UBj,                         &
     &                      IminS, ImaxS, JminS, JmaxS,                 &
     &                      nnew(ng),                                   &
# ifdef MASKING
     &                      GRID(ng) % rmask,                           &
# endif
# ifdef WET_DRY
     &                      GRID(ng) % rmask_wet,                       &
# endif
     &                      GRID(ng) % Hz,                              &
     &                      GRID(ng) % z_r,                             &
     &                      OCEAN(ng) % rho,                            &
     &                      OCEAN(ng) % t,                              &
# ifdef CICE_MODEL
     &                      FORCES(ng) % sustr,                         &
     &                      FORCES(ng) % svstr,                         &
     &                      FORCES(ng) % stflx,                         &
     &                      ICE(ng) % ai,                               &
     &                      ICE(ng) % hi,                               &
     &                      GRID(ng) % z_w,                             &
# endif
     &                      ICE(ng) % wfr)
# ifdef PROFILE
      CALL wclock_off (ng, iNLM, 80, __LINE__, __FILE__)
# endif
      RETURN
      END SUBROUTINE ice_frazil
!
!***********************************************************************
      subroutine ice_frazil_tile (ng, tile,                             &
     &                            LBi, UBi, LBj, UBj,                   &
     &                            IminS, ImaxS, JminS, JmaxS,           &
     &                            nnew,                                 &
# ifdef MASKING
     &                            rmask,                                &
# endif
# ifdef WET_DRY
     &                            rmask_wet,                            &
# endif
     &                            Hz, z_r, rho, t,                      &
# ifdef CICE_MODEL
     &                            sustr, svstr, stflx, ai, hi, z_w,     &
# endif
     &                            wfr)
!***********************************************************************
!
      USE mod_param
      USE mod_scalars
# ifdef CICE_MODEL
      USE shr_const_mod
# endif
!
      USE bc_2d_mod, ONLY : bc_r2d_tile
      USE exchange_3d_mod, ONLY : exchange_r3d_tile
# ifdef DISTRIBUTE
      USE mp_exchange_mod
      USE distribute_mod, ONLY : mp_reduce
# endif
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: nnew

# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: rmask(LBi:,LBj:)
#  endif
#  ifdef WET_DRY
      real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
#  endif
      real(r8), intent(in) :: Hz(LBi:,LBj:,:)
      real(r8), intent(in) :: z_r(LBi:,LBj:,:)
      real(r8), intent(in) :: rho(LBi:,LBj:,:)
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
#  ifdef CICE_MODEL
      real(r8), intent(in) :: sustr(LBi:,LBj:)
      real(r8), intent(in) :: svstr(LBi:,LBj:)
      real(r8), intent(in) :: stflx(LBi:,LBj:,:)
      real(r8), intent(in) :: ai(LBi:,LBj:)
      real(r8), intent(in) :: hi(LBi:,LBj:)
      real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
#  endif
      real(r8), intent(out) :: wfr(LBi:,LBj:)
# else
#  ifdef MASKING
      real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
#  endif
#  ifdef WET_DRY
      real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
#  ifdef CICE_MODEL
      real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: stflx(LBi:UBi,LBj:UBj,NT(ng))
      real(r8), intent(in) :: ai(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: hi(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
#  endif
      real(r8), intent(out) :: wfr(LBi:UBi,LBj:UBj)
# endif
!
! Local variable definitions
!
      integer :: i, j, k, itrc

      real(r8), parameter :: Lhat = 79.2_r8
      real(r8), parameter :: r = 0.5_r8

      real(r8) :: t_freeze
      real(r8) :: s1
      real(r8) :: z1
      real(r8) :: sice
      real(r8) :: gamma_k
      real(r8) :: t_fr
      real(r8) :: ice_dens
      real(r8) :: delta_wfr
# ifdef CICE_MODEL
      real(r8), parameter :: depressT = -0.054_r8
# else
      real(r8), parameter :: depressT = -0.0543_r8
# endif

# ifdef CICE_MODEL
#  ifdef ICE_LOG_LAYER
! Compute heat flux from ocean a la Mellor and Kantha, 1989
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: t0mk
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: s0mk
      real(r8), parameter :: prt = 13._r8
      real(r8), parameter :: prs = 2432._r8
      real(r8), parameter :: tpr = 0.85_r8
      real(r8), parameter :: nu = 1.8E-6_r8
      real(r8), parameter :: z0ii = 0.02_r8
      real(r8), parameter :: kappa = 0.4_r8
      real(r8), parameter :: hfus = 3.347E+5_r8         ! [J kg-1]
      real(r8), parameter :: ykf = 3.14
      real(r8), parameter :: eps = 1.0e-20_r8
      real(r8) :: utau, dztop, rno, termt, terms, cht, chs, z0, zdz0
#  endif
      real(r8) :: dz, heat, tfr
      real(r8), parameter :: bl_thick = 5.0_r8
      real(r8) :: ice_thick, thickness, fract
      real(r8), parameter :: latent_heat_fusion = 3.337e5_r8   ! J/kg
      real(r8), parameter :: cp_sw              = 3996.0_r8    ! J/kg/K
      real(r8), parameter :: rho_sw             = 1026.0_r8    ! kg/m^3
      real(r8), parameter :: rho_fw             = 1000.0_r8    ! kg/m^3
      real(r8), parameter :: cp_over_lhfusion   =                       &
     &                 rho_sw*cp_sw/(latent_heat_fusion*rho_fw)
      real(r8), parameter :: rhocpr = 4093740._r8       ! [W s m-3 K-1]
      real(r8), parameter :: salice = 4.0_r8
      real(r8), parameter :: salref = 34.7_r8
      real(r8) :: potice
# endif

# ifdef DISTRIBUTE
      real(r8), allocatable           :: buffer(:)
      character (len=3), allocatable  :: op_handle(:)
# endif

!  Inline functions
!  Freezing temperature (Gill, 1982)
!     t_freeze(s1,z1) = -0.0575*s1 + 1.710523d-3*sqrt(s1)**3
!    &       - 2.154996d-4*s1*s1 + 0.000753*z1
!  Freezing temperature (Steele et al. 1989)
!      t_freeze(s1,z1) = -0.0543*s1 + 0.000759*z1
# ifdef CICE_MODEL
! Need to match CICE
      t_freeze(s1,z1) = -1.8_r8
!      t_freeze(s1,z1) = depressT*s1
# else
! Temperature or potential temperature??
      t_freeze(s1,z1) = depressT*s1
# endif

# include "set_bounds.h"

# ifdef DISTRIBUTE
      IF (.not. allocated(buffer))    allocate (buffer(1))
      IF (.not. allocated(op_handle)) allocate (op_handle(1))
# endif

# ifdef CICE_MODEL
      ice_dens = SHR_CONST_RHOICE
# else
      ice_dens = rhoice(ng)
# endif
      DO j=Jstr,Jend
        DO i=Istr,Iend
          wfr(i,j) = 0.0_r8
        ENDDO
      ENDDO
# if defined CICE_MODEL
      DO j=Jstr,Jend
        DO i=Istr,Iend
#  ifdef WET_DRY
          IF (rmask_wet(i,j) .ne. 0.0_r8) THEN
#  elif defined MASKING
          IF (rmask(i,j) .ne. 0.0_r8) THEN
#  endif
            thickness = 0.0_r8
            DO k=N(ng),1,-1
              dz = Hz(i,j,k)
              t_fr = t_freeze(t(i,j,k,nnew,isalt),z_r(i,j,k))
              IF ((thickness + dz) > bl_thick) THEN
                fract = (bl_thick - thickness)/dz
                dz = fract*dz
              END IF
              thickness = thickness + dz
              potice = (t_fr - t(i,j,k,nnew,itemp)) * dz
!
!*** if potice < 0, compute again later
!*** if potice > 0, keep on freezing (wfr > 0)
!
              IF ((wfr(i,j) + potice) > 0) THEN
                wfr(i,j) = wfr(i,j) + potice
              ELSE
                potice = -wfr(i,j)
                wfr(i,j) = 0.0_r8
              ENDIF
              t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) +             &
     &                potice/Hz(i,j,k)
              t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt)               &
     &              + (salref-salice)*potice*cp_over_lhfusion/Hz(i,j,k)
!*** jump out if done with surface layer
              IF (z_w(i,j,N(ng))-z_w(i,j,k-1) > bl_thick) EXIT
            END DO
#  if defined WET_DRY || defined MASKING
          END IF
#  endif
        END DO
      END DO
# else
! Old ice model
      DO j=Jstr,Jend
        DO i=Istr,Iend
          DO k=1,N(ng)
#  ifdef WET_DRY
            IF (rmask_wet(i,j) .ne. 0.0_r8) THEN
#  elif defined MASKING
            IF (rmask(i,j) .ne. 0.0_r8) THEN
#  endif
              t_fr = t_freeze(t(i,j,k,nnew,isalt),z_r(i,j,k))
              IF (t(i,j,k,nnew,itemp) .lt. t_fr) THEN
                gamma_k = (t_fr - t(i,j,k,nnew,itemp)) /                &
     &                     (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r)     &
     &                         - depressT * t(i,j,k,nnew,isalt))
                IF (gamma_k .lt. 0.0_r8) THEN
                  print *, 'trouble in ice_frazil', i, j, k,            &
     &             t(i,j,k,nnew,itemp), t(i,j,k,nnew,isalt),            &
     &             t_fr, wfr(i,j), gamma_k, Hz(i,j,k)
                  exit_flag = 10
                END IF
                wfr(i,j) = wfr(i,j) + gamma_k * Hz(i,j,k) *             &
     &                    (rho0 + rho(i,j,k) ) / ice_dens
                t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + gamma_k *   &
     &                 (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r))
                sice = MIN(sice_ref,t(i,j,k,nnew,isalt))
                t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) *             &
     &                  (1.0_r8 + gamma_k) - sice*gamma_k
              ELSE IF (wfr(i,j) > 0 .and.                               &
     &                 t(i,j,k,nnew,itemp) .gt. t_fr) THEN
! Use heat at this level to melt some ice from below.
! gamma_k becomes negative here.
                gamma_k = (t_fr - t(i,j,k,nnew,itemp)) /                &
     &                     (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r)     &
     &                         - depressT * t(i,j,k,nnew,isalt))
                delta_wfr = gamma_k * Hz(i,j,k) *                       &
     &                    (rho0 + rho(i,j,k) ) / rhoice(ng)
                IF ((wfr(i,j) + delta_wfr) > 0) THEN
                  wfr(i,j) = wfr(i,j) + delta_wfr
                ELSE
                  gamma_k = -wfr(i,j) * rhoice(ng) /                    &
     &                     (Hz(i,j,k)*(rho0+rho(i,j,k)))
                  wfr(i,j) = 0.0_r8
                ENDIF
                t(i,j,k,nnew,itemp) = t(i,j,k,nnew,itemp) + gamma_k *   &
     &                 (Lhat + t(i,j,k,nnew,itemp)*(1.0_r8 - r))
                t(i,j,k,nnew,isalt) = t(i,j,k,nnew,isalt) *             &
     &                                  (1.0_r8 + gamma_k)
              END IF
#  if defined WET_DRY || defined MASKING
            END IF
#  endif
          END DO
          wfr(i,j) = wfr(i,j)/dt(ng)
          IF (wfr(i,j) .lt. 0.0_r8) THEN
            print *, 'trouble in ice_frazil', i, j,                     &
     &         t(i,j,N(ng),nnew,itemp), t(i,j,N(ng),nnew,isalt),        &
     &         wfr(i,j), gamma_k, Hz(i,j,N(ng))
            exit_flag = 10
          END IF
        END DO
      END DO
# endif

# if defined CICE_MODEL
! Now compute heat flux from ocean and subtract from wfr.
! Using old salt flux, solve for surface salinity s0mk.
#  ifdef ICE_LOG_LAYER
      DO j=Jstr,Jend
        DO i=Istr,Iend
          dztop=z_w(i,j,N(ng))-z_r(i,j,N(ng))
          ice_thick = 0.05_r8+hi(i,j)/MAX(ai(i,j),eps)
          utau = sqrt(sqrt(                                             &
     &             (0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2                &
     &           + (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2                &
     &                   )    )
          utau = max(utau,1.E-4_r8)
!  Need some roughness estimate here
          z0 = max(z0ii*ice_thick,0.01_r8)
          z0 = min(z0,0.1_r8)
!
!     *** Yaglom and Kader formulation for z0t and z0s
!
          zdz0 = dztop/z0   !WPB
          zdz0 = MAX(zdz0,3._r8)

!  Was this:
          rno = utau*0.09_r8/nu
! Mellor
!          rno = utau*z0/nu
          termt = ykf*sqrt(rno)*prt**0.666667_r8
          terms = ykf*sqrt(rno)*prs**0.666667_r8
          cht = utau/(tpr*log(zdz0)/kappa+termt)
          chs = utau/(tpr*log(zdz0)/kappa+terms)
          IF (ai(i,j) .le. min_a(ng)) THEN
            s0mk(i,j) = t(i,j,N(ng),nnew,isalt)
            t0mk(i,j) = t(i,j,N(ng),nnew,itemp)
          ELSE
            s0mk(i,j) = t(i,j,N(ng),nnew,isalt) -                       &
     &                         stflx(i,j,isalt)/chs
            s0mk(i,j) = max(s0mk(i,j),0._r8)
            s0mk(i,j) = min(s0mk(i,j),40._r8)
            t0mk(i,j) = t_freeze(s0mk(i,j),0.0_r8)
          END IF
! convert units to W/m^2
          wfr(i,j) = (wfr(i,j)                                          &
     &        -cht*(t(i,j,N(ng),nnew,itemp)-t0mk(i,j)))*rhocpr
        END DO
      END DO
#  else
      DO j=Jstr,Jend
        DO i=Istr,Iend
          heat = 0.0_r8
          thickness = 0.0_r8
          DO k=N(ng),1,-1
            dz = Hz(i,j,k)
            tfr = t_freeze(t(i,j,k,nnew,isalt),0.0_r8)
            IF ((thickness + dz) < bl_thick) THEN
              thickness = thickness + dz
              heat = heat + (t(i,j,k,nnew,itemp)-tfr)*dz
            ELSE
              fract = (bl_thick - thickness)/dz
!              thickness = thickness + dz*fract
              heat = heat + (t(i,j,k,nnew,itemp)-tfr)*dz*fract
            END IF
            IF (z_w(i,j,N(ng))-z_w(i,j,k-1) > bl_thick) EXIT
          END DO
! convert units to W/m^2
          wfr(i,j) = (wfr(i,j) - heat)*rhocpr/dt(ng)
        END DO
      END DO
#  endif
# endif

      CALL bc_r2d_tile (ng, tile,                                       &
     &                          LBi, UBi, LBj, UBj,                     &
     &                          wfr)
# ifdef DISTRIBUTE
      buffer(1) = exit_flag
      op_handle(1) = 'MAX'
      CALL mp_reduce (ng, iNLM, 1, buffer, op_handle)
      exit_flag = int(buffer(1))
      CALL mp_exchange2d (ng, tile, iNLM, 1,                            &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    NghostPoints, EWperiodic(ng), NSperiodic(ng), &
     &                    wfr)
# endif
      RETURN
      END SUBROUTINE ice_frazil_tile

#endif
      END MODULE ice_frazil_mod

