The TELEMAC-MASCARET system  trunk
source_icover.f
Go to the documentation of this file.
1 ! ************************
2  SUBROUTINE source_icover
3 ! ************************
4  &(npoin,fu,fv,h,u,v,t1,t2,t3,grav,karman,chestr,dt,at)
5 !
6 !***********************************************************************
7 ! KHIONE V7P2 02/11/2016
8 !***********************************************************************
9 !
10 !brief COMPUTES CONTRIBUTION TO MOMENTUM FORCES AND WATER LEVEL
11 !+ TERMS RESULTING FROM ICE PROCESSES.
12 !
13 !history F. SOUILLE (EDF)
14 !+ 30/10/2019
15 !+ V8P0
16 !+ Updated friction source term
17 !
18 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
19 !| AT |-->| TIME IN SECONDS
20 !| CHESTR |-->| FRICTION COEFFICIENT
21 !| DT |-->| TIME STEP
22 !| FU |<->| SOURCE TERMS ON VELOCITY U
23 !| FV |<->| SOURCE TERMS ON VELOCITY V
24 !| GRAV |-->| GRAVITY
25 !| H |-->| WATER DEPTH
26 !| KARMAN |-->| VON KARMAN'S CONSTANT
27 !| NPOIN |-->| NUMBER OF NODES IN THE MESH
28 !| T1 |<->| WORKING ARRAY
29 !| T2 |<->| WORKING ARRAY
30 !| T3 |<->| WORKING ARRAY
31 !| U |-->| X COMPONENT OF THE VELOCITY
32 !| V |-->| Y COMPONENT OF THE VELOCITY
33 !| WINDX |<->| FIRST COMPONENT OF WIND VELOCITY
34 !| WINDY |<->| SECOND COMPONENT OF WIND VELOCITY
35 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36 !
37  USE bief
40  & fice_max,vz,ifrot,ifice,icestr,
41  & thie,thifems,thifemf,hun,ro0
42 !
44 !
45  IMPLICIT NONE
46 !
47 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
48 !
49  INTEGER, INTENT(IN) :: NPOIN
50  TYPE(bief_obj), INTENT(INOUT) :: T1,T2,T3
51  DOUBLE PRECISION, INTENT(IN) :: DT,AT,GRAV,KARMAN
52  TYPE(bief_obj), INTENT(INOUT) :: FU,FV
53  TYPE(bief_obj), INTENT(IN) :: CHESTR, H,U,V
54 !
55 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
56 !
57  INTEGER :: I
58  DOUBLE PRECISION :: SP_EAU
59  DOUBLE PRECISION :: VMAG,WMAG,CD,CSTAR,CWSTAR,VZ31
60  DOUBLE PRECISION :: VZ32,UNTIER,UNSIX
61  DOUBLE PRECISION, PARAMETER :: EPS=1.d-3
62 !
63 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
64 !
65 !=======================================================================
66 !
67 ! FOR GOOD MEASURE - SHOULD BE DONE AT THE TELEMAC-2D LEVEL
68 !
69 !-----------------------------------------------------------------------
70 !
71 ! NOTE THAT "AT" FROM PROSOU IS AREADY TOO FAR GONE
72  CALL sync_meteo(at-dt)
73 !
74 !
75 !=======================================================================
76 !
77 ! STATIC BORDER ICE GROWTH
78 !
79 !-----------------------------------------------------------------------
80 !
81  IF( bd_ice ) THEN
82 !
83 !-----------------------------------------------------------------------
84 ! PREPARATION TO STATIC BORDER ICE GROWTH
85 !
86  cstar = 0.25d0
87  cwstar = 1.d0
88  cd = 1.3d-3
89  untier = 1.d0/3.d0
90  unsix = 1.d0/6.d0
91 !
92  DO i = 1,npoin
93 !
94 ! ~~> WIND SPEED EFFECTS ON ICE
95  wmag = sqrt( windx%R(i)**2 + windy%R(i)**2 )
96 ! ~~> FLOW SPEED EFFECTS ON ICE
97  vmag = sqrt( u%R(i)**2 + v%R(i)**2 )
98 !
99 ! ~~> VERTICAL TURBULENT INTENSITY /!\ REPLACE BY CF%R
100  vz31 = ( sqrt( cstar*grav ) * vmag*chestr%R(i)
101  & / max( h%R(i), 0.05d0 )**unsix)**3
102  vz32 = cwstar * (cd*rho_air/ro0)**(1.5d0)*wmag**3
103  vz%R(i) = (vz31 + vz32)**untier
104 !
105  ENDDO
106 !
107  ENDIF
108 !
109 !
110 !=======================================================================
111 !
112 ! THERMAL BALANCE AND ICE COVER GROWTH
113 !
114 ! FOLLOW-UP IMPACT ON THE HYDRODYNAMICS
115 !
116 ! TODO: Implement THERMAL_BUDGET=TRUE
117 !
118 !=======================================================================
119 !
120 ! ICE COVER IMPACT ON HYDRODYNAMICS
121 !
122 !-----------------------------------------------------------------------
123 !
124  IF( icover_impact ) THEN
125 !
126 ! ~~> PREPARATORY WORK
127 !
128 ! COMPUTATION OF TOTAL ICE THICKNESS
129  CALL os( 'X=Y+Z ', x=t2,y=thifems,z=thifemf )
130  CALL os( 'X=X+Y ', x=t2,y=hun)
131 !
132 ! COMPUTATION OF SURFACE ICE GRADIENTS (TODO)
133 ! CALL OS( 'X=Y+Z ', X=T1,Y=H,Z=ZF)
134 ! CALL OS( 'X=Y+CZ ', X=T1,Y=T1,Z=T2,C=RHO_ICE/RO0 )
135 ! CALL VECTOR(DCOVX,'=','GRADF X',U%ELM,1.D0,T1,
136 ! & S,S,S,S,S,MESH,MSK,T3)
137 ! CALL VECTOR(DCOVY,'=','GRADF Y',U%ELM,1.D0,T1,
138 ! & S,S,S,S,S,MESH,MSK,T3)
139 ! IF(NCSIZE.GT.1) THEN
140 ! CALL PARCOM(DCOVX,2,MESH)
141 ! CALL PARCOM(DCOVY,2,MESH)
142 ! ENDIF
143 ! CALL OS( 'X=XY ',X=DCOVX,Y=UNSV2D )
144 ! CALL OS( 'X=XY ',X=DCOVY,Y=UNSV2D )
145 !
146 ! UPDATE FRICTION IF LINEAR FRICTION COEF LAW IS SELECTED
147  IF (ifice .EQ. 1) THEN
148  IF(ifrot.EQ.4) THEN
149  CALL os( 'X=CY ', x=icestr, y=t2, c=( fice/thie ) )
150  CALL os( 'X=+(Y,C)', x=icestr, y=icestr, c=fice )
151  CALL os( 'X=-(Y,C)', x=icestr, y=icestr, c=fice_max )
152  ELSE
153  WRITE(lu,*) 'IF LINEAR FRICTION IS SELECTED'
154  WRITE(lu,*) 'PLEASE SELECT MANNING LAW'
155  CALL plante(1)
156  stop
157  ENDIF
158  ENDIF
159 !
160  CALL friction_khione(npoin,ifrot,grav,karman,icestr,t1,h,u,v)
161 !
162 ! COMPUTATION WATER ELEVATION
163  CALL os( 'X=+(Y,C)', x=t3, y=h, c=eps )
164 !
165 ! ~~> EFFECT OF THE ICE COVER ON THE HYDRODYNAMICS
166  DO i = 1,npoin
167 ! T1: ICE FRICTION COEFFICIENT
168 ! T2: TOTAL ICE THICKNESS
169 ! T3: WATER DEPTH
170 !
171  IF(t2%R(i).GT.eps) THEN
172  sp_eau = sqrt( u%R(i)**2 + v%R(i)**2 )
173  IF( h%R(i).LT.eps ) sp_eau =
174  & max( sp_eau, sqrt(grav*(eps-h%R(i))*h%R(i)/eps) )
175 !
176 ! ~~> UPDATE SOURCE TERM WITH ICE COVER FRICTION
177  fu%R(i) = fu%R(i) - 0.5d0 * t1%R(i)*sp_eau*u%R(i) / t3%R(i)
178  fv%R(i) = fv%R(i) - 0.5d0 * t1%R(i)*sp_eau*v%R(i) / t3%R(i)
179 !
180 ! ~~> LOCAL STATIC PRESSURE INCREASE DUE TO ICE THICKNESS (TODO)
181 ! IF( H%R(I).GT.EPS ) THEN
182 ! PATMOS%R(I) = PATMOS%R(I) + GRAV * RHO_ICE * T2%R(I)
183 ! ENDIF
184  ENDIF
185  ENDDO
186 !
187 ! ~~> SEEPAGE (TODO)
188 !
189  ENDIF
190 !
191 !=======================================================================
192 !
193 ! CLOGGING
194 !
195 ! NO IMPACT ON THE HYDRODYNAMICS
196 !
197 ! TODO: Implement CLOGGING=TRUE
198 !
199 !-----------------------------------------------------------------------
200 !
201  RETURN
202 !
203 !-----------------------------------------------------------------------
204 !
205  END
type(bief_obj), target, public windy
subroutine friction_khione(NPOIN, KFROT, GRAV, KARMAN, CHESTR, CF, H, U, V)
type(bief_obj), target, public windx
double precision rho_air
subroutine, public sync_meteo(WHEN)
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine source_icover(NPOIN, FU, FV, H, U, V, T1, T2, T3, GRAV, KARMAN, CHESTR, DT, AT)
Definition: source_icover.f:6
Definition: bief.f:3