border_icover.f

Go to the documentation of this file.

!                     ************************
                      SUBROUTINE border_icover
!                     ************************
!
     &( u,v, mesh )
!
!***********************************************************************
! KHIONE   V8P0
!***********************************************************************
!
!brief    COMPUTES PRESENCE OF STATIC BORDER ICE
!
!+  Use of ICETYPE:
!+    1 - open water
!+    2 - static border ice (within or at the edge of the cover)
!+    3 - border ice (both static and dynamic, and including the edge)
!+  Use of ICELOC:
!+    1 - open water
!+    2 - domain boundary node
!+    3 - edge of the border ice where dynamic border ice accumulates
!+  Note that dynamic border ice (where ice is being accumulated),
!+  TODO: This method still need to be updated to account for the
!+    melting down of border ice (in patches if possible)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| U              |-->| X-COMPONENT OF THE VELOCITY
!| V              |-->| Y-COMPONENT OF THE VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY:bd_ice,anfem,thifemf,thifems,
     &                              anfem0,it1,it2,t1,
     &                              icetype,vcrbor,vcrbom,tcr,vz,
     &                              tc,tmelt,iceloc,icetypep,t2
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(IN)    :: U,V
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          I,J,K,N1,N2,IT,NPOIN,NELEM,NELMAX,IL,ITP
      DOUBLE PRECISION THICK, USTAR,VMAG,VB
!
!-----------------------------------------------------------------------
!
      npoin = mesh%NPOIN
      nelem = mesh%NELEM
      nelmax = mesh%NELMAX
!
!
!=======================================================================
!
!     7 - STATIC BORDER ICE GROWTH
!
!-----------------------------------------------------------------------
!
      IF( bd_ice ) THEN
!
!-----------------------------------------------------------------------
!
!     MAYBE THIS SHOULD BE DONE ONCE /!\
!     TODO: T1 CAN BE REPLACED BY VOLUPAR (ALREADY COMPUTED)
!
        CALL vector(t1,'=','MASBAS          ',u%ELM,1.d0,
     &              t2,t2,t2,t2,t2,t2,mesh,.false.,t2)
!       ASSEMBLING THE SUM OF ALL VALUES
        IF( ncsize.GT.1 ) CALL parcom( t1,2,mesh )
!
!-----------------------------------------------------------------------
!
!     LOOKING FOR NEW BORDER ICE NODES (TO BE SWITCHED))
!
        DO i = 1,npoin
          it = icetype%I(i)
          il = iceloc%I(i)
          itp = icetypep%I(i)
!
!         CURRENT VELOCITY
          vmag = sqrt( u%R(i)**2 + v%R(i)**2 )
!         BUOYANT VELOCITY
          vb = max( 0.d0, -0.025d0*tcr%R(i) + 0.005d0 )
!         BORDER ICE THICKNESS /!\ TODO: DOUBLE CHECK VALIDITY
          thick = thifems%R(i)+thifemf%R(i)
!         VELOCITY CRITERIA FOR DYNAMIC BORDER ICE GROWTH
          ustar = max( 0.175d0, vmag / max(vcrbom,1.d-12) )
!
!         IT1 IS USED AS A STATIC ICE MASK DEFINING NODES FOR WHICH
!         THE THRESHOLDS ARE TRUE (I.E. SUBJECT TO BE BODER ICE NODE).
!         - IT1 = 0 IF THRESHOLDS ARE --NOT-- TRUE
!         - IT1 = 1 IF THRESHOLDS ARE TRUE
!         - IT1 = 2 IF THRESHOLD ARE TRUE AND IS ALSO ON THE DOMAIN
!         BORDER IT IS NOTED THAT IF IT1 = 2, THEN THE NODE IS
!         --NOT-- BORDER ICE ALREADY, BUT ABOUT TO BE SWITCHED.
          it1%I(i) = 0
!
!     > THRESHOLDS FOR STATIC BORDER ICE FORMATION
          IF( tcr%R(i).LE.(tc-tmelt%R(i)) .AND. ! THERMAL PROPERTY
     &        vb.GT.1.1d0*vz%R(i) .AND.         ! BUOYANT VS. TURBULENCE
     &        vmag.LT.vcrbor ) THEN             ! CRITICAL VELOCITY
            IF( it .NE. 2 ) THEN
!           NODE OUTSIDE STATIC BORDER ICE COVER (NOT MULTIPLIER OF 3)
              it1%I(i) = 1
!             BORDER NODE ABOUT TO BE SWITCH TO STATIC BORDER ICE
              IF( il.EQ.2.OR.il.EQ.4 ) it1%I(i) = 2
            ENDIF
          ENDIF
!
!         IT2 IS USED AS A DYNAMIC ICE MASK DEFINING NODES FOR WHICH
!         THE THRESHOLDS ARE TRUE (I.E. SUBJECT TO BE BODER ICE NODE).
!         - IT2 = 0 IF THRESHOLDS ARE --NOT-- TRUE
!         - IT2 = 1 IF THRESHOLDS ARE TRUE
!         - IT2 = 2 IF THRESHOLD ARE TRUE AND IS ALSO ON THE DOMAIN
!         BORDER IT IS NOTED THAT IF IT2 = 2, THEN THE NODE IS
!         --NOT-- BORDER ICE ALREADY, BUT ABOUT TO BE SWITCHED.
          it2%I(i) = 0
!
!     > THRESHOLDS FOR DYNAMIC BORDER ICE FORMATION
          IF( !THICK.GE.0.0D0 .AND.              ! CRITICAL THICKNESS
     &        ustar.LE.1.d0 ) THEN              ! CRITICAL VELOCITY
            IF( il .LT. 3.AND.itp.NE.2  ) THEN
              it2%I(i) = 1
!             BORDER NODE ABOUT TO BE SWITCH TO STATIC BORDER ICE
              IF( il.EQ.2 ) it2%I(i) = 2
!
!            ELSE !IF( IL.GE.3 ) THEN
!     > FURTHER GROWING DYNAMIC ICE
!              R = 14.1 * USTAR ** (-0.93) * ANFEM%R(I) ** 1.08
!              DPHI = LIN_WATAIR*( TWAT%R(I) - TAIR%R(I) )
!              DW = DPHI * R / ( ROEAU * LH_ICE )
!              DWB%R(I) = MIN( 1.D0, DWB%R(I) + DW / SQRT( T1%R(I) ) )
!
            ENDIF
          ENDIF
!
        ENDDO
!
!-----------------------------------------------------------------------
!
!     CHECK BORDER NODES ABOUT TO BE SWITCHED INTO BORDER ICE NODES
!
        DO j = 1,nelem
!
          DO k = 1,3
! RA: why we compute again node indices ? They are stored in IKLE (NELEM,3)
!      1- introduce IKLE as an argument like IKLE(NELEM,3)
!      2- nodes of element IELEM are :
!              I1 = IKLE(IELEM,1)
!              I2 = IKLE(IELEM,2)
!              I3 = IKLE(IELEM,3)
!
            i = mesh%IKLE%I( j + (k-1)*nelmax )
!
            it = icetype%I(i)
            itp = icetypep%I(i)
            il = iceloc%I(i)
!           NODES ABOUT TO BE SWITCHED NEED TO BE BELOW THRESHOLDS AND
!           --NOT-- ALREADY BODER ICE NODES. THEY ALSO NEED TO BE CLOSE
!           TO TWO OTHER (SUFFICIENTLY THICK) BORDER ICE NODES.
            IF( it1%I(i).EQ.1 ) THEN
              n1 = icetype%I( mesh%IKLE%I( j + mod(k,3)*nelmax ) )
              n2 = icetype%I( mesh%IKLE%I( j + mod(k+1,3)*nelmax ) )
              IF( n1 .EQ. 2 .AND. n2 .EQ. 2 ) THEN
                it1%I(i) = 2
              ENDIF
            ENDIF
            IF( it2%I(i).EQ.1 ) THEN
              n1 = icetypep%I( mesh%IKLE%I( j + mod(k,3)*nelmax ) )
              n2 = icetypep%I( mesh%IKLE%I( j + mod(k+1,3)*nelmax ) )
              IF( n1 .EQ. 2 .AND. n2 .EQ. 2 ) THEN
                it2%I(i) = 2
              ENDIF
            ENDIF
          ENDDO
!
        ENDDO
!
!     ASSEMBLING THE MAX OF ALL VALUES - THE TWOS WILL BE SWITCHED
!     CONVERSION BETWEEN INTEGER AND DOUBLE IS NECESSARY FOR PARCOM
        IF(ncsize.GT.1)CALL parcom2i(it1%I,it2%I,it2%I,npoin,1,3,2,mesh)
!
!-----------------------------------------------------------------------
!
!     SWITCHES BORDER NODES INTO STATIC BORDER ICE NODES
!
        DO i = 1,npoin
          it = icetype%I(i)
          itp = icetypep%I(i)
          il = iceloc%I(i)
!
          IF( it1%I(i).EQ.2 ) THEN
!     > STATIC BORDER ICE
            icetype%I(i) = 2           !=> INSTANT SWITCH
!
!         SWITCH FROM SOLID DYNAMIC TO STATIC BORDER ICE
            IF( itp .EQ. 2 ) THEN
!
!         SWITCH FROM GROWING DYNAMIC TO STATIC BORDER ICE
            ELSEIF( il .GE. 3 ) THEN
!           USE DWB TO COMPUTE THIFEMS, THIFEMF AND ANFEM ?
              IF(il.EQ.3) iceloc%I(i) = 1         !=> INSTANT SWITCH
              IF(il.EQ.4) iceloc%I(i) = 2         !=> INSTANT SWITCH
              icetypep%I(i) = 2
!
!         NEW STATIC BORDER ICE
            ELSE
!           NOT PART OF THE DYNAMIC ICE COVER
              icetypep%I(i) = 2         !=> INSTANT SWITCH
              thifems%R(i) = 0.d0
              thifemf%R(i) = 0.d0
              anfem%R(i) = anfem0
!
            ENDIF
!
!     > DYNAMIC BORDER ICE
          ELSEIF( it2%I(i).EQ.2 ) THEN
            IF( il .LT. 3 ) THEN
              thifems%R(i) = 0.d0
              thifemf%R(i) = 0.d0
              anfem%R(i) = anfem0
              IF(il.EQ.1) iceloc%I(i) = 3         !=> INSTANT SWITCH
              IF(il.EQ.2) iceloc%I(i) = 4         !=> INSTANT SWITCH
            ENDIF
          ENDIF
!
        ENDDO
!
      ENDIF
!
!
!=======================================================================
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE border_icover