Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
MODULE icbdyn
!!======================================================================
!! *** MODULE icbdyn ***
!! Iceberg: time stepping routine for iceberg tracking
!!======================================================================
!! History : 3.3 ! 2010-01 (Martin&Adcroft) Original code
!! - ! 2011-03 (Madec) Part conversion to NEMO form
!! - ! Removal of mapping from another grid
!! - ! 2011-04 (Alderson) Split into separate modules
!! - ! 2011-05 (Alderson) Replace broken grounding routine with one of
!! - ! Gurvan's suggestions (just like the broken one)
!!----------------------------------------------------------------------
USE par_oce ! NEMO parameters
USE dom_oce ! NEMO ocean domain
USE phycst ! NEMO physical constants
USE in_out_manager ! IO parameters
!
USE icb_oce ! define iceberg arrays
USE icbutl ! iceberg utility routines
USE icbdia ! iceberg budget routines
IMPLICIT NONE
PRIVATE
PUBLIC icb_dyn ! routine called in icbstp.F90 module
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: icbdyn.F90 15088 2021-07-06 13:03:34Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE icb_dyn( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE icb_dyn ***
!!
!! ** Purpose : iceberg evolution.
!!
!! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt !
!
LOGICAL :: ll_bounced
REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
REAL(wp) :: zuvel_n, zvvel_n, zxi_n , zyj_n
REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2
TYPE(iceberg), POINTER :: berg
TYPE(point) , POINTER :: pt
!!----------------------------------------------------------------------
!
! 4th order Runge-Kutta to solve: d/dt X = V, d/dt V = A
! with I.C.'s: X=X1 and V=V1
!
! ; A1=A(X1,V1)
! X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
! X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
! X4 = X1+ dt*V3 ; V4 = V1+ dt*A3 ; A4=A(X4,V4)
!
! Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
! Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
! time steps
zdt = berg_dt
zdt_2 = zdt * 0.5_wp
zdt_6 = zdt / 6._wp
berg => first_berg ! start from the first berg
!
DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==!
!
pt => berg%current_point
ll_bounced = .FALSE.
! STEP 1 !
! ====== !
zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s
zyj1 = pt%yj ; zvvel1 = pt%vvel
! !** A1 = A(X1,V1)
CALL icb_accel( kt, berg , zxi1, ze1, zuvel1, zuvel1, zax1, &
& zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 0.5_wp )
!
zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt
zv1 = zvvel1 / ze2
! STEP 2 !
! ====== !
! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1
! position using di/dt & djdt ! V2 in m/s
zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1
zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1
!
CALL icb_ground( berg, zxi2, zxi1, zu1, &
& zyj2, zyj1, zv1, ll_bounced )
! !** A2 = A(X2,V2)
CALL icb_accel( kt, berg , zxi2, ze1, zuvel2, zuvel1, zax2, &
& zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 0.5_wp )
!
zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt
zv2 = zvvel2 / ze2
!
! STEP 3 !
! ====== !
! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3)
zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2
zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2
!
CALL icb_ground( berg, zxi3, zxi1, zu2, &
& zyj3, zyj1, zv2, ll_bounced )
! !** A3 = A(X3,V3)
CALL icb_accel( kt, berg , zxi3, ze1, zuvel3, zuvel1, zax3, &
& zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1._wp )
!
zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt
zv3 = zvvel3 / ze2
! STEP 4 !
! ====== !
! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3
zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3
zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3
CALL icb_ground( berg, zxi4, zxi1, zu3, &
& zyj4, zyj1, zv3, ll_bounced )
! !** A4 = A(X4,V4)
CALL icb_accel( kt, berg , zxi4, ze1, zuvel4, zuvel1, zax4, &
& zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1._wp )
zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt
zv4 = zvvel4 / ze2
! FINAL STEP !
! ========== !
! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 )
zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 )
zuvel_n = pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 )
zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 )
CALL icb_ground( berg, zxi_n, zxi1, zuvel_n, &
& zyj_n, zyj1, zvvel_n, ll_bounced )
pt%uvel = zuvel_n !** save in berg structure
pt%vvel = zvvel_n
pt%xi = zxi_n
pt%yj = zyj_n
berg => berg%next ! switch to the next berg
!
END DO !== end loop over all bergs ==!
!
END SUBROUTINE icb_dyn
SUBROUTINE icb_ground( berg, pi, pi0, pu, &
& pj, pj0, pv, ld_bounced )
!!----------------------------------------------------------------------
!! *** ROUTINE icb_ground ***
!!
!! ** Purpose : iceberg grounding.
!!
!! ** Method : - adjust velocity and then put iceberg back to start position
!! NB two possibilities available one of which is hard-coded here
!!----------------------------------------------------------------------
TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg
!
REAL(wp), INTENT(inout) :: pi , pj ! current iceberg position
REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position
REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities
LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator
!
INTEGER :: ii, ii0
INTEGER :: ij, ij0
INTEGER :: ikb
INTEGER :: ibounce_method
!
REAL(wp) :: zD
REAL(wp), DIMENSION(jpk) :: ze3t
!!----------------------------------------------------------------------
!
ld_bounced = .FALSE.
!
ii0 = INT( pi0+0.5 ) + (nn_hls-1) ; ij0 = INT( pj0+0.5 ) + (nn_hls-1) ! initial gridpoint position (T-cell)
ii = INT( pi +0.5 ) + (nn_hls-1) ; ij = INT( pj +0.5 ) + (nn_hls-1) ! current - -
!
IF( ii == ii0 .AND. ij == ij0 ) RETURN ! berg remains in the same cell
!
! map into current processor
ii0 = mi1( ii0 )
ij0 = mj1( ij0 )
ii = mi1( ii )
ij = mj1( ij )
!
! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
IF ( ln_M2016 .AND. ln_icb_grd ) THEN
!
! draught (keel depth)
zD = rho_berg_1_oce * berg%current_point%thickness
!
! interpol needed data
CALL icb_utl_interp( pi, pj, pe3t=ze3t )
!
!compute bottom level
CALL icb_utl_getkb( ikb, ze3t, zD )
!
! berg reach a new t-cell, but an ocean one
! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0)
IF( tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN
!
ELSE
IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one
END IF
!
! From here, berg have reach land: treat grounding/bouncing
! -------------------------------
ld_bounced = .TRUE.
!! not obvious what should happen now
!! if berg tries to enter a land box, the only location we can return it to is the start
!! position (pi0,pj0), since it has to be in a wet box to do any melting;
!! first option is simply to set whole velocity to zero and move back to start point
!! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
!! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast
ibounce_method = 2
SELECT CASE ( ibounce_method )
CASE ( 1 )
pi = pi0
pj = pj0
pu = 0._wp
pv = 0._wp
CASE ( 2 )
IF( ii0 /= ii ) THEN
pi = pi0 ! return back to the initial position
pu = 0._wp ! zeroing of velocity in the direction of the grounding
ENDIF
IF( ij0 /= ij ) THEN
pj = pj0 ! return back to the initial position
pv = 0._wp ! zeroing of velocity in the direction of the grounding
ENDIF
END SELECT
!
END SUBROUTINE icb_ground
SUBROUTINE icb_accel( kt, berg , pxi, pe1, puvel, puvel0, pax, &
& pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale )
!!----------------------------------------------------------------------
!! *** ROUTINE icb_accel ***
!!
!! ** Purpose : compute the iceberg acceleration.
!!
!! ** Method : - sum the terms in the momentum budget
!!----------------------------------------------------------------------
TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg
INTEGER , INTENT(in ) :: kt ! time step
REAL(wp) , INTENT(in ) :: pcfl_scale
REAL(wp) , INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential
REAL(wp) , INTENT(in ) :: puvel , pvvel ! berg velocity [m/s]
REAL(wp) , INTENT(in ) :: puvel0, pvvel0 ! initial berg velocity [m/s]
REAL(wp) , INTENT( out) :: pe1, pe2 ! horizontal scale factor at (xi,yj)
REAL(wp) , INTENT(inout) :: pax, pay ! berg acceleration
REAL(wp) , INTENT(in ) :: pdt ! berg time step
!
REAL(wp), PARAMETER :: pp_alpha = 0._wp !
REAL(wp), PARAMETER :: pp_beta = 1._wp !
REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed
REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration
REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp !
!
INTEGER :: itloop, ikb, jk
REAL(wp) :: zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi
REAL(wp) :: zvo, zssv, zvi, zva, zvwave, zssh_y
REAL(wp) :: zff, zT, zD, zW, zL, zM, zF
REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
REAL(wp) :: z_ocn, z_atm, z_ice, zdep
REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw
!!----------------------------------------------------------------------
! Interpolate gridded fields to berg
nknberg = berg%number(1)
CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor
& pssu=zssu, pui=zui, pua=zua, & ! oce/ice/atm velocities
& pssv=zssv, pvi=zvi, pva=zva, & ! oce/ice/atm velocities
& pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient
& phi=zhi, pff=zff) ! ice thickness and coriolis
zM = berg%current_point%mass
zT = berg%current_point%thickness ! total thickness
zD = rho_berg_1_oce * zT ! draught (keel depth)
zF = zT - zD ! freeboard
zW = berg%current_point%width
zL = berg%current_point%length
zhi = MIN( zhi , zD )
zD_hi = MAX( 0._wp, zD-zhi )
! Wave radiation
zuwave = zua - zssu ; zvwave = zva - zssv ! Use wind speed rel. to ocean for wave model
zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current;
! ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
zampl = 0.5_wp * 0.02025_wp * zwmod ! This is "a", the wave amplitude
zLwavelength = 0.32_wp * zwmod ! Surface wave length fitted to data in table at
! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
zLcutoff = 0.125_wp * zLwavelength
zLtop = 0.25_wp * zLwavelength
zCr = pp_Cr0 * MIN( MAX( 0._wp, (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1._wp) ! Wave radiation coefficient
! ! fitted to graph from Carrieres et al., POAC Drift Model.
zwave_rad = 0.5_wp * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2._wp*zW*zL) / (zW+zL)
zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed
IF( zwmod /= 0._wp ) THEN
zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ?
zvwave = zva/zwmod
ELSE
zuwave = 0._wp ; zvwave=0._wp ; zwave_rad=0._wp ! ... and only when wind is present. !!gm wave_rad=0. is useless
ENDIF
! Weighted drag coefficients
z_ocn = pp_rho_seawater / zM * (0.5_wp*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
z_atm = pp_rho_air / zM * (0.5_wp*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL)
z_ice = pp_rho_ice / zM * (0.5_wp*pp_Cd_iv*zW*zhi )
IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp
! lateral velocities
! default ssu and ssv
! ln_M2016: mean velocity along the profile
IF ( ln_M2016 ) THEN
! interpol needed data
CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) ! 3d velocities
!compute bottom level
CALL icb_utl_getkb( ikb, ze3t, zD )
! compute mean velocity
CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb)
CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb)
ELSE
zuo = zssu
zvo = zssv
END IF
zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel
!
DO itloop = 1, 2 ! Iterate on drag coefficients
!
zus = 0.5_wp * ( zuveln + puvel )
zvs = 0.5_wp * ( zvveln + pvvel )
zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
!
! Explicit accelerations
!zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
! -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
!zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
! -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
zaxe = -grav * zssh_x + zwave_rad * zuwave
zaye = -grav * zssh_y + zwave_rad * zvwave
IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest
zaxe = zaxe + zff*pvvel0
zaye = zaye - zff*puvel0
ELSE
zaxe = zaxe + zff*pvvel
zaye = zaye - zff*puvel
ENDIF
IF( pp_beta > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest
zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
ELSE
zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
ENDIF
! Solve for implicit accelerations
IF( pp_alpha + pp_beta > 0._wp ) THEN
zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
zA11 = 1._wp + pp_beta *pdt*zlambda
zA12 = pp_alpha*pdt*zff
zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 )
pax = zdetA * ( zA11*zaxe + zA12*zaye )
pay = zdetA * ( zA11*zaye - zA12*zaxe )
ELSE
pax = zaxe ; pay = zaye
ENDIF
zuveln = puvel0 + pdt*pax
zvveln = pvvel0 + pdt*pay
!
END DO ! itloop
IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked)
zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg
IF( zspeed > 0._wp ) THEN
zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing
! cfl scale is function of the RK4 step
zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_scale ! Speed limit as a factor of dx / dt
IF( zspeed_new < zspeed ) THEN
zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed
zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction
pax = (zuveln - puvel0)/pdt
pay = (zvveln - pvvel0)/pdt
!
! print speeding ticket
IF (nn_verbose_level > 0) THEN
WRITE(numicb, 9200) 'icb speeding : ',kt, nknberg, zspeed, &
& pxi, pyj, zuo, zvo, zua, zva, zui, zvi
9200 FORMAT(a,i9,i6,f9.2,1x,4(1x,2f9.2))
END IF
!
CALL icb_dia_speed()
ENDIF
ENDIF
ENDIF
! ! check the speed and acceleration limits
IF (nn_verbose_level > 0) THEN
IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) &
WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) &
WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
ENDIF
!
END SUBROUTINE icb_accel
!!======================================================================
END MODULE icbdyn