freezeup_khione.f

Go to the documentation of this file.

!                     **********************
                      MODULE freezeup_khione
!                     **********************
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    Module containing all subroutines to deal with the physics
!+        of freeze-up and frazil ice growth, at a node level.
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
      PRIVATE
      PUBLIC :: thermal_growth,melting_point,seeding,
     &          secondary_nucleation, flocculation_breakup,
     &          nusselt,turbulent_parameters,buoyancy_velocity,
     &          erosion_deposition,clogged_on_bar,precipitation
!
!=======================================================================
!
      CONTAINS
!
!=======================================================================
!
!                   *************************
                    SUBROUTINE thermal_growth
!                   *************************
!
     &  ( i,tn,tfrz,dt,srcgm_exp,srcgm_imp,sum_srcgm,sum_frzl,
     &    eps,alpha,nut,constss,limflag)
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Source term for frazil cristals thermal growth/decay
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| I          |-->| CELL INDEX
!| TN         |-->| TRACERS
!| TFRZ       |-->| FREEZING TEMPERATURE
!| DT         |-->| TIME STEP
!| SRCGM_EXP  |<->| EXPLICIT SOURCE OF FRAZIL
!| SRCGM_IMP  |<->| IMPLICIT SOURCE OF FRAZIL
!| SUM_SRCGM  |<->| SUM OF THERMAL GROWTH SOURCES (FOR TEMPERATURE SRC)
!| SUM_FRZL   |<->| SUM OF FRAZIL (TOTAL FRAIL CONC. FOR TEMPERATURE SRC)
!| EPS        |-->| TURBULENT DISSIPATION RATE
!| ALPHA      |-->| TURBULENT INTENSITY
!| NUT        |-->| TURBULENT VISCOSITY
!| CONSTSS    |-->| 1.D0/(RO0*CP_EAU)
!| LIMFLAG    |<->| LIMITOR FLAG
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY : nc_fra,ind_fra,ind_t,
     &  rho_ice,lh_ice,tc_wt,rk_frzl,ek_frzl,vk_frzl,
     &  nusst,nussi,nuss,itgm,minnk,iseed
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(IN)       :: TN
      INTEGER, INTENT(IN)              :: I
      LOGICAL, INTENT(INOUT)           :: LIMFLAG
      DOUBLE PRECISION, INTENT(IN)     :: DT,TFRZ,CONSTSS
      DOUBLE PRECISION, INTENT(IN)     :: EPS,ALPHA,NUT
      DOUBLE PRECISION, INTENT(INOUT)  :: SUM_SRCGM,SUM_FRZL
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCGM_EXP(nc_fra)
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCGM_IMP(nc_fra)
      DOUBLE PRECISION, PARAMETER      :: EPS0 = 1.d-8
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                          :: K,J
      DOUBLE PRECISION                 :: TW,QK
      DOUBLE PRECISION                 :: CF(nc_fra)
      DOUBLE PRECISION                 :: SRCGM(nc_fra)
      DOUBLE PRECISION                 :: RM1,RP1,FACT
      DOUBLE PRECISION                 :: MAXDT0,SRGMSTB
!
!-----------------------------------------------------------------------
!
!     TEMPERATURE AND FRAZIL CONC.
      tw = tn%ADR(ind_t)%P%R(i)
      DO k=1,nc_fra
        j = ind_fra+k-1
        cf(k) = tn%ADR(j)%P%R(i)
      ENDDO
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     COMPUTING THERMAL GROWTH/DECAY
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! ~~> LOOP ON ALL CLASS OF FRAZIL CRISTALS
      DO k=1,nc_fra
        j = ind_fra+k-1
!
! ~~>   COMPUTE NUSSLET NUMBER OF FRAZIL CLASS K
        IF (nussi.EQ.1) THEN
          nusst(k) = nuss
        ELSE IF (nussi.EQ.2) THEN
          nusst(k) = nusselt(rk_frzl(k), eps, alpha, nut)
        ENDIF
!
! ~~>   GROWTH
        IF(tw.LT.tfrz)THEN
          qk = nusst(k)*tc_wt/rk_frzl(k)
          srcgm(k) = 2.d0*qk/(rk_frzl(k)*rho_ice*lh_ice)
! ~~>   MELTING
        ELSE IF(tw.GT.tfrz)THEN
          qk = nusst(k)*tc_wt/rk_frzl(k)
          fact = (1.d0/ek_frzl(k)) + (1.d0/rk_frzl(k))
          srcgm(k) = 2.d0*qk*fact/(rho_ice*lh_ice)
          IF ((iseed.EQ.1).AND.
     &        (tn%ADR(j)%P%R(i).LE.minnk*vk_frzl(k)/nc_fra)) THEN
            srcgm(k) = 0.d0
          ENDIF
! ~~>   EQUILIBIRUM
        ELSE
          srcgm(k) = 0.d0
        ENDIF
!
        IF(itgm.EQ.1) THEN
          srcgm(k) = srcgm(k)*(tfrz-tw)*cf(k)
        ENDIF
!
      ENDDO
!
      srgmstb = srcgm(1)
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     COMPUTING MASS EXCHANGES BETWEEN CLASSES
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      limflag = .false.
      sum_srcgm = 0.d0
      sum_frzl = 0.d0
!
! ~~> LOOP ON ALL CLASSES OF FRAZIL CRISTALS
      DO k=1,nc_fra
        j = ind_fra+k-1
!
!       ~~~~~~~~~~~~~~~~~~~~~~~
!       EXPLICIT IMPLEMENTATION
!       ~~~~~~~~~~~~~~~~~~~~~~~
        IF(itgm.EQ.1) THEN
!
! ~~>     COMPUTING SOURCE TERM OF FRAZIL CONC. (MONO-CLASS)
          IF(nc_fra.EQ.1) THEN
            srcgm_exp(k) = max(-tn%ADR(j)%P%R(i)/dt, srcgm(k))
!
! ~~>     COMPUTING SOURCE TERM OF FRAZIL CONC. OF CLASS K
          ELSE
            rm1 = (rk_frzl(max(k-1, 1))/rk_frzl(k))**3
            rp1 = (rk_frzl(min(k+1,nc_fra))/rk_frzl(k))**3
!
!           GROWTH
            IF(tw.LT.tfrz)THEN
              IF(k.EQ.1) THEN
                srcgm_exp(k) =-srcgm(k)/(rp1-1.d0)
              ELSEIF(k.EQ.nc_fra) THEN
                srcgm_exp(k) = srcgm(k-1)/(1.d0-rm1)
              ELSE
                srcgm_exp(k) = srcgm(k-1)/(1.d0-rm1)-srcgm(k)/(rp1-1.d0)
              ENDIF
!
!           MELTING
            ELSEIF(tw.GT.tfrz)THEN
              IF(k.EQ.1) THEN
                srcgm_exp(k) = srcgm(k)-srcgm(k+1)/(rp1-1.d0)
              ELSEIF(k.EQ.nc_fra) THEN
                srcgm_exp(k) = srcgm(k)/(1.d0-rm1)
              ELSE
                srcgm_exp(k) = srcgm(k)/(1.d0-rm1)-srcgm(k+1)/(rp1-1.d0)
              ENDIF
!
!           EQUILIBRIUM
            ELSE
              srcgm_exp(k) = 0.d0
            ENDIF
!
!           LIMITOR
            IF (srcgm_exp(k).LE.-tn%ADR(j)%P%R(i)/dt)THEN
              limflag = .true.
              srcgm_exp(k) = -tn%ADR(j)%P%R(i)/dt
            ENDIF
          ENDIF
!
!         RATE OF TOTAL VOLUME FRACTION CHANGE
          sum_srcgm = sum_srcgm + srcgm_exp(k)
!
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       SEMI-IMPLICIT IMPLEMENTATION
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ELSEIF(itgm.EQ.2) THEN
!
! ~~>     COMPUTING SOURCE TERM OF FRAZIL CONC. (MONO-CLASS)
          IF(nc_fra.EQ.1) THEN
            srcgm_exp(k) = 0.d0
            srcgm_imp(k) = srcgm(k)*(tfrz-tw)
            srcgm(k) = srcgm(k)*cf(k)
!
! ~~>     COMPUTING SOURCE TERM OF FRAZIL CONC. OF CLASS K
          ELSE
            rm1 = (rk_frzl(max(k-1, 1))/rk_frzl(k))**3
            rp1 = (rk_frzl(min(k+1,nc_fra))/rk_frzl(k))**3
!
!           GROWTH
            IF(tw.LT.tfrz)THEN
              IF(k.EQ.1) THEN
                srcgm_exp(k) = 0.d0
                srcgm_imp(k) =-srcgm(k)*(tfrz-tw)/(rp1-1.d0)
              ELSEIF(k.EQ.nc_fra) THEN
                srcgm_exp(k) = srcgm(k-1)*(tfrz-tw)*cf(k-1)/(1.d0-rm1)
                srcgm_imp(k) = 0.d0
              ELSE
                srcgm_exp(k) = srcgm(k-1)*(tfrz-tw)*cf(k-1)/(1.d0-rm1)
                srcgm_imp(k) =-srcgm(k)*(tfrz-tw)/(rp1-1.d0)
              ENDIF
!
!           MELTING
            ELSEIF(tw.GT.tfrz)THEN
              IF(k.EQ.1) THEN
                srcgm_exp(k) =-srcgm(k+1)*(tfrz-tw)*cf(k+1)/(rp1-1.d0)
                srcgm_imp(k) = srcgm(k)*(tfrz-tw)
              ELSEIF(k.EQ.nc_fra) THEN
                srcgm_exp(k) = 0.d0
                srcgm_imp(k) = srcgm(k)*(tfrz-tw)/(1.d0-rm1)
              ELSE
                srcgm_exp(k) =-srcgm(k+1)*(tfrz-tw)*cf(k+1)/(rp1-1.d0)
                srcgm_imp(k) = srcgm(k)*(tfrz-tw)/(1.d0-rm1)
              ENDIF
!
!           EQUILIBRIUM
            ELSE
              srcgm_exp(k) = 0.d0
              srcgm_imp(k) = 0.d0
            ENDIF
!
!           LIMITOR
            IF (srcgm_exp(k).LE.-tn%ADR(j)%P%R(i)/dt)THEN
              limflag = .true.
              srcgm_exp(k) = -tn%ADR(j)%P%R(i)/dt
            ENDIF
          ENDIF
!
!         RATE OF TOTAL VOLUME FRACTION CHANGE
          IF (tfrz.NE.tw) THEN
            sum_srcgm = sum_srcgm + srcgm_exp(k)/(tfrz-tw)
     &                            + srcgm_imp(k)*cf(k)/(tfrz-tw)
          ENDIF
!
        ENDIF
!
!       TOTAL FRAZIL CONC.
        sum_frzl = sum_frzl + cf(k)
!
      ENDDO
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     STABILITY/POSITIVITY WARNING
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF(sum_srcgm.GT.eps0)THEN
        IF(itgm.EQ.1) THEN
          maxdt0 = min(2.d0/(rho_ice*lh_ice*constss*srgmstb
     &            *max(eps0, sum_frzl)),
     &             1.d0/(srgmstb*max(eps0, abs(tw))))
        ELSEIF(itgm.EQ.2) THEN
          maxdt0 = 1.d0/(srgmstb*max(eps0, abs(tw)))
        ENDIF
!
        IF( dt.GT.maxdt0) THEN
          WRITE(lu,*) ''
          WRITE(lu,*) 'STABILITY OF THERMAL GROWTH '
          WRITE(lu,*) 'OR POSITIVITY OF FRAZIL CONC. '
          WRITE(lu,*) 'MAY BE COMPROMISED '
          WRITE(lu,*) ''
          WRITE(lu,*) 'MAKE SURE TIME STEP '
          WRITE(lu,*) 'IS LOWER THAN ', maxdt0
          WRITE(lu,*) 'OR CHECK MODEL PARAMETERS '
          WRITE(lu,*) ''
        ENDIF
      ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE thermal_growth
!
!=======================================================================
!
!                     ******************
                      SUBROUTINE seeding
!                     ******************
!
     &  ( i, tn, tfrz, srcse, minflag )
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes frazil ice seeding
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| I          |-->| CELL INDEX
!| TN         |<->| TRACERS
!| TFRZ       |-->| FREEZING TEMPERATURE
!| SRCSE      |<->| EXPLICIT SOURCE OF FRAZIL
!| MINFLAG    |<->| LOGICAL FLAG TO DISABLE SRCSN WHEN MIN CONC. REACHED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY : nc_fra,ind_fra,ind_t,vk_frzl,
     &  iseed,seedr,minnk
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(INOUT)    :: TN
      LOGICAL, INTENT(INOUT)           :: MINFLAG
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCSE(nc_fra)
      DOUBLE PRECISION, INTENT(IN)     :: TFRZ
      INTEGER, INTENT(IN)              :: I
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                          :: K,J,FROZNB
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     ~~~~~~~~~~
!     NO SEEDING
!     ~~~~~~~~~~
      IF(iseed.EQ.0) THEN
        DO k=1,nc_fra
          srcse(k) = 0.d0
        ENDDO

!     ~~~~~~~~~~~~~~~~~~~~~~~
!     MINIMUM CONC. THRESHOLD
!     ~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF((iseed.EQ.1).OR.(iseed.EQ.3)) THEN
        froznb = 0
        DO k=1,nc_fra
          j = ind_fra+k-1
          srcse(k) = 0.d0
!
!         RMK: TRACER VALUES NEED TO BE MODIFIED DIRECTLY OTHERWISE
!         SEMI-IMPLICIT THERMAL GROWTH WOULD NOT WORK AS INTENDED
          tn%ADR(j)%P%R(i) = max(minnk*vk_frzl(k)/nc_fra,
     &                           tn%ADR(j)%P%R(i))
!
!         IN CASE OF MELTING SECONDARY NUCLEATION AND FLOCCULATION
!         ARE NOT COMPUTED IF ALL FRAZIL CLASS CONC. ARE LE MIN CONC.
!         (ONLY WHEN MIN CONC THRESHOLD IS SELECTED AS SEEDING METHOD)
          IF ((tn%ADR(ind_t)%P%R(i).GT.tfrz).AND.
     &        (tn%ADR(j)%P%R(i).LE.minnk*vk_frzl(k)/nc_fra)) THEN
            froznb = froznb + 1
          ENDIF
        ENDDO
        IF (froznb.EQ.nc_fra) THEN
          minflag = .true.
        ELSE
          minflag = .false.
        ENDIF
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     CONSTANT SEEDING RATE PER CLASS
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF((iseed.EQ.2).OR.(iseed.EQ.3)) THEN
        DO k=1,nc_fra
          srcse(k) = vk_frzl(k)*seedr
        ENDDO
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~
!     SEEDING RATE COMPUTED
!     FROM ATMOSPHERIC DRIVERS
!     ~~~~~~~~~~~~~~~~~~~~~~~~
!     TODO SRCSE(K) = SEED%R(I)
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE seeding
!
!=======================================================================
!
!             *******************************
              SUBROUTINE secondary_nucleation
!             *******************************
!
     &  ( i, tn, dt, srcsn, eps, limflag, minflag )
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes secondary nucleation
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| I          |-->| CELL INDEX
!| TN         |-->| TRACERS
!| DT         |-->| TIME STEP
!| SRCSN      |<->| EXPLICIT SOURCE OF FRAZIL
!| EPS        |-->| TURBULENT DISSIPATION RATE
!| LIMFLAG    |-->| LOGICAL FLAG TO DISABLE SRCSN WHEN MAX SRCGM REACHED
!| MINFLAG    |-->| LOGICAL FLAG TO DISABLE SRCSN WHEN MIN CONC. REACHED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY : nc_fra, ind_fra, isnuc, xnu,
     &  rk_frzl, vk_frzl, vbk, snnmax
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(IN)       :: TN
      INTEGER, INTENT(IN)              :: I
      LOGICAL, INTENT(IN)              :: LIMFLAG, MINFLAG
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCSN(nc_fra)
      DOUBLE PRECISION, INTENT(IN)     :: DT,EPS
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                          :: K,J
      DOUBLE PRECISION                 :: SUMSN,WTSQ,CF,WRK,NTILDE
      DOUBLE PRECISION, PARAMETER      :: PI = 4.d0*atan(1.d0)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     ~~~~~~~~~~~~~~~~~~~~~~~
!     NO SECONDARY NUCLEATION
!     ~~~~~~~~~~~~~~~~~~~~~~~
      IF((isnuc.EQ.0).OR.(nc_fra.EQ.1).OR.limflag.OR.minflag) THEN
        DO k=1,nc_fra
          srcsn(k) = 0.d0
        ENDDO
!
!     ~~~~~~~~~~~~~~~~~~~~~~~
!     SVENSSON & OMSTEDT 1994
!     ~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF(isnuc.EQ.1) THEN
!
!       COMPUTATION OF AVERAGE NUMBER OF ICE CRYSTALS
        ntilde = 0.d0
        DO k=1,nc_fra
          j = ind_fra+k-1
          cf = tn%ADR(j)%P%R(i)
          ntilde = ntilde + cf/vk_frzl(k)
        ENDDO
        ntilde = min(snnmax, ntilde)
!
        IF(ntilde.GT.1.d0) THEN
!         LOOP ON ALL CLASSES EXCEPT 1
          sumsn = 0.d0
          DO k=2,nc_fra
            j = ind_fra+k-1
            cf = tn%ADR(j)%P%R(i)
!           TURBULENT COMPONENT OF RELATIVE VELOCITY
            wtsq = (4.d0/15.d0)*(eps/xnu)*rk_frzl(k)**2
!           RELATIVE VELOCITY
            wrk = sqrt(wtsq + vbk(k)**2)
!           NUCLEATION SINK TERM
            srcsn(k) = -pi*ntilde*wrk*cf*rk_frzl(k)**2
!           LIMITOR
            srcsn(k) = max(-cf/dt, srcsn(k))
!           SUM FOR FIRST CLASS SOURCE
            sumsn = sumsn + srcsn(k)
          ENDDO
!         NUCLEATION SOURCE TERM FOR FIRST CLASS
          srcsn(1) = -sumsn
        ELSE
          DO k=1,nc_fra
            srcsn(k) = 0.d0
          ENDDO
        ENDIF
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~
!     WANG & DOERING 2005
!     WARNING: NOT CONSERVATIVE
!     ~~~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF(isnuc.EQ.2) THEN
!
!       COMPUTATION OF AVERAGE NUMBER OF ICE CRYSTALS
        ntilde = 0.d0
        DO k=1,nc_fra
          j = ind_fra+k-1
          cf = tn%ADR(j)%P%R(i)
          ntilde = ntilde + cf/vk_frzl(k)
        ENDDO
        ntilde = min(snnmax, ntilde)
!
        IF(ntilde.GT.1.d0) THEN
!         LOOP ON ALL CLASSES EXCEPT 1
          sumsn = 0.d0
          DO k=2,nc_fra
            j = ind_fra+k-1
            cf = tn%ADR(j)%P%R(i)
!           TURBULENT COMPONENT OF RELATIVE VELOCITY
            wtsq = (4.d0/15.d0)*(eps/xnu)*rk_frzl(k)**2
!           RELATIVE VELOCITY
            wrk = sqrt(wtsq + vbk(k)**2)
!           NUCLEATION SINK TERM
            srcsn(k) = -pi*ntilde*wrk*cf*rk_frzl(k)**2
!           LIMITOR
            srcsn(k) = max(-cf/dt, srcsn(k))
!           SUM FOR FIRST CLASS SOURCE
            sumsn = sumsn + srcsn(k)
            srcsn(k) = srcsn(k)*(rk_frzl(1)/rk_frzl(k))**3
          ENDDO
!         NUCLEATION SOURCE TERM FOR FIRST CLASS
          srcsn(1) = -sumsn
        ELSE
          DO k=1,nc_fra
            srcsn(k) = 0.d0
          ENDDO
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE secondary_nucleation
!
!=======================================================================
!
!             *******************************
              SUBROUTINE flocculation_breakup
!             *******************************
!
     &  ( i, tn, dt, srcfb, limflag, minflag )
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes floculation & breaking
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| I          |-->| CELL INDEX
!| TN         |-->| TRACERS
!| DT         |-->| TIME STEP
!| SRCFB      |<->| EXPLICIT SOURCE OF FRAZIL
!| LIMFLAG    |-->| LOGICAL FLAG TO DISABLE SRCFB WHEN MAX SRCGM REACHED
!| MINFLAG    |-->| LOGICAL FLAG TO DISABLE SRCFB WHEN MIN CONC. REACHED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY : nc_fra, ind_fra, rk_frzl,
     &  ifloc, afloc
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(IN)       :: TN
      LOGICAL, INTENT(IN)              :: LIMFLAG, MINFLAG
      INTEGER, INTENT(IN)              :: I
      DOUBLE PRECISION, INTENT(IN)     :: DT
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCFB(nc_fra)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                          :: K,J
      DOUBLE PRECISION                 :: CF, CF0
      DOUBLE PRECISION                 :: BETA, BETA0
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     ~~~~~~~~~~~~~~~~~~~~~~~
!     NO FLOCULATION & BREAUP
!     ~~~~~~~~~~~~~~~~~~~~~~~
      IF((ifloc.EQ.0).OR.(nc_fra.EQ.1).OR.limflag.OR.minflag) THEN
        DO k=1,nc_fra
          srcfb(k) = 0.d0
        ENDDO
!
!     ~~~~~~~~~~~~~~~~~~~~~~~
!     SVENSSON & OMSTEDT 1994
!     ~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF(ifloc.EQ.1) THEN
        DO k=1,nc_fra
          j = ind_fra+k-1
!
          IF(k.EQ.1) THEN
            cf = tn%ADR(j)%P%R(i)
            beta = afloc*rk_frzl(k)/rk_frzl(1)
            beta = min(1.d0/dt, beta) !LIMITOR
            srcfb(k) = -beta*cf
!
          ELSE IF(k.EQ.nc_fra) THEN
            beta0 = beta
            cf0 = tn%ADR(j-1)%P%R(i)
            srcfb(k) = beta0*cf0
!
          ELSE
            cf = tn%ADR(j)%P%R(i)
            cf0 = tn%ADR(j-1)%P%R(i)
            beta0 = beta
            beta = afloc*rk_frzl(k)/rk_frzl(1)
            beta = min(1.d0/dt, beta) !LIMITOR
            srcfb(k) = -beta*cf + beta0*cf0
          ENDIF
        ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE flocculation_breakup
!
!=======================================================================
!
!             *******************************
              SUBROUTINE turbulent_parameters
!             *******************************
!
     &  ( i, vmag, h, kt, eps, alpha, nut, cf, ak, ep, iturb_tel)
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes turbulent parameters for frazil ice model
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AK       |-->| TURBULENT KINETIC ENERGY K AT TIME T(N+1)
!| EP       |-->| TURBULENT ENERGY DISSIPATION AT TIME T(N+1)
!| ITURB_TEL|-->| TURBULENT MODEL IN T2D
!| VMAG     |-->| VELOCITY MAGNITUDE
!| H        |-->| WATER DEPTH
!| KT       |<->| MEAN TURBULENT KINETIC ENERGY
!| EPS      |<->| MEAN TURBULENT DISSIPATION RATE
!| ALPHA    |<->| TURBULENT INTENSITY
!| NUT      |<->| TURBULENT VISCOSITY
!| CF       |<->| QUADRATIC FRICTION COEFFICIENT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY : xnu, iturb
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj), INTENT(IN)       :: AK,EP
      INTEGER, INTENT(IN)              :: I,ITURB_TEL
      DOUBLE PRECISION, INTENT(IN)     :: VMAG,H,CF
      DOUBLE PRECISION, INTENT(INOUT)  :: KT,EPS,ALPHA,NUT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, PARAMETER      :: KARMAN_CONST = 0.4d0
      DOUBLE PRECISION, PARAMETER      :: HEPS = 1.d-3
      DOUBLE PRECISION, PARAMETER      :: UEPS = 1.d-6
      DOUBLE PRECISION                 :: USTAR,FACT0,FACT1,HTEMP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     CONSTANTS (FOR TEST PURPOSES)
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF(iturb.EQ.0) THEN
        kt = 9.6d-4
        eps = 12.d-4
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     DEPTH AVERAGED K, EPS
!     FROM MIXING LENGTH VERTICAL PROFILE
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF(iturb.EQ.1) THEN
        ustar = 0.5d0 * cf * vmag * vmag
        ustar = sqrt(ustar)
        htemp = max(h, heps)
        fact0 = (htemp/2.d0 - xnu/ustar + xnu**2/(2.d0*htemp*ustar**2))
        fact1 = (log(htemp*ustar/xnu) - 1.d0 + xnu/(htemp*ustar))
        kt = (ustar**2)*fact0/0.3d0
        eps = (ustar**3)*fact1/karman_const
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~
!     DEPTH AVERAGED K, EPS
!     FROM TELEMAC2D K-EPS MODEL
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~
      ELSE IF(iturb.EQ.2) THEN
        IF(iturb_tel.NE.3) THEN
          WRITE(lu,*) ''
          WRITE(lu,*) 'TURBULENCE ESTIMATION MODEL', iturb
          WRITE(lu,*) 'MUST BE USED WITH K-EPS MODEL'
          WRITE(lu,*) ''
          WRITE(lu,*) 'PLEASE SET TURBULENCE MODEL = 3'
          WRITE(lu,*) 'IN T2D STEERING FILE'
          WRITE(lu,*) ''
          CALL plante(1)
          stop
        ENDIF
        kt = ak%R(i)
        eps = ep%R(i)
!
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~
      ELSE
        WRITE(lu,*) ''
        WRITE(lu,*) 'TURBULENCE ESTIMATION MODEL', iturb
        WRITE(lu,*) 'NOT IMPLEMENTED'
        WRITE(lu,*) ''
        CALL plante(1)
        stop
      ENDIF
!
      nut = 0.09d0*kt**2/eps
      alpha = sqrt(2.d0*kt)/max(ueps, vmag)
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE turbulent_parameters
!
!=======================================================================
!
!             *********************************
              DOUBLE PRECISION FUNCTION nusselt
!             *********************************
!
     &  ( radius, eps, alpha, nut )
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes the Nusselt number for thermal growth
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| RADIUS   |-->| RADIUS OF FRAZIL PARTICLE
!| EPS      |-->| TURBULENT DISSIPATION RATE
!| ALPHA    |-->| TURBULENT INTENSITY
!| NUT      |-->| TURBULENT VISCOSITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY : tc_wt, cp_eau, ro0, xnu
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)     :: RADIUS, ALPHA, EPS, NUT
      DOUBLE PRECISION                 :: MSTAR, KOLML, PR, PR12, PR13
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      kolml = (xnu**3/eps)**(0.25d0)
      mstar = radius/kolml
      pr = ro0*cp_eau*nut/tc_wt
      pr12 = sqrt(pr)
      pr13 = pr**(1.d0/3.d0)
!
      IF (mstar .LE. 1.d0/pr12) THEN
        nusselt = 1.d0 + 0.17d0*mstar*pr12
      ELSE IF ((1.d0/pr12 .LT. mstar) .AND. (mstar .LE. 1.d0)) THEN
        nusselt = 1.d0 + 0.55d0*(mstar**(2.d0/3.d0))*pr13
      ELSE IF (mstar .GT. 1.d0) THEN
        IF (alpha*mstar**(4.d0/3.d0) .LE. 1.d3) THEN
          nusselt = 1.1d0 + 0.77d0*(alpha**0.035d0)
     &              *(mstar**(2.d0/3.d0))*pr13
        ELSE
          nusselt = 1.1d0 + 0.77d0*(alpha**0.25d0)*mstar*pr13
        ENDIF
      ELSE
        nusselt = 1.d0
      ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END FUNCTION nusselt
!
!=======================================================================
!
!           *******************************************
            DOUBLE PRECISION FUNCTION buoyancy_velocity
!           *******************************************
!
     &  ( radius, thickness )
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes the buoyancy velocity of frazil particles
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| RADIUS    |-->| RADIUS OF FRAZIL PARTICLE
!| THICKNESS |-->| THICKNESS OF FRAZIL PARTICLE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY : ibuoy, xnu, de,
     &                                ro0, rho_ice
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)     :: RADIUS, THICKNESS
      DOUBLE PRECISION, PARAMETER      :: CKD = 2.d0
      DOUBLE PRECISION, PARAMETER      :: GRAV = 9.81d0
      DOUBLE PRECISION                 :: KV,GP,PI
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      IF(ibuoy.EQ.1) THEN
!       DALY (1984)
        pi = 4.d0*atan(1.d0)
        kv = 2.d0*pi/de
        gp = 2.d0*kv*grav*(ro0-rho_ice)/(pi*ro0)
        IF (radius .LE. 3.e-4) THEN
            buoyancy_velocity = 0.08d0*gp*(radius**2)/xnu
        ELSEIF ((3.e-4 .LT. radius).AND.(radius .LE. 1.4e-3)) THEN
            buoyancy_velocity = 0.16d0*(gp**0.715d0)
     &                        *(xnu**(-0.428d0))*(radius**1.14d0)
        ELSEIF (1.4e-3 .LT. radius) THEN
            buoyancy_velocity = sqrt(gp*radius/2.d0)
        ENDIF
!
      ELSEIF(ibuoy.EQ.2) THEN
!       DALY (1984) INTERMEDIATE RANGE
!       RMK: VALID FOR 0.03 < R < 0.14 CM
        buoyancy_velocity = 1.d-2*30.d0*(radius*1.d2)**1.2d0

      ELSEIF(ibuoy.EQ.3) THEN
!       MATOUSEK (1992)
        buoyancy_velocity = (1.31d-5)*((2.d0*radius)**0.29d0)
     &                               *(thickness**0.61d0)/xnu

      ELSEIF(ibuoy.EQ.4) THEN
!       GOSIK & OSTERKAMP (1983)
!       RMK1: STOKES REGIME (ONLY VALID FOR R < 0.03 CM)
!       RMK2: DRAG COEFFICIENT CKD NEEDS TO BE ADJUSTED
        buoyancy_velocity = sqrt(4.d0*(ro0-rho_ice)
     &                      *grav*radius/(de*ro0*ckd))
!
      ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END FUNCTION buoyancy_velocity
!
!=======================================================================
!
!             ***************************************
              DOUBLE PRECISION FUNCTION melting_point
!             ***************************************
!
     &  ( sal )
!
!***********************************************************************
! KHIONE   V7P2
!***********************************************************************
!
!brief    Computes the melting point as a function of salinity
!
!reference:
!+   Millero, F.J. “Freezing point of seawater”. In “Eighth Report of
!+     the Joint Panel on Oceanographic Tables and Standards”,
!+     UNESCO Tech. Pap. Mar. Sci. No.28, Annex 6. UNESCO, Paris, 1978
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| SAL   |-->| SALINITY CONCENTRATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)     :: SAL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      melting_point =
     &  -0.0575d0*sal + 0.001710523d0*( sal**1.5d0 )
     &  -0.0002154996d0*( sal**2 )
!    &  -0.00753 * PRESSURE BELOW SURFACE
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END FUNCTION melting_point
!
!=======================================================================
!
!                   *****************************
                    SUBROUTINE erosion_deposition
!                   *****************************
!
     &  ( frzl,srcgm,theta0,theta1,beta1,vbb,thifemf,hun,
     &    anfem,vmag,depth,isbar )
!
!***********************************************************************
! KHIONE   V7P2
!***********************************************************************
!
!brief    Computes frazil ice sink/source term due to deposition/erosion
!         under the ice cover
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ANFEM   |-->| CONCENTRATION OF SURFACE ICE PARTICLES
!| BETA1   |-->| RATE OF REENTRAINMENT OF SURFACE ICE PER UNIT AREA
!| DEPTH   |-->| WATER DEPTH
!| FRZL    |-->| FRAZIL CONCENTRATION
!| HUN     |-->| UNCOVERED ICE THICKNESS
!| ISBAR   |-->| WHETHER OR NOT ICE SETTLES AT THE SURFACE OR ON A BAR
!| SRCGM    |<->| EXPLICIT SOURCE OF FRAZIL
!| THETA0  |-->| PROBABILITY OF DEPOSITION OF FRAZIL - OPEN WATER
!| THETA1  |-->| PROBABILITY OF DEPOSITION OF FRAZIL - ICE
!| THIFEMF |-->| FRAZIL ICE THICKNESS
!| VBB     |-->| BUOYANCY VELOCITY OF FRAZIL GRANULES
!| VMAG    |-->| CURRENT VELOCITY MAGNITUDE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY : af, surf_ef
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)     :: FRZL,ANFEM
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCGM
      DOUBLE PRECISION, INTENT(IN)     :: THETA0,THETA1,BETA1
      DOUBLE PRECISION, INTENT(IN)     :: VMAG,DEPTH,VBB,THIFEMF,HUN
      LOGICAL, INTENT(IN)              :: ISBAR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION HE,PV
!
!-----------------------------------------------------------------------
!
!     NODES ON THE BAR: NO DEPOSITION ON THE WATER SURFACE
      IF( isbar ) THEN
        pv = ( af*vmag )/depth * frzl
!     NO BAR: DEPOSITION ON WATER SURFACE
      ELSE
        pv = ( (theta0*vbb*(1.d0-anfem) ) +
     &    ( theta1*vbb*anfem) )/depth * frzl
      ENDIF
!
!     EROSION / DEPOSITION SOURCE TERM
      he = ( 1.d0-surf_ef ) * ( thifemf + hun )
      srcgm = (beta1*he*anfem)/depth - pv
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE erosion_deposition
!
!=======================================================================
!
!                   ************************
                    SUBROUTINE precipitation
!                   ************************
!
     &  (vbb, frzl, srcp, vmag, g, rk, it)
!
!***********************************************************************
! KHIONE   V8P2
!***********************************************************************
!
!brief    Computes frazil precipitation on the surface
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FRZL    |-->| FRAZIL CONCENTRATION
!| G       |-->| GRAVITY CONSTANT
!| IT      |-->| LOGICAL FOR ICE COVER
!| RK      |-->| FRAZIL RADIUS
!| SRCP    |<->| EXPLICIT SOURCE OF FRAZIL
!| VBB     |-->| BUOYANCY VELOCITY OF FRAZIL GRANULES
!| VMAG    |-->| CURRENT VELOCITY MAGNITUDE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY : ro0, rho_ice
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL, INTENT(IN)              :: IT
      DOUBLE PRECISION, INTENT(IN)     :: FRZL,G,RK
      DOUBLE PRECISION, INTENT(INOUT)  :: SRCP
      DOUBLE PRECISION, INTENT(IN)     :: VMAG,VBB
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION UC2,TC,M
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     COMPUTATION OF THE CRITICAL VELOCITY (USEFUL IF ICE COVER IS
!     PRESENT)
      IF(it) THEN
        tc = 2.d2/3.d0
        uc2 = tc*(ro0-rho_ice)*g*rk/ro0
        m = max(0.d0,1.d0-(vmag*vmag/uc2))
        srcp = -vbb*frzl*m
      ELSE
        srcp = -vbb*frzl
      ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE precipitation
!
!=======================================================================
!
!                   *************************
                    SUBROUTINE clogged_on_bar
!                   *************************
!
     &  ( rfr0,rfr1,db,bar,nbar,ang1,fm1,fmt )
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    Computes the accumulated mass on vertical and/or
!         transverse bars.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ANG1   |-->| ACCUMULATION ANGLE
!| BAR    |-->| PHYSICAL SIZE OF LENGTHS OF THE BAR
!| DB     |-->| BAR DIAMETER
!| FM1    |<->| MASS OF FRAZIL ICE ACCULATED ON ONE BAR
!| FMT    |<->| TOTAL MASS ACCUMULATED
!| NBAR   |-->| REPRESENTATIVE NUMBER OF BAR
!| RFR0   |-->| INITIAL RADIUS BEFORE THE FRAZIL ICE ACCUMULATION
!| RFR1   |-->| RADIUS OF THE FRAZIL ICE ACCUMULATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY : rho_ice,clog_ef
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)    :: NBAR
      DOUBLE PRECISION, INTENT(IN)    :: RFR0,RFR1,BAR,ANG1,DB
      DOUBLE PRECISION, INTENT(INOUT) :: FM1,FMT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION :: RB,DR,AMNU,ASH,RATIO,ADONT,AICE
!
!-----------------------------------------------------------------------
!
!      RB = RFR0/2.0/COS(ANG1)
      rb = db/2.d0
      dr = 2.d0*rb - rfr0
!
      amnu = rfr0**2*cos(ang1)*sin(ang1) + 2.d0*ang1*rb**2 -
     &       rb**2*cos(2.d0*ang1)*sin(2.d0*ang1)
      ash = 2.d0*ang1*rb**2 - (ang1*rfr0**2 - rb*rfr0*sin(ang1)) ! SHADE AREA
!
      IF (rfr1.GT.2.d0*rb) THEN
        aice = ang1*rfr1**2 - amnu
      ELSE
        adont = ( rfr1**2 - rfr0**2 )*ang1
        ratio = ( rfr1 - rfr0 ) / dr
        aice = adont - ratio*ash
      ENDIF
!
      fm1 = aice*bar*rho_ice*( 1.d0-clog_ef )
      fmt = fm1*nbar
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE clogged_on_bar
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      END MODULE freezeup_khione