The TELEMAC-MASCARET system  trunk
border_icover.f
Go to the documentation of this file.
1 ! ************************
2  SUBROUTINE border_icover
3 ! ************************
4 !
5  &( u,v, mesh )
6 !
7 !***********************************************************************
8 ! KHIONE V8P0
9 !***********************************************************************
10 !
11 !brief COMPUTES PRESENCE OF STATIC BORDER ICE
12 !
13 !+ Use of ICETYPE:
14 !+ 1 - open water
15 !+ 2 - static border ice (within or at the edge of the cover)
16 !+ 3 - border ice (both static and dynamic, and including the edge)
17 !+ Use of ICELOC:
18 !+ 1 - open water
19 !+ 2 - domain boundary node
20 !+ 3 - edge of the border ice where dynamic border ice accumulates
21 !+ Note that dynamic border ice (where ice is being accumulated),
22 !+ TODO: This method still need to be updated to account for the
23 !+ melting down of border ice (in patches if possible)
24 !
25 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
26 !| U |-->| X-COMPONENT OF THE VELOCITY
27 !| V |-->| Y-COMPONENT OF THE VELOCITY
28 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29 !
30  USE bief
32  & anfem0,it1,it2,t1,
33  & icetype,vcrbor,vcrbom,tcr,vz,
34  & tc,tmelt,iceloc,icetypep,t2
35  IMPLICIT NONE
36 !
37 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
38 !
39  TYPE(bief_mesh), INTENT(INOUT) :: MESH
40  TYPE(bief_obj), INTENT(IN) :: U,V
41 !
42 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
43 !
44  INTEGER I,J,K,N1,N2,IT,NPOIN,NELEM,NELMAX,IL,ITP
45  DOUBLE PRECISION THICK, USTAR,VMAG,VB
46 !
47 !-----------------------------------------------------------------------
48 !
49  npoin = mesh%NPOIN
50  nelem = mesh%NELEM
51  nelmax = mesh%NELMAX
52 !
53 !
54 !=======================================================================
55 !
56 ! 7 - STATIC BORDER ICE GROWTH
57 !
58 !-----------------------------------------------------------------------
59 !
60  IF( bd_ice ) THEN
61 !
62 !-----------------------------------------------------------------------
63 !
64 ! MAYBE THIS SHOULD BE DONE ONCE /!\
65 ! TODO: T1 CAN BE REPLACED BY VOLUPAR (ALREADY COMPUTED)
66 !
67  CALL vector(t1,'=','MASBAS ',u%ELM,1.d0,
68  & t2,t2,t2,t2,t2,t2,mesh,.false.,t2)
69 ! ASSEMBLING THE SUM OF ALL VALUES
70  IF( ncsize.GT.1 ) CALL parcom( t1,2,mesh )
71 !
72 !-----------------------------------------------------------------------
73 !
74 ! LOOKING FOR NEW BORDER ICE NODES (TO BE SWITCHED))
75 !
76  DO i = 1,npoin
77  it = icetype%I(i)
78  il = iceloc%I(i)
79  itp = icetypep%I(i)
80 !
81 ! CURRENT VELOCITY
82  vmag = sqrt( u%R(i)**2 + v%R(i)**2 )
83 ! BUOYANT VELOCITY
84  vb = max( 0.d0, -0.025d0*tcr%R(i) + 0.005d0 )
85 ! BORDER ICE THICKNESS /!\ TODO: DOUBLE CHECK VALIDITY
86  thick = thifems%R(i)+thifemf%R(i)
87 ! VELOCITY CRITERIA FOR DYNAMIC BORDER ICE GROWTH
88  ustar = max( 0.175d0, vmag / max(vcrbom,1.d-12) )
89 !
90 ! IT1 IS USED AS A STATIC ICE MASK DEFINING NODES FOR WHICH
91 ! THE THRESHOLDS ARE TRUE (I.E. SUBJECT TO BE BODER ICE NODE).
92 ! - IT1 = 0 IF THRESHOLDS ARE --NOT-- TRUE
93 ! - IT1 = 1 IF THRESHOLDS ARE TRUE
94 ! - IT1 = 2 IF THRESHOLD ARE TRUE AND IS ALSO ON THE DOMAIN
95 ! BORDER IT IS NOTED THAT IF IT1 = 2, THEN THE NODE IS
96 ! --NOT-- BORDER ICE ALREADY, BUT ABOUT TO BE SWITCHED.
97  it1%I(i) = 0
98 !
99 ! > THRESHOLDS FOR STATIC BORDER ICE FORMATION
100  IF( tcr%R(i).LE.(tc-tmelt%R(i)) .AND. ! THERMAL PROPERTY
101  & vb.GT.1.1d0*vz%R(i) .AND. ! BUOYANT VS. TURBULENCE
102  & vmag.LT.vcrbor ) THEN ! CRITICAL VELOCITY
103  IF( it .NE. 2 ) THEN
104 ! NODE OUTSIDE STATIC BORDER ICE COVER (NOT MULTIPLIER OF 3)
105  it1%I(i) = 1
106 ! BORDER NODE ABOUT TO BE SWITCH TO STATIC BORDER ICE
107  IF( il.EQ.2.OR.il.EQ.4 ) it1%I(i) = 2
108  ENDIF
109  ENDIF
110 !
111 ! IT2 IS USED AS A DYNAMIC ICE MASK DEFINING NODES FOR WHICH
112 ! THE THRESHOLDS ARE TRUE (I.E. SUBJECT TO BE BODER ICE NODE).
113 ! - IT2 = 0 IF THRESHOLDS ARE --NOT-- TRUE
114 ! - IT2 = 1 IF THRESHOLDS ARE TRUE
115 ! - IT2 = 2 IF THRESHOLD ARE TRUE AND IS ALSO ON THE DOMAIN
116 ! BORDER IT IS NOTED THAT IF IT2 = 2, THEN THE NODE IS
117 ! --NOT-- BORDER ICE ALREADY, BUT ABOUT TO BE SWITCHED.
118  it2%I(i) = 0
119 !
120 ! > THRESHOLDS FOR DYNAMIC BORDER ICE FORMATION
121  IF( !THICK.GE.0.0D0 .AND. ! CRITICAL THICKNESS
122  & ustar.LE.1.d0 ) THEN ! CRITICAL VELOCITY
123  IF( il .LT. 3.AND.itp.NE.2 ) THEN
124  it2%I(i) = 1
125 ! BORDER NODE ABOUT TO BE SWITCH TO STATIC BORDER ICE
126  IF( il.EQ.2 ) it2%I(i) = 2
127 !
128 ! ELSE !IF( IL.GE.3 ) THEN
129 ! > FURTHER GROWING DYNAMIC ICE
130 ! R = 14.1 * USTAR ** (-0.93) * ANFEM%R(I) ** 1.08
131 ! DPHI = LIN_WATAIR*( TWAT%R(I) - TAIR%R(I) )
132 ! DW = DPHI * R / ( ROEAU * LH_ICE )
133 ! DWB%R(I) = MIN( 1.D0, DWB%R(I) + DW / SQRT( T1%R(I) ) )
134 !
135  ENDIF
136  ENDIF
137 !
138  ENDDO
139 !
140 !-----------------------------------------------------------------------
141 !
142 ! CHECK BORDER NODES ABOUT TO BE SWITCHED INTO BORDER ICE NODES
143 !
144  DO j = 1,nelem
145 !
146  DO k = 1,3
147 ! RA: why we compute again node indices ? They are stored in IKLE (NELEM,3)
148 ! 1- introduce IKLE as an argument like IKLE(NELEM,3)
149 ! 2- nodes of element IELEM are :
150 ! I1 = IKLE(IELEM,1)
151 ! I2 = IKLE(IELEM,2)
152 ! I3 = IKLE(IELEM,3)
153 !
154  i = mesh%IKLE%I( j + (k-1)*nelmax )
155 !
156  it = icetype%I(i)
157  itp = icetypep%I(i)
158  il = iceloc%I(i)
159 ! NODES ABOUT TO BE SWITCHED NEED TO BE BELOW THRESHOLDS AND
160 ! --NOT-- ALREADY BODER ICE NODES. THEY ALSO NEED TO BE CLOSE
161 ! TO TWO OTHER (SUFFICIENTLY THICK) BORDER ICE NODES.
162  IF( it1%I(i).EQ.1 ) THEN
163  n1 = icetype%I( mesh%IKLE%I( j + mod(k,3)*nelmax ) )
164  n2 = icetype%I( mesh%IKLE%I( j + mod(k+1,3)*nelmax ) )
165  IF( n1 .EQ. 2 .AND. n2 .EQ. 2 ) THEN
166  it1%I(i) = 2
167  ENDIF
168  ENDIF
169  IF( it2%I(i).EQ.1 ) THEN
170  n1 = icetypep%I( mesh%IKLE%I( j + mod(k,3)*nelmax ) )
171  n2 = icetypep%I( mesh%IKLE%I( j + mod(k+1,3)*nelmax ) )
172  IF( n1 .EQ. 2 .AND. n2 .EQ. 2 ) THEN
173  it2%I(i) = 2
174  ENDIF
175  ENDIF
176  ENDDO
177 !
178  ENDDO
179 !
180 ! ASSEMBLING THE MAX OF ALL VALUES - THE TWOS WILL BE SWITCHED
181 ! CONVERSION BETWEEN INTEGER AND DOUBLE IS NECESSARY FOR PARCOM
182  IF(ncsize.GT.1)CALL parcom2i(it1%I,it2%I,it2%I,npoin,1,3,2,mesh)
183 !
184 !-----------------------------------------------------------------------
185 !
186 ! SWITCHES BORDER NODES INTO STATIC BORDER ICE NODES
187 !
188  DO i = 1,npoin
189  it = icetype%I(i)
190  itp = icetypep%I(i)
191  il = iceloc%I(i)
192 !
193  IF( it1%I(i).EQ.2 ) THEN
194 ! > STATIC BORDER ICE
195  icetype%I(i) = 2 !=> INSTANT SWITCH
196 !
197 ! SWITCH FROM SOLID DYNAMIC TO STATIC BORDER ICE
198  IF( itp .EQ. 2 ) THEN
199 !
200 ! SWITCH FROM GROWING DYNAMIC TO STATIC BORDER ICE
201  ELSEIF( il .GE. 3 ) THEN
202 ! USE DWB TO COMPUTE THIFEMS, THIFEMF AND ANFEM ?
203  IF(il.EQ.3) iceloc%I(i) = 1 !=> INSTANT SWITCH
204  IF(il.EQ.4) iceloc%I(i) = 2 !=> INSTANT SWITCH
205  icetypep%I(i) = 2
206 !
207 ! NEW STATIC BORDER ICE
208  ELSE
209 ! NOT PART OF THE DYNAMIC ICE COVER
210  icetypep%I(i) = 2 !=> INSTANT SWITCH
211  thifems%R(i) = 0.d0
212  thifemf%R(i) = 0.d0
213  anfem%R(i) = anfem0
214 !
215  ENDIF
216 !
217 ! > DYNAMIC BORDER ICE
218  ELSEIF( it2%I(i).EQ.2 ) THEN
219  IF( il .LT. 3 ) THEN
220  thifems%R(i) = 0.d0
221  thifemf%R(i) = 0.d0
222  anfem%R(i) = anfem0
223  IF(il.EQ.1) iceloc%I(i) = 3 !=> INSTANT SWITCH
224  IF(il.EQ.2) iceloc%I(i) = 4 !=> INSTANT SWITCH
225  ENDIF
226  ENDIF
227 !
228  ENDDO
229 !
230  ENDIF
231 !
232 !
233 !=======================================================================
234 !
235 !-----------------------------------------------------------------------
236 !
237  RETURN
238  END SUBROUTINE border_icover
subroutine parcom2i(X1, X2, X3, NPOIN, NPLAN, ICOM, IAN, MESH)
Definition: parcom2i.f:7
subroutine border_icover(U, V, MESH)
Definition: border_icover.f:7
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
type(bief_obj), target anfem
type(bief_obj), target thifemf
type(bief_obj), target thifems
Definition: bief.f:3