source: NEMO/branches/UKMO/NEMO_4.0.4_bouncing_icebergs/src/OCE/ICB/icbdyn.F90 @ 14711

Last change on this file since 14711 was 14711, checked in by dancopsey, 2 months ago

Add option for icebergs to bounce away from the coast

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