The TELEMAC-MASCARET system  trunk
freezeup_khione.f
Go to the documentation of this file.
1 ! **********************
2  MODULE freezeup_khione
3 ! **********************
4 !
5 !***********************************************************************
6 ! KHIONE V7P3
7 !***********************************************************************
8 !
9 !brief Module containing all subroutines to deal with the physics
10 !+ of freeze-up and frazil ice growth, at a node level.
11 !+
12 !
13 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
14 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
15 !
17  IMPLICIT NONE
18 !
19  PRIVATE
24 !
25 !=======================================================================
26 !
27  CONTAINS
28 !
29 !=======================================================================
30 !
31 ! *************************
32  SUBROUTINE thermal_growth
33 ! *************************
34 !
35  & ( i,tn,tfrz,dt,srcgm_exp,srcgm_imp,sum_srcgm,sum_frzl,
36  & eps,alpha,nut,constss,limflag)
37 !
38 !***********************************************************************
39 ! KHIONE V8P2
40 !***********************************************************************
41 !
42 !brief Source term for frazil cristals thermal growth/decay
43 !
44 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45 !| I |-->| CELL INDEX
46 !| TN |-->| TRACERS
47 !| TFRZ |-->| FREEZING TEMPERATURE
48 !| DT |-->| TIME STEP
49 !| SRCGM_EXP |<->| EXPLICIT SOURCE OF FRAZIL
50 !| SRCGM_IMP |<->| IMPLICIT SOURCE OF FRAZIL
51 !| SUM_SRCGM |<->| SUM OF THERMAL GROWTH SOURCES (FOR TEMPERATURE SRC)
52 !| SUM_FRZL |<->| SUM OF FRAZIL (TOTAL FRAIL CONC. FOR TEMPERATURE SRC)
53 !| EPS |-->| TURBULENT DISSIPATION RATE
54 !| ALPHA |-->| TURBULENT INTENSITY
55 !| NUT |-->| TURBULENT VISCOSITY
56 !| CONSTSS |-->| 1.D0/(RO0*CP_EAU)
57 !| LIMFLAG |<->| LIMITOR FLAG
58 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59 !
60  USE bief
62  & rho_ice,lh_ice,tc_wt,rk_frzl,ek_frzl,vk_frzl,
63  & nusst,nussi,nuss,itgm,minnk,iseed
64  IMPLICIT NONE
65 !
66 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
67 !
68  TYPE(bief_obj), INTENT(IN) :: TN
69  INTEGER, INTENT(IN) :: I
70  LOGICAL, INTENT(INOUT) :: LIMFLAG
71  DOUBLE PRECISION, INTENT(IN) :: DT,TFRZ,CONSTSS
72  DOUBLE PRECISION, INTENT(IN) :: EPS,ALPHA,NUT
73  DOUBLE PRECISION, INTENT(INOUT) :: SUM_SRCGM,SUM_FRZL
74  DOUBLE PRECISION, INTENT(INOUT) :: SRCGM_EXP(nc_fra)
75  DOUBLE PRECISION, INTENT(INOUT) :: SRCGM_IMP(nc_fra)
76  DOUBLE PRECISION, PARAMETER :: EPS0 = 1.d-8
77 !
78 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
79 !
80  INTEGER :: K,J
81  DOUBLE PRECISION :: TW,QK
82  DOUBLE PRECISION :: CF(nc_fra)
83  DOUBLE PRECISION :: SRCGM(nc_fra)
84  DOUBLE PRECISION :: RM1,RP1,FACT
85  DOUBLE PRECISION :: MAXDT0,SRGMSTB
86 !
87 !-----------------------------------------------------------------------
88 !
89 ! TEMPERATURE AND FRAZIL CONC.
90  tw = tn%ADR(ind_t)%P%R(i)
91  DO k=1,nc_fra
92  j = ind_fra+k-1
93  cf(k) = tn%ADR(j)%P%R(i)
94  ENDDO
95 !
96 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
97 ! COMPUTING THERMAL GROWTH/DECAY
98 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99 !
100 ! ~~> LOOP ON ALL CLASS OF FRAZIL CRISTALS
101  DO k=1,nc_fra
102  j = ind_fra+k-1
103 !
104 ! ~~> COMPUTE NUSSLET NUMBER OF FRAZIL CLASS K
105  IF (nussi.EQ.1) THEN
106  nusst(k) = nuss
107  ELSE IF (nussi.EQ.2) THEN
108  nusst(k) = nusselt(rk_frzl(k), eps, alpha, nut)
109  ENDIF
110 !
111 ! ~~> GROWTH
112  IF(tw.LT.tfrz)THEN
113  qk = nusst(k)*tc_wt/rk_frzl(k)
114  srcgm(k) = 2.d0*qk/(rk_frzl(k)*rho_ice*lh_ice)
115 ! ~~> MELTING
116  ELSE IF(tw.GT.tfrz)THEN
117  qk = nusst(k)*tc_wt/rk_frzl(k)
118  fact = (1.d0/ek_frzl(k)) + (1.d0/rk_frzl(k))
119  srcgm(k) = 2.d0*qk*fact/(rho_ice*lh_ice)
120  IF ((iseed.EQ.1).AND.
121  & (tn%ADR(j)%P%R(i).LE.minnk*vk_frzl(k)/nc_fra)) THEN
122  srcgm(k) = 0.d0
123  ENDIF
124 ! ~~> EQUILIBIRUM
125  ELSE
126  srcgm(k) = 0.d0
127  ENDIF
128 !
129  IF(itgm.EQ.1) THEN
130  srcgm(k) = srcgm(k)*(tfrz-tw)*cf(k)
131  ENDIF
132 !
133  ENDDO
134 !
135  srgmstb = srcgm(1)
136 !
137 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
138 ! COMPUTING MASS EXCHANGES BETWEEN CLASSES
139 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140 !
141  limflag = .false.
142  sum_srcgm = 0.d0
143  sum_frzl = 0.d0
144 !
145 ! ~~> LOOP ON ALL CLASSES OF FRAZIL CRISTALS
146  DO k=1,nc_fra
147  j = ind_fra+k-1
148 !
149 ! ~~~~~~~~~~~~~~~~~~~~~~~
150 ! EXPLICIT IMPLEMENTATION
151 ! ~~~~~~~~~~~~~~~~~~~~~~~
152  IF(itgm.EQ.1) THEN
153 !
154 ! ~~> COMPUTING SOURCE TERM OF FRAZIL CONC. (MONO-CLASS)
155  IF(nc_fra.EQ.1) THEN
156  srcgm_exp(k) = max(-tn%ADR(j)%P%R(i)/dt, srcgm(k))
157 !
158 ! ~~> COMPUTING SOURCE TERM OF FRAZIL CONC. OF CLASS K
159  ELSE
160  rm1 = (rk_frzl(max(k-1, 1))/rk_frzl(k))**3
161  rp1 = (rk_frzl(min(k+1,nc_fra))/rk_frzl(k))**3
162 !
163 ! GROWTH
164  IF(tw.LT.tfrz)THEN
165  IF(k.EQ.1) THEN
166  srcgm_exp(k) =-srcgm(k)/(rp1-1.d0)
167  ELSEIF(k.EQ.nc_fra) THEN
168  srcgm_exp(k) = srcgm(k-1)/(1.d0-rm1)
169  ELSE
170  srcgm_exp(k) = srcgm(k-1)/(1.d0-rm1)-srcgm(k)/(rp1-1.d0)
171  ENDIF
172 !
173 ! MELTING
174  ELSEIF(tw.GT.tfrz)THEN
175  IF(k.EQ.1) THEN
176  srcgm_exp(k) = srcgm(k)-srcgm(k+1)/(rp1-1.d0)
177  ELSEIF(k.EQ.nc_fra) THEN
178  srcgm_exp(k) = srcgm(k)/(1.d0-rm1)
179  ELSE
180  srcgm_exp(k) = srcgm(k)/(1.d0-rm1)-srcgm(k+1)/(rp1-1.d0)
181  ENDIF
182 !
183 ! EQUILIBRIUM
184  ELSE
185  srcgm_exp(k) = 0.d0
186  ENDIF
187 !
188 ! LIMITOR
189  IF (srcgm_exp(k).LE.-tn%ADR(j)%P%R(i)/dt)THEN
190  limflag = .true.
191  srcgm_exp(k) = -tn%ADR(j)%P%R(i)/dt
192  ENDIF
193  ENDIF
194 !
195 ! RATE OF TOTAL VOLUME FRACTION CHANGE
196  sum_srcgm = sum_srcgm + srcgm_exp(k)
197 !
198 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
199 ! SEMI-IMPLICIT IMPLEMENTATION
200 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201  ELSEIF(itgm.EQ.2) THEN
202 !
203 ! ~~> COMPUTING SOURCE TERM OF FRAZIL CONC. (MONO-CLASS)
204  IF(nc_fra.EQ.1) THEN
205  srcgm_exp(k) = 0.d0
206  srcgm_imp(k) = srcgm(k)*(tfrz-tw)
207  srcgm(k) = srcgm(k)*cf(k)
208 !
209 ! ~~> COMPUTING SOURCE TERM OF FRAZIL CONC. OF CLASS K
210  ELSE
211  rm1 = (rk_frzl(max(k-1, 1))/rk_frzl(k))**3
212  rp1 = (rk_frzl(min(k+1,nc_fra))/rk_frzl(k))**3
213 !
214 ! GROWTH
215  IF(tw.LT.tfrz)THEN
216  IF(k.EQ.1) THEN
217  srcgm_exp(k) = 0.d0
218  srcgm_imp(k) =-srcgm(k)*(tfrz-tw)/(rp1-1.d0)
219  ELSEIF(k.EQ.nc_fra) THEN
220  srcgm_exp(k) = srcgm(k-1)*(tfrz-tw)*cf(k-1)/(1.d0-rm1)
221  srcgm_imp(k) = 0.d0
222  ELSE
223  srcgm_exp(k) = srcgm(k-1)*(tfrz-tw)*cf(k-1)/(1.d0-rm1)
224  srcgm_imp(k) =-srcgm(k)*(tfrz-tw)/(rp1-1.d0)
225  ENDIF
226 !
227 ! MELTING
228  ELSEIF(tw.GT.tfrz)THEN
229  IF(k.EQ.1) THEN
230  srcgm_exp(k) =-srcgm(k+1)*(tfrz-tw)*cf(k+1)/(rp1-1.d0)
231  srcgm_imp(k) = srcgm(k)*(tfrz-tw)
232  ELSEIF(k.EQ.nc_fra) THEN
233  srcgm_exp(k) = 0.d0
234  srcgm_imp(k) = srcgm(k)*(tfrz-tw)/(1.d0-rm1)
235  ELSE
236  srcgm_exp(k) =-srcgm(k+1)*(tfrz-tw)*cf(k+1)/(rp1-1.d0)
237  srcgm_imp(k) = srcgm(k)*(tfrz-tw)/(1.d0-rm1)
238  ENDIF
239 !
240 ! EQUILIBRIUM
241  ELSE
242  srcgm_exp(k) = 0.d0
243  srcgm_imp(k) = 0.d0
244  ENDIF
245 !
246 ! LIMITOR
247  IF (srcgm_exp(k).LE.-tn%ADR(j)%P%R(i)/dt)THEN
248  limflag = .true.
249  srcgm_exp(k) = -tn%ADR(j)%P%R(i)/dt
250  ENDIF
251  ENDIF
252 !
253 ! RATE OF TOTAL VOLUME FRACTION CHANGE
254  IF (tfrz.NE.tw) THEN
255  sum_srcgm = sum_srcgm + srcgm_exp(k)/(tfrz-tw)
256  & + srcgm_imp(k)*cf(k)/(tfrz-tw)
257  ENDIF
258 !
259  ENDIF
260 !
261 ! TOTAL FRAZIL CONC.
262  sum_frzl = sum_frzl + cf(k)
263 !
264  ENDDO
265 !
266 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267 ! STABILITY/POSITIVITY WARNING
268 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
269  IF(sum_srcgm.GT.eps0)THEN
270  IF(itgm.EQ.1) THEN
271  maxdt0 = min(2.d0/(rho_ice*lh_ice*constss*srgmstb
272  & *max(eps0, sum_frzl)),
273  & 1.d0/(srgmstb*max(eps0, abs(tw))))
274  ELSEIF(itgm.EQ.2) THEN
275  maxdt0 = 1.d0/(srgmstb*max(eps0, abs(tw)))
276  ENDIF
277 !
278  IF( dt.GT.maxdt0) THEN
279  WRITE(lu,*) ''
280  WRITE(lu,*) 'STABILITY OF THERMAL GROWTH '
281  WRITE(lu,*) 'OR POSITIVITY OF FRAZIL CONC. '
282  WRITE(lu,*) 'MAY BE COMPROMISED '
283  WRITE(lu,*) ''
284  WRITE(lu,*) 'MAKE SURE TIME STEP '
285  WRITE(lu,*) 'IS LOWER THAN ', maxdt0
286  WRITE(lu,*) 'OR CHECK MODEL PARAMETERS '
287  WRITE(lu,*) ''
288  ENDIF
289  ENDIF
290 !
291  RETURN
292 !
293 !-----------------------------------------------------------------------
294 !
295  END SUBROUTINE thermal_growth
296 !
297 !=======================================================================
298 !
299 ! ******************
300  SUBROUTINE seeding
301 ! ******************
302 !
303  & ( i, tn, tfrz, srcse, minflag )
304 !
305 !***********************************************************************
306 ! KHIONE V8P2
307 !***********************************************************************
308 !
309 !brief Computes frazil ice seeding
310 !
311 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
312 !| I |-->| CELL INDEX
313 !| TN |<->| TRACERS
314 !| TFRZ |-->| FREEZING TEMPERATURE
315 !| SRCSE |<->| EXPLICIT SOURCE OF FRAZIL
316 !| MINFLAG |<->| LOGICAL FLAG TO DISABLE SRCSN WHEN MIN CONC. REACHED
317 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318 !
319  USE bief
321  & iseed,seedr,minnk
322  IMPLICIT NONE
323 !
324 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
325 !
326  TYPE(bief_obj), INTENT(INOUT) :: TN
327  LOGICAL, INTENT(INOUT) :: MINFLAG
328  DOUBLE PRECISION, INTENT(INOUT) :: SRCSE(nc_fra)
329  DOUBLE PRECISION, INTENT(IN) :: TFRZ
330  INTEGER, INTENT(IN) :: I
331 !
332 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
333 !
334  INTEGER :: K,J,FROZNB
335 !
336 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
337 !
338 ! ~~~~~~~~~~
339 ! NO SEEDING
340 ! ~~~~~~~~~~
341  IF(iseed.EQ.0) THEN
342  DO k=1,nc_fra
343  srcse(k) = 0.d0
344  ENDDO
345 
346 ! ~~~~~~~~~~~~~~~~~~~~~~~
347 ! MINIMUM CONC. THRESHOLD
348 ! ~~~~~~~~~~~~~~~~~~~~~~~
349  ELSE IF((iseed.EQ.1).OR.(iseed.EQ.3)) THEN
350  froznb = 0
351  DO k=1,nc_fra
352  j = ind_fra+k-1
353  srcse(k) = 0.d0
354 !
355 ! RMK: TRACER VALUES NEED TO BE MODIFIED DIRECTLY OTHERWISE
356 ! SEMI-IMPLICIT THERMAL GROWTH WOULD NOT WORK AS INTENDED
357  tn%ADR(j)%P%R(i) = max(minnk*vk_frzl(k)/nc_fra,
358  & tn%ADR(j)%P%R(i))
359 !
360 ! IN CASE OF MELTING SECONDARY NUCLEATION AND FLOCCULATION
361 ! ARE NOT COMPUTED IF ALL FRAZIL CLASS CONC. ARE LE MIN CONC.
362 ! (ONLY WHEN MIN CONC THRESHOLD IS SELECTED AS SEEDING METHOD)
363  IF ((tn%ADR(ind_t)%P%R(i).GT.tfrz).AND.
364  & (tn%ADR(j)%P%R(i).LE.minnk*vk_frzl(k)/nc_fra)) THEN
365  froznb = froznb + 1
366  ENDIF
367  ENDDO
368  IF (froznb.EQ.nc_fra) THEN
369  minflag = .true.
370  ELSE
371  minflag = .false.
372  ENDIF
373 !
374 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
375 ! CONSTANT SEEDING RATE PER CLASS
376 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
377  ELSE IF((iseed.EQ.2).OR.(iseed.EQ.3)) THEN
378  DO k=1,nc_fra
379  srcse(k) = vk_frzl(k)*seedr
380  ENDDO
381 !
382 ! ~~~~~~~~~~~~~~~~~~~~~~~~
383 ! SEEDING RATE COMPUTED
384 ! FROM ATMOSPHERIC DRIVERS
385 ! ~~~~~~~~~~~~~~~~~~~~~~~~
386 ! TODO SRCSE(K) = SEED%R(I)
387 !
388  ENDIF
389 !
390 !-----------------------------------------------------------------------
391 !
392  RETURN
393 !
394 !-----------------------------------------------------------------------
395 !
396  END SUBROUTINE seeding
397 !
398 !=======================================================================
399 !
400 ! *******************************
401  SUBROUTINE secondary_nucleation
402 ! *******************************
403 !
404  & ( i, tn, dt, srcsn, eps, limflag, minflag )
405 !
406 !***********************************************************************
407 ! KHIONE V8P2
408 !***********************************************************************
409 !
410 !brief Computes secondary nucleation
411 !
412 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
413 !| I |-->| CELL INDEX
414 !| TN |-->| TRACERS
415 !| DT |-->| TIME STEP
416 !| SRCSN |<->| EXPLICIT SOURCE OF FRAZIL
417 !| EPS |-->| TURBULENT DISSIPATION RATE
418 !| LIMFLAG |-->| LOGICAL FLAG TO DISABLE SRCSN WHEN MAX SRCGM REACHED
419 !| MINFLAG |-->| LOGICAL FLAG TO DISABLE SRCSN WHEN MIN CONC. REACHED
420 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
421 !
422  USE bief
423  USE declarations_khione, ONLY : nc_fra, ind_fra, isnuc, xnu,
424  & rk_frzl, vk_frzl, vbk, snnmax
425  IMPLICIT NONE
426 !
427 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
428 !
429  TYPE(bief_obj), INTENT(IN) :: TN
430  INTEGER, INTENT(IN) :: I
431  LOGICAL, INTENT(IN) :: LIMFLAG, MINFLAG
432  DOUBLE PRECISION, INTENT(INOUT) :: SRCSN(nc_fra)
433  DOUBLE PRECISION, INTENT(IN) :: DT,EPS
434 !
435 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
436 !
437  INTEGER :: K,J
438  DOUBLE PRECISION :: SUMSN,WTSQ,CF,WRK,NTILDE
439  DOUBLE PRECISION, PARAMETER :: PI = 4.d0*atan(1.d0)
440 !
441 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
442 !
443 ! ~~~~~~~~~~~~~~~~~~~~~~~
444 ! NO SECONDARY NUCLEATION
445 ! ~~~~~~~~~~~~~~~~~~~~~~~
446  IF((isnuc.EQ.0).OR.(nc_fra.EQ.1).OR.limflag.OR.minflag) THEN
447  DO k=1,nc_fra
448  srcsn(k) = 0.d0
449  ENDDO
450 !
451 ! ~~~~~~~~~~~~~~~~~~~~~~~
452 ! SVENSSON & OMSTEDT 1994
453 ! ~~~~~~~~~~~~~~~~~~~~~~~
454  ELSE IF(isnuc.EQ.1) THEN
455 !
456 ! COMPUTATION OF AVERAGE NUMBER OF ICE CRYSTALS
457  ntilde = 0.d0
458  DO k=1,nc_fra
459  j = ind_fra+k-1
460  cf = tn%ADR(j)%P%R(i)
461  ntilde = ntilde + cf/vk_frzl(k)
462  ENDDO
463  ntilde = min(snnmax, ntilde)
464 !
465  IF(ntilde.GT.1.d0) THEN
466 ! LOOP ON ALL CLASSES EXCEPT 1
467  sumsn = 0.d0
468  DO k=2,nc_fra
469  j = ind_fra+k-1
470  cf = tn%ADR(j)%P%R(i)
471 ! TURBULENT COMPONENT OF RELATIVE VELOCITY
472  wtsq = (4.d0/15.d0)*(eps/xnu)*rk_frzl(k)**2
473 ! RELATIVE VELOCITY
474  wrk = sqrt(wtsq + vbk(k)**2)
475 ! NUCLEATION SINK TERM
476  srcsn(k) = -pi*ntilde*wrk*cf*rk_frzl(k)**2
477 ! LIMITOR
478  srcsn(k) = max(-cf/dt, srcsn(k))
479 ! SUM FOR FIRST CLASS SOURCE
480  sumsn = sumsn + srcsn(k)
481  ENDDO
482 ! NUCLEATION SOURCE TERM FOR FIRST CLASS
483  srcsn(1) = -sumsn
484  ELSE
485  DO k=1,nc_fra
486  srcsn(k) = 0.d0
487  ENDDO
488  ENDIF
489 !
490 ! ~~~~~~~~~~~~~~~~~~~~~~~~~
491 ! WANG & DOERING 2005
492 ! WARNING: NOT CONSERVATIVE
493 ! ~~~~~~~~~~~~~~~~~~~~~~~~~
494  ELSE IF(isnuc.EQ.2) THEN
495 !
496 ! COMPUTATION OF AVERAGE NUMBER OF ICE CRYSTALS
497  ntilde = 0.d0
498  DO k=1,nc_fra
499  j = ind_fra+k-1
500  cf = tn%ADR(j)%P%R(i)
501  ntilde = ntilde + cf/vk_frzl(k)
502  ENDDO
503  ntilde = min(snnmax, ntilde)
504 !
505  IF(ntilde.GT.1.d0) THEN
506 ! LOOP ON ALL CLASSES EXCEPT 1
507  sumsn = 0.d0
508  DO k=2,nc_fra
509  j = ind_fra+k-1
510  cf = tn%ADR(j)%P%R(i)
511 ! TURBULENT COMPONENT OF RELATIVE VELOCITY
512  wtsq = (4.d0/15.d0)*(eps/xnu)*rk_frzl(k)**2
513 ! RELATIVE VELOCITY
514  wrk = sqrt(wtsq + vbk(k)**2)
515 ! NUCLEATION SINK TERM
516  srcsn(k) = -pi*ntilde*wrk*cf*rk_frzl(k)**2
517 ! LIMITOR
518  srcsn(k) = max(-cf/dt, srcsn(k))
519 ! SUM FOR FIRST CLASS SOURCE
520  sumsn = sumsn + srcsn(k)
521  srcsn(k) = srcsn(k)*(rk_frzl(1)/rk_frzl(k))**3
522  ENDDO
523 ! NUCLEATION SOURCE TERM FOR FIRST CLASS
524  srcsn(1) = -sumsn
525  ELSE
526  DO k=1,nc_fra
527  srcsn(k) = 0.d0
528  ENDDO
529  ENDIF
530 !
531  ENDIF
532 !
533 !-----------------------------------------------------------------------
534 !
535  RETURN
536 !
537 !-----------------------------------------------------------------------
538 !
539  END SUBROUTINE secondary_nucleation
540 !
541 !=======================================================================
542 !
543 ! *******************************
544  SUBROUTINE flocculation_breakup
545 ! *******************************
546 !
547  & ( i, tn, dt, srcfb, limflag, minflag )
548 !
549 !***********************************************************************
550 ! KHIONE V8P2
551 !***********************************************************************
552 !
553 !brief Computes floculation & breaking
554 !
555 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
556 !| I |-->| CELL INDEX
557 !| TN |-->| TRACERS
558 !| DT |-->| TIME STEP
559 !| SRCFB |<->| EXPLICIT SOURCE OF FRAZIL
560 !| LIMFLAG |-->| LOGICAL FLAG TO DISABLE SRCFB WHEN MAX SRCGM REACHED
561 !| MINFLAG |-->| LOGICAL FLAG TO DISABLE SRCFB WHEN MIN CONC. REACHED
562 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
563 !
564  USE bief
566  & ifloc, afloc
567  IMPLICIT NONE
568 !
569 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
570 !
571  TYPE(bief_obj), INTENT(IN) :: TN
572  LOGICAL, INTENT(IN) :: LIMFLAG, MINFLAG
573  INTEGER, INTENT(IN) :: I
574  DOUBLE PRECISION, INTENT(IN) :: DT
575  DOUBLE PRECISION, INTENT(INOUT) :: SRCFB(nc_fra)
576 !
577 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
578 !
579  INTEGER :: K,J
580  DOUBLE PRECISION :: CF, CF0
581  DOUBLE PRECISION :: BETA, BETA0
582 !
583 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
584 !
585 ! ~~~~~~~~~~~~~~~~~~~~~~~
586 ! NO FLOCULATION & BREAUP
587 ! ~~~~~~~~~~~~~~~~~~~~~~~
588  IF((ifloc.EQ.0).OR.(nc_fra.EQ.1).OR.limflag.OR.minflag) THEN
589  DO k=1,nc_fra
590  srcfb(k) = 0.d0
591  ENDDO
592 !
593 ! ~~~~~~~~~~~~~~~~~~~~~~~
594 ! SVENSSON & OMSTEDT 1994
595 ! ~~~~~~~~~~~~~~~~~~~~~~~
596  ELSE IF(ifloc.EQ.1) THEN
597  DO k=1,nc_fra
598  j = ind_fra+k-1
599 !
600  IF(k.EQ.1) THEN
601  cf = tn%ADR(j)%P%R(i)
602  beta = afloc*rk_frzl(k)/rk_frzl(1)
603  beta = min(1.d0/dt, beta) !LIMITOR
604  srcfb(k) = -beta*cf
605 !
606  ELSE IF(k.EQ.nc_fra) THEN
607  beta0 = beta
608  cf0 = tn%ADR(j-1)%P%R(i)
609  srcfb(k) = beta0*cf0
610 !
611  ELSE
612  cf = tn%ADR(j)%P%R(i)
613  cf0 = tn%ADR(j-1)%P%R(i)
614  beta0 = beta
615  beta = afloc*rk_frzl(k)/rk_frzl(1)
616  beta = min(1.d0/dt, beta) !LIMITOR
617  srcfb(k) = -beta*cf + beta0*cf0
618  ENDIF
619  ENDDO
620 !
621  ENDIF
622 !
623 !-----------------------------------------------------------------------
624 !
625  RETURN
626 !
627 !-----------------------------------------------------------------------
628 !
629  END SUBROUTINE flocculation_breakup
630 !
631 !=======================================================================
632 !
633 ! *******************************
634  SUBROUTINE turbulent_parameters
635 ! *******************************
636 !
637  & ( i, vmag, h, kt, eps, alpha, nut, cf, ak, ep, iturb_tel)
638 !
639 !***********************************************************************
640 ! KHIONE V8P2
641 !***********************************************************************
642 !
643 !brief Computes turbulent parameters for frazil ice model
644 !
645 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
646 !| AK |-->| TURBULENT KINETIC ENERGY K AT TIME T(N+1)
647 !| EP |-->| TURBULENT ENERGY DISSIPATION AT TIME T(N+1)
648 !| ITURB_TEL|-->| TURBULENT MODEL IN T2D
649 !| VMAG |-->| VELOCITY MAGNITUDE
650 !| H |-->| WATER DEPTH
651 !| KT |<->| MEAN TURBULENT KINETIC ENERGY
652 !| EPS |<->| MEAN TURBULENT DISSIPATION RATE
653 !| ALPHA |<->| TURBULENT INTENSITY
654 !| NUT |<->| TURBULENT VISCOSITY
655 !| CF |<->| QUADRATIC FRICTION COEFFICIENT
656 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
657 !
658  USE bief
659  USE declarations_khione, ONLY : xnu, iturb
660  IMPLICIT NONE
661 !
662 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
663 !
664  TYPE(bief_obj), INTENT(IN) :: AK,EP
665  INTEGER, INTENT(IN) :: I,ITURB_TEL
666  DOUBLE PRECISION, INTENT(IN) :: VMAG,H,CF
667  DOUBLE PRECISION, INTENT(INOUT) :: KT,EPS,ALPHA,NUT
668 !
669 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
670 !
671  DOUBLE PRECISION, PARAMETER :: KARMAN_CONST = 0.4d0
672  DOUBLE PRECISION, PARAMETER :: HEPS = 1.d-3
673  DOUBLE PRECISION, PARAMETER :: UEPS = 1.d-6
674  DOUBLE PRECISION :: USTAR,FACT0,FACT1,HTEMP
675 !
676 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
677 !
678 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
679 ! CONSTANTS (FOR TEST PURPOSES)
680 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
681  IF(iturb.EQ.0) THEN
682  kt = 9.6d-4
683  eps = 12.d-4
684 !
685 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
686 ! DEPTH AVERAGED K, EPS
687 ! FROM MIXING LENGTH VERTICAL PROFILE
688 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689  ELSE IF(iturb.EQ.1) THEN
690  ustar = 0.5d0 * cf * vmag * vmag
691  ustar = sqrt(ustar)
692  htemp = max(h, heps)
693  fact0 = (htemp/2.d0 - xnu/ustar + xnu**2/(2.d0*htemp*ustar**2))
694  fact1 = (log(htemp*ustar/xnu) - 1.d0 + xnu/(htemp*ustar))
695  kt = (ustar**2)*fact0/0.3d0
696  eps = (ustar**3)*fact1/karman_const
697 !
698 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
699 ! DEPTH AVERAGED K, EPS
700 ! FROM TELEMAC2D K-EPS MODEL
701 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
702  ELSE IF(iturb.EQ.2) THEN
703  IF(iturb_tel.NE.3) THEN
704  WRITE(lu,*) ''
705  WRITE(lu,*) 'TURBULENCE ESTIMATION MODEL', iturb
706  WRITE(lu,*) 'MUST BE USED WITH K-EPS MODEL'
707  WRITE(lu,*) ''
708  WRITE(lu,*) 'PLEASE SET TURBULENCE MODEL = 3'
709  WRITE(lu,*) 'IN T2D STEERING FILE'
710  WRITE(lu,*) ''
711  CALL plante(1)
712  stop
713  ENDIF
714  kt = ak%R(i)
715  eps = ep%R(i)
716 !
717 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
718  ELSE
719  WRITE(lu,*) ''
720  WRITE(lu,*) 'TURBULENCE ESTIMATION MODEL', iturb
721  WRITE(lu,*) 'NOT IMPLEMENTED'
722  WRITE(lu,*) ''
723  CALL plante(1)
724  stop
725  ENDIF
726 !
727  nut = 0.09d0*kt**2/eps
728  alpha = sqrt(2.d0*kt)/max(ueps, vmag)
729 !
730 !-----------------------------------------------------------------------
731 !
732  RETURN
733 !
734 !-----------------------------------------------------------------------
735 !
736  END SUBROUTINE turbulent_parameters
737 !
738 !=======================================================================
739 !
740 ! *********************************
741  DOUBLE PRECISION FUNCTION nusselt
742 ! *********************************
743 !
744  & ( radius, eps, alpha, nut )
745 !
746 !***********************************************************************
747 ! KHIONE V8P2
748 !***********************************************************************
749 !
750 !brief Computes the Nusselt number for thermal growth
751 !
752 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
753 !| RADIUS |-->| RADIUS OF FRAZIL PARTICLE
754 !| EPS |-->| TURBULENT DISSIPATION RATE
755 !| ALPHA |-->| TURBULENT INTENSITY
756 !| NUT |-->| TURBULENT VISCOSITY
757 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
758 !
759  USE declarations_khione, ONLY : tc_wt, cp_eau, ro0, xnu
760  IMPLICIT NONE
761 !
762 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
763 !
764  DOUBLE PRECISION, INTENT(IN) :: RADIUS, ALPHA, EPS, NUT
765  DOUBLE PRECISION :: MSTAR, KOLML, PR, PR12, PR13
766 !
767 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
768 !
769  kolml = (xnu**3/eps)**(0.25d0)
770  mstar = radius/kolml
771  pr = ro0*cp_eau*nut/tc_wt
772  pr12 = sqrt(pr)
773  pr13 = pr**(1.d0/3.d0)
774 !
775  IF (mstar .LE. 1.d0/pr12) THEN
776  nusselt = 1.d0 + 0.17d0*mstar*pr12
777  ELSE IF ((1.d0/pr12 .LT. mstar) .AND. (mstar .LE. 1.d0)) THEN
778  nusselt = 1.d0 + 0.55d0*(mstar**(2.d0/3.d0))*pr13
779  ELSE IF (mstar .GT. 1.d0) THEN
780  IF (alpha*mstar**(4.d0/3.d0) .LE. 1.d3) THEN
781  nusselt = 1.1d0 + 0.77d0*(alpha**0.035d0)
782  & *(mstar**(2.d0/3.d0))*pr13
783  ELSE
784  nusselt = 1.1d0 + 0.77d0*(alpha**0.25d0)*mstar*pr13
785  ENDIF
786  ELSE
787  nusselt = 1.d0
788  ENDIF
789 !
790  RETURN
791 !
792 !-----------------------------------------------------------------------
793 !
794  END FUNCTION nusselt
795 !
796 !=======================================================================
797 !
798 ! *******************************************
799  DOUBLE PRECISION FUNCTION buoyancy_velocity
800 ! *******************************************
801 !
802  & ( radius, thickness )
803 !
804 !***********************************************************************
805 ! KHIONE V8P2
806 !***********************************************************************
807 !
808 !brief Computes the buoyancy velocity of frazil particles
809 !
810 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811 !| RADIUS |-->| RADIUS OF FRAZIL PARTICLE
812 !| THICKNESS |-->| THICKNESS OF FRAZIL PARTICLE
813 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
814 !
815  USE declarations_khione, ONLY : ibuoy, xnu, de,
816  & ro0, rho_ice
817  IMPLICIT NONE
818 !
819 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
820 !
821  DOUBLE PRECISION, INTENT(IN) :: RADIUS, THICKNESS
822  DOUBLE PRECISION, PARAMETER :: CKD = 2.d0
823  DOUBLE PRECISION, PARAMETER :: GRAV = 9.81d0
824  DOUBLE PRECISION :: KV,GP,PI
825 !
826 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
827 !
828  IF(ibuoy.EQ.1) THEN
829 ! DALY (1984)
830  pi = 4.d0*atan(1.d0)
831  kv = 2.d0*pi/de
832  gp = 2.d0*kv*grav*(ro0-rho_ice)/(pi*ro0)
833  IF (radius .LE. 3.e-4) THEN
834  buoyancy_velocity = 0.08d0*gp*(radius**2)/xnu
835  ELSEIF ((3.e-4 .LT. radius).AND.(radius .LE. 1.4e-3)) THEN
836  buoyancy_velocity = 0.16d0*(gp**0.715d0)
837  & *(xnu**(-0.428d0))*(radius**1.14d0)
838  ELSEIF (1.4e-3 .LT. radius) THEN
839  buoyancy_velocity = sqrt(gp*radius/2.d0)
840  ENDIF
841 !
842  ELSEIF(ibuoy.EQ.2) THEN
843 ! DALY (1984) INTERMEDIATE RANGE
844 ! RMK: VALID FOR 0.03 < R < 0.14 CM
845  buoyancy_velocity = 1.d-2*30.d0*(radius*1.d2)**1.2d0
846 
847  ELSEIF(ibuoy.EQ.3) THEN
848 ! MATOUSEK (1992)
849  buoyancy_velocity = (1.31d-5)*((2.d0*radius)**0.29d0)
850  & *(thickness**0.61d0)/xnu
851 
852  ELSEIF(ibuoy.EQ.4) THEN
853 ! GOSIK & OSTERKAMP (1983)
854 ! RMK1: STOKES REGIME (ONLY VALID FOR R < 0.03 CM)
855 ! RMK2: DRAG COEFFICIENT CKD NEEDS TO BE ADJUSTED
856  buoyancy_velocity = sqrt(4.d0*(ro0-rho_ice)
857  & *grav*radius/(de*ro0*ckd))
858 !
859  ENDIF
860 !
861  RETURN
862 !
863 !-----------------------------------------------------------------------
864 !
865  END FUNCTION buoyancy_velocity
866 !
867 !=======================================================================
868 !
869 ! ***************************************
870  DOUBLE PRECISION FUNCTION melting_point
871 ! ***************************************
872 !
873  & ( sal )
874 !
875 !***********************************************************************
876 ! KHIONE V7P2
877 !***********************************************************************
878 !
879 !brief Computes the melting point as a function of salinity
880 !
881 !reference:
882 !+ Millero, F.J. “Freezing point of seawater”. In “Eighth Report of
883 !+ the Joint Panel on Oceanographic Tables and Standards”,
884 !+ UNESCO Tech. Pap. Mar. Sci. No.28, Annex 6. UNESCO, Paris, 1978
885 !
886 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
887 !| SAL |-->| SALINITY CONCENTRATION
888 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
889 !
890  IMPLICIT NONE
891 !
892 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
893 !
894  DOUBLE PRECISION, INTENT(IN) :: SAL
895 !
896 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
897 !
898 !
899  melting_point =
900  & -0.0575d0*sal + 0.001710523d0*( sal**1.5d0 )
901  & -0.0002154996d0*( sal**2 )
902 ! & -0.00753 * PRESSURE BELOW SURFACE
903 !
904  RETURN
905 !
906 !-----------------------------------------------------------------------
907 !
908  END FUNCTION melting_point
909 !
910 !=======================================================================
911 !
912 ! *****************************
913  SUBROUTINE erosion_deposition
914 ! *****************************
915 !
916  & ( frzl,srcgm,theta0,theta1,beta1,vbb,thifemf,hun,
917  & anfem,vmag,depth,isbar )
918 !
919 !***********************************************************************
920 ! KHIONE V7P2
921 !***********************************************************************
922 !
923 !brief Computes frazil ice sink/source term due to deposition/erosion
924 ! under the ice cover
925 !
926 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
927 !| ANFEM |-->| CONCENTRATION OF SURFACE ICE PARTICLES
928 !| BETA1 |-->| RATE OF REENTRAINMENT OF SURFACE ICE PER UNIT AREA
929 !| DEPTH |-->| WATER DEPTH
930 !| FRZL |-->| FRAZIL CONCENTRATION
931 !| HUN |-->| UNCOVERED ICE THICKNESS
932 !| ISBAR |-->| WHETHER OR NOT ICE SETTLES AT THE SURFACE OR ON A BAR
933 !| SRCGM |<->| EXPLICIT SOURCE OF FRAZIL
934 !| THETA0 |-->| PROBABILITY OF DEPOSITION OF FRAZIL - OPEN WATER
935 !| THETA1 |-->| PROBABILITY OF DEPOSITION OF FRAZIL - ICE
936 !| THIFEMF |-->| FRAZIL ICE THICKNESS
937 !| VBB |-->| BUOYANCY VELOCITY OF FRAZIL GRANULES
938 !| VMAG |-->| CURRENT VELOCITY MAGNITUDE
939 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
940 !
941  USE declarations_khione, ONLY : af, surf_ef
942  IMPLICIT NONE
943 !
944 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
945 !
946  DOUBLE PRECISION, INTENT(IN) :: FRZL,ANFEM
947  DOUBLE PRECISION, INTENT(INOUT) :: SRCGM
948  DOUBLE PRECISION, INTENT(IN) :: THETA0,THETA1,BETA1
949  DOUBLE PRECISION, INTENT(IN) :: VMAG,DEPTH,VBB,THIFEMF,HUN
950  LOGICAL, INTENT(IN) :: ISBAR
951 !
952 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
953 !
954  DOUBLE PRECISION HE,PV
955 !
956 !-----------------------------------------------------------------------
957 !
958 ! NODES ON THE BAR: NO DEPOSITION ON THE WATER SURFACE
959  IF( isbar ) THEN
960  pv = ( af*vmag )/depth * frzl
961 ! NO BAR: DEPOSITION ON WATER SURFACE
962  ELSE
963  pv = ( (theta0*vbb*(1.d0-anfem) ) +
964  & ( theta1*vbb*anfem) )/depth * frzl
965  ENDIF
966 !
967 ! EROSION / DEPOSITION SOURCE TERM
968  he = ( 1.d0-surf_ef ) * ( thifemf + hun )
969  srcgm = (beta1*he*anfem)/depth - pv
970 !
971  RETURN
972 !
973 !-----------------------------------------------------------------------
974 !
975  END SUBROUTINE erosion_deposition
976 !
977 !=======================================================================
978 !
979 ! ************************
980  SUBROUTINE precipitation
981 ! ************************
982 !
983  & (vbb, frzl, srcp, vmag, g, rk, it)
984 !
985 !***********************************************************************
986 ! KHIONE V8P2
987 !***********************************************************************
988 !
989 !brief Computes frazil precipitation on the surface
990 !
991 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
992 !| FRZL |-->| FRAZIL CONCENTRATION
993 !| G |-->| GRAVITY CONSTANT
994 !| IT |-->| LOGICAL FOR ICE COVER
995 !| RK |-->| FRAZIL RADIUS
996 !| SRCP |<->| EXPLICIT SOURCE OF FRAZIL
997 !| VBB |-->| BUOYANCY VELOCITY OF FRAZIL GRANULES
998 !| VMAG |-->| CURRENT VELOCITY MAGNITUDE
999 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1000 !
1001  USE declarations_khione, ONLY : ro0, rho_ice
1002  IMPLICIT NONE
1003 !
1004 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1005 !
1006  LOGICAL, INTENT(IN) :: IT
1007  DOUBLE PRECISION, INTENT(IN) :: FRZL,G,RK
1008  DOUBLE PRECISION, INTENT(INOUT) :: SRCP
1009  DOUBLE PRECISION, INTENT(IN) :: VMAG,VBB
1010 !
1011 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1012 !
1013  DOUBLE PRECISION UC2,TC,M
1014 !
1015 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1016 !
1017 ! COMPUTATION OF THE CRITICAL VELOCITY (USEFUL IF ICE COVER IS
1018 ! PRESENT)
1019  IF(it) THEN
1020  tc = 2.d2/3.d0
1021  uc2 = tc*(ro0-rho_ice)*g*rk/ro0
1022  m = max(0.d0,1.d0-(vmag*vmag/uc2))
1023  srcp = -vbb*frzl*m
1024  ELSE
1025  srcp = -vbb*frzl
1026  ENDIF
1027 !
1028  RETURN
1029 !
1030 !-----------------------------------------------------------------------
1031 !
1032  END SUBROUTINE precipitation
1033 !
1034 !=======================================================================
1035 !
1036 ! *************************
1037  SUBROUTINE clogged_on_bar
1038 ! *************************
1039 !
1040  & ( rfr0,rfr1,db,bar,nbar,ang1,fm1,fmt )
1042 !***********************************************************************
1043 ! KHIONE V7P3
1044 !***********************************************************************
1045 !
1046 !brief Computes the accumulated mass on vertical and/or
1047 ! transverse bars.
1048 !
1049 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1050 !| ANG1 |-->| ACCUMULATION ANGLE
1051 !| BAR |-->| PHYSICAL SIZE OF LENGTHS OF THE BAR
1052 !| DB |-->| BAR DIAMETER
1053 !| FM1 |<->| MASS OF FRAZIL ICE ACCULATED ON ONE BAR
1054 !| FMT |<->| TOTAL MASS ACCUMULATED
1055 !| NBAR |-->| REPRESENTATIVE NUMBER OF BAR
1056 !| RFR0 |-->| INITIAL RADIUS BEFORE THE FRAZIL ICE ACCUMULATION
1057 !| RFR1 |-->| RADIUS OF THE FRAZIL ICE ACCUMULATION
1058 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1059 !
1060  USE declarations_khione, ONLY : rho_ice,clog_ef
1061  IMPLICIT NONE
1062 !
1063 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1064 !
1065  DOUBLE PRECISION, INTENT(IN) :: NBAR
1066  DOUBLE PRECISION, INTENT(IN) :: RFR0,RFR1,BAR,ANG1,DB
1067  DOUBLE PRECISION, INTENT(INOUT) :: FM1,FMT
1068 !
1069 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1070 !
1071  DOUBLE PRECISION :: RB,DR,AMNU,ASH,RATIO,ADONT,AICE
1072 !
1073 !-----------------------------------------------------------------------
1074 !
1075 ! RB = RFR0/2.0/COS(ANG1)
1076  rb = db/2.d0
1077  dr = 2.d0*rb - rfr0
1078 !
1079  amnu = rfr0**2*cos(ang1)*sin(ang1) + 2.d0*ang1*rb**2 -
1080  & rb**2*cos(2.d0*ang1)*sin(2.d0*ang1)
1081  ash = 2.d0*ang1*rb**2 - (ang1*rfr0**2 - rb*rfr0*sin(ang1)) ! SHADE AREA
1082 !
1083  IF (rfr1.GT.2.d0*rb) THEN
1084  aice = ang1*rfr1**2 - amnu
1085  ELSE
1086  adont = ( rfr1**2 - rfr0**2 )*ang1
1087  ratio = ( rfr1 - rfr0 ) / dr
1088  aice = adont - ratio*ash
1089  ENDIF
1090 !
1091  fm1 = aice*bar*rho_ice*( 1.d0-clog_ef )
1092  fmt = fm1*nbar
1093 !
1094  RETURN
1095 !
1096 !-----------------------------------------------------------------------
1097 !
1098  END SUBROUTINE clogged_on_bar
1099 !
1100 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1101 !
1102  END MODULE freezeup_khione
subroutine, public secondary_nucleation(I, TN, DT, SRCSN, EPS, LIMFLAG, MINFLAG)
subroutine, public erosion_deposition(FRZL, SRCGM, THETA0, THETA1, BETA1, VBB, THIFEMF, HUN, ANFEM, VMAG, DEPTH, ISBAR)
double precision function, public buoyancy_velocity(RADIUS, THICKNESS)
double precision surf_ef
double precision function, public melting_point(SAL)
subroutine, public clogged_on_bar(RFR0, RFR1, DB, BAR, NBAR, ANG1, FM1, FMT)
double precision, dimension(:), allocatable vk_frzl
double precision tc_wt
double precision cp_eau
double precision, dimension(:), allocatable rk_frzl
subroutine, public seeding(I, TN, TFRZ, SRCSE, MINFLAG)
subroutine, public flocculation_breakup(I, TN, DT, SRCFB, LIMFLAG, MINFLAG)
subroutine, public precipitation(VBB, FRZL, SRCP, VMAG, G, RK, IT)
subroutine, public turbulent_parameters(I, VMAG, H, KT, EPS, ALPHA, NUT, CF, AK, EP, ITURB_TEL)
subroutine, public thermal_growth(I, TN, TFRZ, DT, SRCGM_EXP, SRCGM_IMP, SUM_SRCGM, SUM_FRZL, EPS, ALPHA, NUT, CONSTSS, LIMFLAG)
double precision rho_ice
double precision function, public nusselt(RADIUS, EPS, ALPHA, NUT)
double precision clog_ef
Definition: bief.f:3