New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icbdyn.F90 in NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdyn.F90 @ 14102

Last change on this file since 14102 was 14102, checked in by davestorkey, 4 years ago

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