source_icover.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE source_icover
!                   ************************
     &(npoin,fu,fv,h,u,v,t1,t2,t3,grav,karman,chestr,dt,at)
!
!***********************************************************************
! KHIONE   V7P2                                             02/11/2016
!***********************************************************************
!
!brief    COMPUTES CONTRIBUTION TO MOMENTUM FORCES AND WATER LEVEL
!+        TERMS RESULTING FROM ICE PROCESSES.
!
!history  F. SOUILLE (EDF)
!+        30/10/2019
!+        V8P0
!+        Updated friction source term
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AT         |-->| TIME IN SECONDS
!| CHESTR     |-->| FRICTION COEFFICIENT
!| DT         |-->| TIME STEP
!| FU         |<->| SOURCE TERMS ON VELOCITY U
!| FV         |<->| SOURCE TERMS ON VELOCITY V
!| GRAV       |-->| GRAVITY
!| H          |-->| WATER DEPTH
!| KARMAN     |-->| VON KARMAN'S CONSTANT
!| NPOIN      |-->| NUMBER OF NODES IN THE MESH
!| T1         |<->| WORKING ARRAY
!| T2         |<->| WORKING ARRAY
!| T3         |<->| WORKING ARRAY
!| U          |-->| X COMPONENT OF THE VELOCITY
!| V          |-->| Y COMPONENT OF THE VELOCITY
!| WINDX      |<->| FIRST COMPONENT OF WIND VELOCITY
!| WINDY      |<->| SECOND COMPONENT OF WIND VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_special
      USE declarations_khione, ONLY: icover_impact,bd_ice,rho_air,fice,
     &                               fice_max,vz,ifrot,ifice,icestr,
     &                               thie,thifems,thifemf,hun,ro0
!
      USE meteo_telemac,       ONLY: sync_meteo,windx,windy
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,          INTENT(IN)    :: NPOIN
      TYPE(bief_obj),   INTENT(INOUT) :: T1,T2,T3
      DOUBLE PRECISION, INTENT(IN)    :: DT,AT,GRAV,KARMAN
      TYPE(bief_obj),   INTENT(INOUT) :: FU,FV
      TYPE(bief_obj),   INTENT(IN)    :: CHESTR, H,U,V
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                     :: I
      DOUBLE PRECISION            :: SP_EAU
      DOUBLE PRECISION            :: VMAG,WMAG,CD,CSTAR,CWSTAR,VZ31
      DOUBLE PRECISION            :: VZ32,UNTIER,UNSIX
      DOUBLE PRECISION, PARAMETER :: EPS=1.d-3
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!=======================================================================
!
!     FOR GOOD MEASURE - SHOULD BE DONE AT THE TELEMAC-2D LEVEL
!
!-----------------------------------------------------------------------
!
!     NOTE THAT "AT" FROM PROSOU IS AREADY TOO FAR GONE
      CALL sync_meteo(at-dt)
!
!
!=======================================================================
!
!     STATIC BORDER ICE GROWTH
!
!-----------------------------------------------------------------------
!
      IF( bd_ice ) THEN
!
!-----------------------------------------------------------------------
!       PREPARATION TO STATIC BORDER ICE GROWTH
!
        cstar  = 0.25d0
        cwstar = 1.d0
        cd     = 1.3d-3
        untier = 1.d0/3.d0
        unsix  = 1.d0/6.d0
!
        DO i = 1,npoin
!
! ~~>     WIND SPEED EFFECTS ON ICE
          wmag = sqrt( windx%R(i)**2 + windy%R(i)**2 )
! ~~>     FLOW SPEED EFFECTS ON ICE
          vmag = sqrt( u%R(i)**2 + v%R(i)**2 )
!
! ~~>     VERTICAL TURBULENT INTENSITY /!\ REPLACE BY CF%R
          vz31 = ( sqrt( cstar*grav ) * vmag*chestr%R(i)
     &             / max( h%R(i), 0.05d0 )**unsix)**3
          vz32 = cwstar * (cd*rho_air/ro0)**(1.5d0)*wmag**3
          vz%R(i) = (vz31 + vz32)**untier
!
        ENDDO
!
      ENDIF
!
!
!=======================================================================
!
!     THERMAL BALANCE AND ICE COVER GROWTH
!
!     FOLLOW-UP IMPACT ON THE HYDRODYNAMICS
!
!     TODO: Implement THERMAL_BUDGET=TRUE
!
!=======================================================================
!
!     ICE COVER IMPACT ON HYDRODYNAMICS
!
!-----------------------------------------------------------------------
!
      IF( icover_impact ) THEN
!
! ~~>   PREPARATORY WORK
!
!       COMPUTATION OF TOTAL ICE THICKNESS
        CALL os( 'X=Y+Z   ', x=t2,y=thifems,z=thifemf )
        CALL os( 'X=X+Y   ', x=t2,y=hun)
!
!       COMPUTATION OF SURFACE ICE GRADIENTS (TODO)
!        CALL OS( 'X=Y+Z   ', X=T1,Y=H,Z=ZF)
!        CALL OS( 'X=Y+CZ  ', X=T1,Y=T1,Z=T2,C=RHO_ICE/RO0 )
!        CALL VECTOR(DCOVX,'=','GRADF          X',U%ELM,1.D0,T1,
!     &    S,S,S,S,S,MESH,MSK,T3)
!        CALL VECTOR(DCOVY,'=','GRADF          Y',U%ELM,1.D0,T1,
!     &    S,S,S,S,S,MESH,MSK,T3)
!        IF(NCSIZE.GT.1) THEN
!          CALL PARCOM(DCOVX,2,MESH)
!          CALL PARCOM(DCOVY,2,MESH)
!        ENDIF
!        CALL OS( 'X=XY    ',X=DCOVX,Y=UNSV2D )
!        CALL OS( 'X=XY    ',X=DCOVY,Y=UNSV2D )
!
!       UPDATE FRICTION IF LINEAR FRICTION COEF LAW IS SELECTED
        IF (ifice .EQ. 1) THEN
          IF(ifrot.EQ.4) THEN
            CALL os( 'X=CY    ', x=icestr, y=t2, c=( fice/thie ) )
            CALL os( 'X=+(Y,C)', x=icestr, y=icestr, c=fice )
            CALL os( 'X=-(Y,C)', x=icestr, y=icestr, c=fice_max )
          ELSE
            WRITE(lu,*) 'IF LINEAR FRICTION IS SELECTED'
            WRITE(lu,*) 'PLEASE SELECT MANNING LAW'
            CALL plante(1)
            stop
          ENDIF
        ENDIF
!
        CALL friction_khione(npoin,ifrot,grav,karman,icestr,t1,h,u,v)
!
!       COMPUTATION WATER ELEVATION
        CALL os( 'X=+(Y,C)', x=t3, y=h, c=eps )
!
! ~~>   EFFECT OF THE ICE COVER ON THE HYDRODYNAMICS
        DO i = 1,npoin
!         T1: ICE FRICTION COEFFICIENT
!         T2: TOTAL ICE THICKNESS
!         T3: WATER DEPTH
!
          IF(t2%R(i).GT.eps) THEN
            sp_eau = sqrt( u%R(i)**2 + v%R(i)**2 )
            IF( h%R(i).LT.eps ) sp_eau =
     &        max( sp_eau, sqrt(grav*(eps-h%R(i))*h%R(i)/eps) )
!
! ~~>       UPDATE SOURCE TERM WITH ICE COVER FRICTION
            fu%R(i) = fu%R(i) - 0.5d0 * t1%R(i)*sp_eau*u%R(i) / t3%R(i)
            fv%R(i) = fv%R(i) - 0.5d0 * t1%R(i)*sp_eau*v%R(i) / t3%R(i)
!
! ~~>       LOCAL STATIC PRESSURE INCREASE DUE TO ICE THICKNESS (TODO)
!            IF( H%R(I).GT.EPS ) THEN
!              PATMOS%R(I) = PATMOS%R(I) + GRAV * RHO_ICE * T2%R(I)
!            ENDIF
          ENDIF
        ENDDO
!
! ~~>   SEEPAGE (TODO)
!
      ENDIF
!
!=======================================================================
!
!     CLOGGING
!
!     NO IMPACT ON THE HYDRODYNAMICS
!
!     TODO: Implement CLOGGING=TRUE
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END