The TELEMAC-MASCARET system  trunk
borice.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE borice
3 ! *****************
4 !
5  &(h,u,v,f,at,lt,dt,tra05,tra06,liubor,nptfr,numliq,klog,mask,mesh
6  & ,s)
7 !
8 !***********************************************************************
9 ! KHIONE V7P3
10 !***********************************************************************
11 !
12 !brief INCLUDE ICE PROCESSES TO THE BOUNDARIES SETTINGS.
13 !
14 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
15 !| AT |-->| TIME OF THE DATASET
16 !| DT |-->| TIME STEP
17 !| F |-->| TRACER VALUES
18 !| H |-->| WATER DEPTH
19 !| KLOG |-->| CONVENTION FOR SOLID BOUNDARY
20 !| LIUBOR |-->| TYPE OF BOUNDARY CONDITION FOR THE VELOCITY
21 !| LT |-->| NUMBER OF TIME STEPS
22 !| MASK |-->| BLOCK OF MASKS FOR DIFFERENT BOUNDARY CONDITIONS
23 !| MESH |-->| MESH STRUCTURE
24 !| NPTFR |-->| NUMBER OF BOUNDARY POINT
25 !| NUMLIQ |-->| LIQUID BOUNDARY NUMBER
26 !| S |-->| VOID STRUCTURE
27 !| TRA05 |-->| WORK ARRAY IN A BIEF_OBJ STRUCTURE
28 !| TRA06 |-->| WORK ARRAY IN A BIEF_OBJ STRUCTURE
29 !| U |-->| X-COMPONENT OF THE VELOCITY
30 !| V |-->| Y-COMPONENT OF THE VELOCITY
31 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
32 !
33  USE bief
34  USE declarations_khione, ONLY:
44  USE interface_parallel, ONLY : p_dsum
45  IMPLICIT NONE
46 !
47 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
48 !
49  INTEGER, INTENT(IN) :: KLOG, NPTFR
50 !
51  INTEGER, INTENT(IN) :: NUMLIQ(nptfr)
52  INTEGER, INTENT(IN) :: LIUBOR(nptfr)
53  TYPE(bief_mesh), INTENT(INOUT) :: MESH
54  TYPE(bief_obj), INTENT(IN) :: H,U,V,F,S
55  TYPE(bief_obj), INTENT(INOUT) :: MASK, TRA05,TRA06
56  INTEGER, INTENT(IN) :: LT
57  DOUBLE PRECISION, INTENT(IN) :: AT,DT
58 !
59 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
60 !
61  INTEGER K,N,I,IFRLIQ,IELEB,IFILE
62  LOGICAL OUTPUT
63  DOUBLE PRECISION PI, ANG,RV0,RV1,RT0,RT1, RATE, QCV,CV,CA,
64  & mfbar,mfvb,mftb, mfgrd,mfvt,mftt,mass0, dw,dv,taw, nv,nt
65  DOUBLE PRECISION, ALLOCATABLE :: QF(:),QA(:),QV(:)
66  CHARACTER(LEN=250) CELL
67 !
68 !-----------------------------------------------------------------------
69 !
70  INTRINSIC adjustl,trim
71 !
72 !=======================================================================
73 !
74 ! 5 - CLOGGING
75 !
76 ! NO IMPACT ON THE HYDRODYNAMICS FOR NOW
77 ! WARNING: ONLY WORKS FOR MONO-CLASS FRAZIL ICE
78 !
79  IF( clogging ) THEN
80 !
81  IF(nseclog.GT.0) THEN
82  ALLOCATE(qf(nseclog))
83  ALLOCATE(qa(nseclog))
84  ALLOCATE(qv(nseclog))
85  CALL section_khione(mesh,mesh%IKLE%I,mesh%NELMAX,
86  & mesh%IFABOR%I,f,u,v,qf,qa,qv,h,s,lt)
87  ENDIF
88 !
89  IF( lt.LE.1 .AND. ntotclog.GT.0 ) THEN
90 !
91  IF(nfrclog.GT.0) THEN
92 !
93 !-----------------------------------------------------------------------
94 ! COMPUTES THE LENGTH OF THE OPEN BOUNDARY (CLOG_TLGTH) ONCE
95  CALL os('X=C ', x=t3, c=1.d0)
96  DO ifrliq = 1,nfrclog
97  n = numclog(ifrliq)
98  CALL os( 'X=0 ', x=tra05 )
99  DO ieleb = 1,mesh%NELEB
100  k = mesh%IKLBOR%I(ieleb)
101  IF( numliq(k).EQ.n )
102  & tra05%R(ieleb) = mask%ADR(8)%P%R(ieleb)
103  ENDDO
104  CALL vector(tra06,'=','MASVEC ',
105  & ielbor(11,1),1.d0,t3,t3,t3,t3,t3,t3,mesh,.true.,tra05)
106  clog_tlgth(ifrliq) = bief_sum(tra06)
107  IF( ncsize.GT.1 )
108  & clog_tlgth(ifrliq) = p_dsum(clog_tlgth(ifrliq))
109  ENDDO
110  ENDIF
111 !
112 !-----------------------------------------------------------------------
113 ! HEADER OF ASCII CLOGGING FILE
114 !
115  IF( ice_files(clgrfo)%NAME.NE.' ' ) THEN
116  IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
117  ifile = ice_files(clgrfo)%LU
118  WRITE(ifile,301) ''
119  WRITE(ifile,301) titicecas
120  WRITE(ifile,301) ''
121  301 FORMAT('# ',a)
122  lines%HTXT = 'T'
123  WRITE(ifile,300,advance='NO') trim(adjustl(lines%HTXT))
124  lines%HUNT = 'S'
125  DO ifrliq = 1,ntotclog
126  IF( lines%CELLS(ifrliq)%NVAL.EQ.0 ) THEN
127  lines%CELLS(ifrliq)%NVAL = 8
128  ALLOCATE(lines%CELLS(ifrliq)%TXT(8))
129  ALLOCATE(lines%CELLS(ifrliq)%UNT(8))
130  ALLOCATE(lines%CELLS(ifrliq)%VAL(8))
131  ENDIF
132  WRITE(cell,'(I12)') ifrliq
133  lines%CELLS(ifrliq)%TXT(1) = 'FRAZIL('
134  & // trim(adjustl(cell)) // ')'
135  lines%CELLS(ifrliq)%UNT(1) = 'SI'
136  lines%CELLS(ifrliq)%TXT(2) = 'WIDTH_V('
137  & // trim(adjustl(cell)) // ')'
138  lines%CELLS(ifrliq)%UNT(2) = 'M'
139  lines%CELLS(ifrliq)%TXT(3) = 'WIDTH_H('
140  & // trim(adjustl(cell)) // ')'
141  lines%CELLS(ifrliq)%UNT(3) = 'M'
142  lines%CELLS(ifrliq)%TXT(4) = 'CLG_AREA('
143  & // trim(adjustl(cell)) // ')'
144  lines%CELLS(ifrliq)%UNT(4) = 'M^2'
145  lines%CELLS(ifrliq)%TXT(5) = 'IMASS_V('
146  & // trim(adjustl(cell)) // ')'
147  lines%CELLS(ifrliq)%UNT(5) = 'G'
148  lines%CELLS(ifrliq)%TXT(6) = 'IMASS_H('
149  & // trim(adjustl(cell)) // ')'
150  lines%CELLS(ifrliq)%UNT(6) = 'G'
151  lines%CELLS(ifrliq)%TXT(7) = 'CLG_VOL('
152  & // trim(adjustl(cell)) // ')'
153  lines%CELLS(ifrliq)%UNT(7) = 'M^3'
154  lines%CELLS(ifrliq)%TXT(8) = 'CLG_MAS('
155  & // trim(adjustl(cell)) // ')'
156  lines%CELLS(ifrliq)%UNT(8) = 'SI'
157  WRITE(ifile,302,advance='NO') ( ' ' //
158  & trim(adjustl(lines%CELLS(ifrliq)%TXT(k))), k=1,8)
159  ENDDO
160  WRITE(ifile,300) ''
161  WRITE(ifile,300,advance='NO') trim(adjustl(lines%HUNT))
162  DO i = 1,ntotclog
163  WRITE(ifile,302,advance='NO')
164  & ( ' ' // trim(adjustl(lines%CELLS(i)%UNT(k)))
165  & , k=1,8 )
166  ENDDO
167  WRITE(ifile,300) ''
168  300 FORMAT(a)
169  302 FORMAT(8a)
170  ENDIF
171  ENDIF
172  ENDIF
173 !
174 !-----------------------------------------------------------------------
175 ! PREPARE COMPUTATION OF THE FLUXES AND THE AVERAGE BAR HEIGHTS
176 !
177 !
178  CALL os('X=Y ', x=t2, y=h )
179  CALL os('X=YZ ', x=t1, y=h, z=f%ADR(ind_fra)%P )
180 
181 !
182  pi = 4.d0*atan(1.d0)
183  output = .false.
184 !
185 ! ACUMULATION ANGLE
186  ang = pi * clog_theta / 180.d0
187 ! INITIAL FRAZIL RADIUS
188  rv0 = clog_vdiam*cos(ang)
189  rt0 = clog_tdiam*cos(ang)
190 !
191 ! CLOGGING RATES ( AF = 1.D0 )
192  rate = 2.d0*sin(ang)*sin(ang)*(1.d0)/( 1.d0-clog_ef )/ang
193 !
194 !-----------------------------------------------------------------------
195 !
196  IF(nfrclog.GT.0) THEN
197  DO ifrliq = 1,nfrclog
198  n = numclog(ifrliq)
199 !
200 ! MASKING ONE LIQUID BOUNDARY AT A TIME
201  CALL os( 'X=0 ', x=tra05 )
202  DO ieleb = 1,mesh%NELEB
203  k = mesh%IKLBOR%I(ieleb)
204  IF( numliq(k).EQ.n )
205  & tra05%R(ieleb) = mask%ADR(8)%P%R(ieleb)
206  ENDDO
207 !
208 ! FRAZIL FLUX = INTEGRAL( F * H * UV ) ACROSS BOUNDARY N
209  CALL vector(tra06,'=','FLUBDF ',
210  & ielbor(11,1),1.d0,t1,t1,t1,u,v,v,mesh,.true.,tra05)
211  qcv = bief_sum(tra06)
212  IF( ncsize.GT.1 ) qcv = p_dsum(qcv)
213 ! FRAZIL FLUX BY UNIT OF WIDTH
214  qcv = abs(qcv) / clog_tlgth(ifrliq)
215 !
216 ! CROSS SECTIONAL AREA = INTEGRAL( H ) ACROSS BOUNDARY N
217  CALL vector(tra06,'=','MASVEC ',
218  & ielbor(11,1),1.d0,t2,t2,t2,t2,t2,t2,mesh,.true.,tra05)
219  ca = bief_sum(tra06)
220  IF( ncsize.GT.1 ) ca = p_dsum(ca)
221 ! AVERAGE BAR HEIGHT
222  IF( clog_tlgth(ifrliq).GT.0.d0 ) THEN
223  clog_vlgth(ifrliq) = ca/clog_tlgth(ifrliq)
224  ELSE
225  clog_vlgth(ifrliq) = 0.d0
226  ENDIF
227 !
228 ! AVERAGE CONCENTRATION = INTEGRAL( F * H ) ACROSS BOUNDARY N
229  CALL vector(tra06,'=','MASVEC ',
230  & ielbor(11,1),1.d0,t1,t1,t1,t1,t1,t1,mesh,.true.,tra05)
231  cv = bief_sum(tra06)
232  IF( ncsize.GT.1 ) cv = p_dsum(cv)
233 ! ... DIVIDED BY THE AREA OF THE CROSS SECTION
234  IF( ca.GT.0.d0 ) THEN
235  cv = cv/ca
236  ELSE
237  cv = 0.d0
238  ENDIF
239 !
240 ! REPRESENTATIVE NUMBERS OF BARS
241  IF( clog_vdist*clog_vdiam.GT.0.d0 ) THEN
242 ! THIS COULD BE AN ODD NUMBER
243  nv = clog_tlgth(ifrliq)/clog_vdist
244  ELSE
245  nv = 0.d0
246  ENDIF
247  IF( clog_tdist*clog_tdiam.GT.0.d0 ) THEN
248 ! COUNTIN ONLY THE FULLY SUBMERGED
249  nt = 1.d0*int( clog_vlgth(ifrliq)/clog_tdist )
250  ELSE
251  nt = 0.d0
252  ENDIF
253 !
254 ! VERTICAL BARS
255  IF( clog_vlgth(ifrliq)*nv.GT.0.d0 ) THEN
256  clog_vwdth(ifrliq) = clog_vwdth(ifrliq) +
257  & dt*rate*qcv/clog_vlgth(ifrliq)
258  rv1 = clog_vwdth(ifrliq)/2.d0/sin(ang)
259  CALL clogged_on_bar( rv0,rv1,
260  & clog_vdiam,clog_vlgth(ifrliq),nv,ang,mfvb,mfvt)
261  ELSE
262  rv1 = 0.d0
263  mfvb = 0.d0
264  mfvt = 0.d0
265  ENDIF
266 !
267 ! TRANSVERSE BARS
268  IF( clog_tlgth(ifrliq)*nt.GT.0.d0 ) THEN
269  clog_twdth(ifrliq) = clog_twdth(ifrliq) +
270  & dt*rate*qcv/clog_vlgth(ifrliq)
271  rt1 = clog_twdth(ifrliq)/2.d0/sin(ang)
272  CALL clogged_on_bar( rt0,rt1,
273  & clog_tdiam,clog_tlgth(ifrliq),nt,ang,mftb,mftt)
274  ELSE
275  rt1 = 0.d0
276  mftb = 0.d0
277  mftt = 0.d0
278  ENDIF
279 !
280  IF( ice_files(clgrfo)%NAME.EQ.' ' ) cycle
281  IF( lt.NE.leoprd*int(lt/leoprd) .AND. lt.GT.1 ) cycle
282 !
283 ! ADDITIONAL PRINTOUTS
284  output = .true.
285 !
286  IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
287  WRITE(lines%CELLS(ifrliq)%VAL(1),303) cv
288  WRITE(lines%CELLS(ifrliq)%VAL(2),303) clog_vwdth(ifrliq)
289  WRITE(lines%CELLS(ifrliq)%VAL(3),303) clog_twdth(ifrliq)
290 ! WIDTH OF THE VERTICAL BARS EXCLUDED FROM THE TRANSVERSE BARS
291  dw = clog_vdist - clog_vwdth(ifrliq)
292  dv = max( clog_vwdth(ifrliq),clog_vdiam )
293 ! AVAILABLE AREA FOR THE RACK
294  taw = clog_vlgth(ifrliq)*dw*(nv-1.d0)
295  & - clog_twdth(ifrliq) * ( clog_tlgth(ifrliq)-nv*dv ) * nt
296  taw = max(taw,0.d0)
297  WRITE(lines%CELLS(ifrliq)%VAL(4),303) taw
298 ! FRAZIL VOLUME FOR THE GRID
299  clog_volum(ifrliq) = clog_volum(ifrliq) +
300  & qcv * ( 1.d0-clog_ef ) * dt *
301  & ( clog_vwdth(ifrliq)*nv + clog_twdth(ifrliq)*nt )
302 ! FRAZIL MASS FOR THE GRID
303  mass0 = clog_volum(ifrliq)*rho_ice
304 !
305  mfbar = mfvb + mftb
306  WRITE(lines%CELLS(ifrliq)%VAL(5),303) mfbar
307  mfgrd = mfvt + mftt
308  WRITE(lines%CELLS(ifrliq)%VAL(6),303) mfgrd
309  WRITE(lines%CELLS(ifrliq)%VAL(7),303) clog_volum(ifrliq)
310  WRITE(lines%CELLS(ifrliq)%VAL(8),303) mass0
311  ENDIF
312 !
313  ENDDO
314  ENDIF
315 !
316  IF(nseclog.GT.0) THEN
317  DO i=nfrclog+1,ntotclog
318  clog_vlgth(i) = qa(i-nfrclog)/clog_tlgth(i)
319  qcv = qf(i-nfrclog)
320  IF(qa(i-nfrclog).GT.0.d0) THEN
321  cv = qv(i-nfrclog) / qa(i-nfrclog)
322  ELSE
323  cv=0.d0
324  ENDIF
325 ! REPRESENTATIVE NUMBERS OF BARS
326  IF( clog_vdist*clog_vdiam.GT.0.d0 ) THEN
327 ! THIS COULD BE AN ODD NUMBER
328  nv = clog_tlgth(i)/clog_vdist
329  ELSE
330  nv = 0.d0
331  ENDIF
332  IF( clog_tdist*clog_tdiam.GT.0.d0 ) THEN
333 ! COUNTIN ONLY THE FULLY SUBMERGED
334  nt = 1.d0*int( clog_vlgth(i)/clog_tdist )
335  ELSE
336  nt = 0.d0
337  ENDIF
338 !
339 ! VERTICAL BARS
340  IF( clog_vlgth(i)*nv.GT.0.d0 ) THEN
341  clog_vwdth(i) = clog_vwdth(i) +
342  & dt*rate*qcv/clog_vlgth(i)
343  rv1 = clog_vwdth(i)/2.d0/sin(ang)
344  CALL clogged_on_bar( rv0,rv1,
345  & clog_vdiam,clog_vlgth(i),nv,ang,mfvb,mfvt)
346  ELSE
347  rv1 = 0.d0
348  mfvb = 0.d0
349  mfvt = 0.d0
350  ENDIF
351 !
352 ! TRANSVERSE BARS
353  IF( clog_tlgth(i)*nt.GT.0.d0 ) THEN
354  clog_twdth(i) = clog_twdth(i) +
355  & dt*rate*qcv/clog_vlgth(i)
356  rt1 = clog_twdth(i)/2.d0/sin(ang)
357  CALL clogged_on_bar( rt0,rt1,
358  & clog_tdiam,clog_tlgth(i),nt,ang,mftb,mftt)
359  ELSE
360  rt1 = 0.d0
361  mftb = 0.d0
362  mftt = 0.d0
363  ENDIF
364 !
365  IF( ice_files(clgrfo)%NAME.EQ.' ' ) cycle
366  IF( lt.NE.leoprd*int(lt/leoprd) .AND. lt.GT.1 ) cycle
367 !
368 ! ADDITIONAL PRINTOUTS
369  output = .true.
370 !
371  IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
372  WRITE(lines%CELLS(i)%VAL(1),303) cv
373  WRITE(lines%CELLS(i)%VAL(2),303) clog_vwdth(i)
374  WRITE(lines%CELLS(i)%VAL(3),303) clog_twdth(i)
375 ! WIDTH OF THE VERTICAL BARS EXCLUDED FROM THE TRANSVERSE BARS
376  dw = clog_vdist - clog_vwdth(i)
377  dv = max( clog_vwdth(i),clog_vdiam )
378 ! AVAILABLE AREA FOR THE RACK
379  taw = clog_vlgth(i)*dw*(nv-1.d0)
380  & - clog_twdth(i) * ( clog_tlgth(i)-nv*dv ) * nt
381  taw = max(taw,0.d0)
382  WRITE(lines%CELLS(i)%VAL(4),303) taw
383 ! FRAZIL VOLUME FOR THE GRID
384  clog_volum(i) = clog_volum(i) +
385  & qcv * ( 1.d0-clog_ef ) * dt *
386  & ( clog_vwdth(i)*nv + clog_twdth(i)*nt )
387 ! FRAZIL MASS FOR THE GRID
388  mass0 = clog_volum(i)*rho_ice
389 !
390  mfbar = mfvb + mftb
391  WRITE(lines%CELLS(i)%VAL(5),303) mfbar
392  mfgrd = mfvt + mftt
393  WRITE(lines%CELLS(i)%VAL(6),303) mfgrd
394  WRITE(lines%CELLS(i)%VAL(7),303) clog_volum(i)
395  WRITE(lines%CELLS(i)%VAL(8),303) mass0
396  ENDIF
397 
398  ENDDO
399  ENDIF
400 !
401 !-----------------------------------------------------------------------
402 !
403  303 FORMAT(g16.9e3)
404  IF( ice_files(clgrfo)%NAME.NE.' ' .AND. output ) THEN
405  IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
406  WRITE(lines%HVAL,303) at
407  WRITE(ice_files(clgrfo)%LU,302,advance='NO')
408  & trim(adjustl(lines%HVAL))
409  DO i=1,ntotclog
410  WRITE(ice_files(clgrfo)%LU,302,advance='NO')
411  & ( ' ' // trim(adjustl(lines%CELLS(i)%VAL(k)))
412  & , k=1,8 )
413  ENDDO
414  WRITE(ice_files(clgrfo)%LU,300) ''
415  ENDIF
416  ENDIF
417 !
418  ENDIF
419 !
420 !=======================================================================
421 !
422 ! 7 - STATIC BORDER ICE TRACKING
423 !
424  IF( bd_ice ) THEN
425 !
426  IF( lt.LE.1 ) THEN
427 !
428  DO k = 1,nptfr
429  IF( liubor(k).EQ.klog ) THEN
430  n = mesh%NBOR%I(k)
431  IF( icetype%I(n).EQ.1 ) iceloc%I(n) = 2
432  ENDIF
433  ENDDO
434 !
435  ENDIF
436 !
437  ENDIF
438 !
439 !=======================================================================
440 !
441 ! ? - OTHER PROCESSES
442 !
443  IF(nseclog.GT.0) THEN
444  DEALLOCATE(qf)
445  DEALLOCATE(qa)
446  DEALLOCATE(qv)
447  ENDIF
448 !
449 !-----------------------------------------------------------------------
450 !
451  RETURN
452  END SUBROUTINE
double precision clog_vdist
double precision, dimension(:), allocatable clog_volum
type(bief_obj), pointer t2
type(str_line_type) lines
double precision clog_vdiam
double precision clog_tdist
double precision clog_tdiam
double precision, dimension(:), allocatable clog_vlgth
character(len=72) titicecas
type(bief_obj), target iceloc
subroutine section_khione(MESH, IKLE, NELMAX, IFABOR, F, U, V, QVC, CA, CV, H, S, LT)
Definition: section_khione.f:6
integer, dimension(:), allocatable numclog
double precision, dimension(:), allocatable clog_tlgth
subroutine, public clogged_on_bar(RFR0, RFR1, DB, BAR, NBAR, ANG1, FM1, FMT)
subroutine borice(H, U, V, F, AT, LT, DT, TRA05, TRA06, LIUBOR, NPTFR, NUMLIQ, KLOG, MASK, MESH, S)
Definition: borice.f:8
type(bief_obj), target icetype
integer function ielbor(IELM, I)
Definition: ielbor.f:7
type(bief_file), dimension(maxlu_ice) ice_files
double precision, dimension(:), allocatable clog_twdth
type(bief_obj), pointer t1
double precision function p_dsum(MYPART)
Definition: p_dsum.F:7
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
double precision, dimension(:), allocatable clog_vwdth
double precision clog_theta
double precision function bief_sum(X)
Definition: bief_sum.f:7
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
double precision rho_ice
double precision clog_ef
type(bief_obj), pointer t3
Definition: bief.f:3