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 branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90 @ 3339

Last change on this file since 3339 was 3339, checked in by sga, 12 years ago

NEMO branch dev_r3337_NOCS10_ICB: add new iceberg sub-directory ICB
File size: 17.3 KB
Line 
1MODULE icbdyn
2
3   !!======================================================================
4   !!                       ***  MODULE  icbdyn  ***
5   !! Ocean physics:  time stepping routine for iceberg tracking
6   !!======================================================================
7   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
8   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
9   !!            -    !                            Removal of mapping from another grid
10   !!            -    !  2011-04  (Alderson)       Split into separate modules
11   !!            -    !  2011-05  (Alderson)       Replace broken grounding routine
12   !!            -    !                            with one of Gurvan's suggestions (just like
13   !!            -    !                            the broken one)
14   !!----------------------------------------------------------------------
15   !!----------------------------------------------------------------------
16   !!   icb_init      : initialise
17   !!   icb_gen       : generate test icebergs
18   !!   icb_nam       : read iceberg namelist
19   !!----------------------------------------------------------------------
20   USE par_oce        ! NEMO parameters
21   USE dom_oce        ! NEMO ocean domain
22   USE phycst         ! NEMO physical constants
23
24   USE icb_oce        ! define iceberg arrays
25   USE icbutl         ! iceberg utility routines
26   USE icbdia         ! iceberg budget routines
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   evolve_icebergs  ! routine called in xxx.F90 module
32
33CONTAINS
34
35   SUBROUTINE evolve_icebergs()
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE evolve_icebergs  ***
38      !!
39      !! ** Purpose :   iceberg evolution.
40      !!
41      !! ** Method  : - blah blah
42      !!----------------------------------------------------------------------
43      !
44      REAL(wp)                        ::   uvel1 , vvel1 , u1, v1, ax1, ay1, xi1 , yj1
45      REAL(wp)                        ::   uvel2 , vvel2 , u2, v2, ax2, ay2, xi2 , yj2
46      REAL(wp)                        ::   uvel3 , vvel3 , u3, v3, ax3, ay3, xi3 , yj3
47      REAL(wp)                        ::   uvel4 , vvel4 , u4, v4, ax4, ay4, xi4 , yj4
48      REAL(wp)                        ::   uvel_n, vvel_n                  , xi_n, yj_n
49      REAL(wp)                        ::   zdt, zdt_2, zdt_6, e1, e2
50      LOGICAL                         ::   bounced
51      TYPE(iceberg), POINTER          ::   berg
52      TYPE(point)  , POINTER          ::   pt
53      !!----------------------------------------------------------------------
54
55      ! 4th order Runge-Kutta to solve:
56      !    d/dt X = V,  d/dt V = A
57      ! with I.C.'s:
58      !    X=X1 and V=V1
59      !
60      !                                    ; A1=A(X1,V1)
61      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
62      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
63      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
64      !
65      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
66      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
67
68      ! time steps
69      zdt   = berg_dt
70      zdt_2 = zdt * 0.5_wp
71      zdt_6 = zdt / 6._wp
72
73      berg => first_berg                    ! start from the first berg
74      !
75      DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
76         !
77         pt => berg%current_point
78
79         bounced = .false.
80
81
82         ! STEP 1 !
83         ! ====== !
84         xi1 = pt%xi   ;   uvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
85         yj1 = pt%yj   ;   vvel1 = pt%vvel
86
87
88         !                                         !**   A1 = A(X1,V1)
89         CALL accel(        xi1, e1, uvel1, uvel1, ax1,     &
90            &        berg , yj1, e2, vvel1, vvel1, ay1, zdt_2 )
91         !
92         u1 = uvel1 / e1                           !**   V1 in d(i,j)/dt
93         v1 = vvel1 / e2
94
95         ! STEP 2 !
96         ! ====== !
97         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
98         ! position using di/dt & djdt   !   V2  in m/s
99         xi2 = xi1 + zdt_2 * u1          ;   uvel2 = uvel1 + zdt_2 * ax1
100         yj2 = yj1 + zdt_2 * v1          ;   vvel2 = vvel1 + zdt_2 * ay1
101         !
102         CALL adjust_to_ground( xi2, xi1, u1, yj2, yj1, v1, bounced )
103
104         !                                         !**   A2 = A(X2,V2)
105         CALL accel(        xi2, e1, uvel2, uvel1, ax2,    &
106            &        berg , yj2, e2, vvel2, vvel1, ay2, zdt_2 )
107         !
108         u2 = uvel2 / e1                           !**   V2 in d(i,j)/dt
109         v2 = vvel2 / e2
110         !
111         ! STEP 3 !
112         ! ====== !
113         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
114         xi3  = xi1  + zdt_2 * u2   ;   uvel3 = uvel1 + zdt_2 * ax2
115         yj3  = yj1  + zdt_2 * v2   ;   vvel3 = vvel1 + zdt_2 * ay2
116         !
117         CALL adjust_to_ground( xi3, xi1, u3, yj3, yj1, v3, bounced )
118
119         !                                         !**   A3 = A(X3,V3)
120         CALL accel(        xi3, e1, uvel3, uvel1, ax3,    &
121            &        berg , yj3, e2, vvel3, vvel1, ay3, zdt )
122         !
123         u3 = uvel3 / e1                           !**   V3 in d(i,j)/dt
124         v3 = vvel3 / e2
125
126         ! STEP 4 !
127         ! ====== !
128         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
129         xi4 = xi1 + zdt * u3   ;   uvel4 = uvel1 + zdt * ax3
130         yj4 = yj1 + zdt * v3   ;   vvel4 = vvel1 + zdt * ay3
131
132         CALL adjust_to_ground( xi4, xi1, u4, yj4, yj1, v4, bounced )
133
134         !                                         !**   A4 = A(X4,V4)
135         CALL accel(        xi4, e1, uvel4, uvel1, ax4,    &
136            &        berg , yj4, e2, vvel4, vvel1, ay4, zdt )
137
138         u4 = uvel4 / e1                           !**   V4 in d(i,j)/dt
139         v4 = vvel4 / e2
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         xi_n   = pt%xi   + zdt_6 * (  u1  + 2.*(u2  + u3 ) + u4  )
146         yj_n   = pt%yj   + zdt_6 * (  v1  + 2.*(v2  + v3 ) + v4  )
147         uvel_n = pt%uvel + zdt_6 * (  ax1 + 2.*(ax2 + ax3) + ax4 )
148         vvel_n = pt%vvel + zdt_6 * (  ay1 + 2.*(ay2 + ay3) + ay4 )
149
150         CALL adjust_to_ground( xi_n, xi1, uvel_n, yj_n, yj1, vvel_n, bounced )
151
152         pt%uvel = uvel_n                        !** save in berg structure
153         pt%vvel = vvel_n
154         pt%xi   = xi_n
155         pt%yj   = yj_n
156
157         ! sga - update actual position
158         pt%lon  = bilin_x(glamt, pt%xi, pt%yj )
159         pt%lat  = bilin(gphit, pt%xi, pt%yj, 'T', 0, 0 )
160
161         berg => berg%next                         ! switch to the next berg
162         !
163      END DO                                  !==  end loop over all bergs  ==!
164      !
165   END SUBROUTINE evolve_icebergs
166
167
168   SUBROUTINE adjust_to_ground( pi, pi0, pu, pj, pj0, pv, bounced )
169      !!----------------------------------------------------------------------
170      !!                  ***  ROUTINE adjust_to_ground  ***
171      !!
172      !! ** Purpose :   iceberg grounding.
173      !!
174      !! ** Method  : - blah blah
175      !!----------------------------------------------------------------------
176      REAL(wp), INTENT(inout)         ::   pi , pj    ! current iceberg position
177      REAL(wp), INTENT(in   )         ::   pi0, pj0   ! previous iceberg position
178      REAL(wp), INTENT(inout)         ::   pu  , pv   ! current iceberg velocities
179      LOGICAL , INTENT(  out)         ::   bounced    ! bounced indicator
180      !
181      REAL(wp), PARAMETER ::   posn_eps = 0.05_wp     ! bouncing distance in (i,j) referential
182      !
183      INTEGER  ::   ii, ii0
184      INTEGER  ::   ij, ij0
185      INTEGER  ::   bounce_method
186      !!----------------------------------------------------------------------
187      !
188      bounced = .false.
189      bounce_method = 2
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 = ii0 - nimpp + 1
198      ij0 = ij0 - njmpp + 1
199      ii  = ii  - nimpp + 1
200      ij  = ij  - njmpp + 1
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      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
215      SELECT CASE ( bounce_method )
216         CASE ( 1 )
217            pi = pi0
218            pj = pj0
219            pu = 0._wp
220            pv = 0._wp
221         CASE ( 2 )
222            IF( ii0 /= ii ) THEN
223               pi = pi0                   ! return back to the initial position
224               pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
225            ENDIF
226            IF( ij0 /= ij ) THEN
227               pj = pj0                   ! return back to the initial position
228               pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
229            ENDIF
230      END SELECT
231      !
232   END SUBROUTINE adjust_to_ground
233
234
235   SUBROUTINE accel(       xi, e1, uvel, uvel0, ax,                &
236      &             berg , yj, e2, vvel, vvel0, ay, dt )
237      !!----------------------------------------------------------------------
238      !!                  ***  ROUTINE accel  ***
239      !!
240      !! ** Purpose :   compute the iceberg acceleration.
241      !!
242      !! ** Method  : - blah blah
243      !!----------------------------------------------------------------------
244      TYPE(iceberg ), POINTER           ::   berg
245      REAL(wp), INTENT(in   )           ::   xi   , yj      ! berg position in (i,j) referential
246      REAL(wp), INTENT(in   )           ::   uvel , vvel    ! berg velocity [m/s]
247      REAL(wp), INTENT(in   )           ::   uvel0, vvel0   ! initial berg velocity [m/s]
248      REAL(wp), INTENT(  out)           ::   e1, e2         ! horizontal scale factor at (xi,yj)
249      REAL(wp), INTENT(inout)           ::   ax, ay         ! berg acceleration
250      REAL(wp), INTENT(in   )           ::   dt             ! berg time step
251      !
252      REAL(wp), PARAMETER ::   alpha     = 0._wp      !
253      REAL(wp), PARAMETER ::   beta      = 1._wp      !
254      REAL(wp), PARAMETER ::   vel_lim   =15._wp      ! max allowed berg speed
255      REAL(wp), PARAMETER ::   accel_lim = 1.e-2_wp   ! max allowed berg acceleration
256      REAL(wp), PARAMETER ::   Cr0       = 0.06_wp    !
257      !
258      INTEGER  ::   itloop
259      REAL(wp) ::   uo, vo, ui, vi, ua, va, uwave, vwave, ssh_x, ssh_y, sst, cn, hi
260      REAL(wp) ::   zff, T, D, W, L, M, F
261      REAL(wp) ::   drag_ocn, drag_atm, drag_ice, wave_rad
262      REAL(wp) ::   c_ocn, c_atm, c_ice
263      REAL(wp) ::   ampl, wmod, Cr, Lwavelength, Lcutoff, Ltop
264      REAL(wp) ::   lambda, detA, A11, A12, axe, aye, D_hi
265      REAL(wp) ::   uveln, vveln, us, vs, speed, loc_dx, new_speed
266      !!----------------------------------------------------------------------
267
268      ! Interpolate gridded fields to berg
269      knberg = berg%number(1)
270      CALL interp_flds( xi, e1, uo, ui, ua, ssh_x,                     &
271         &              yj, e2, vo, vi, va, ssh_y, sst, cn, hi, zff )
272
273      M = berg%current_point%mass
274      T = berg%current_point%thickness                           ! total thickness
275      D = ( rn_rho_bergs / rho_seawater ) * T   ! draught (keel depth)
276      F = T - D                                    ! freeboard
277      W = berg%current_point%width
278      L = berg%current_point%length
279
280      hi   = MIN( hi   , D    )
281      D_hi = MAX( 0._wp, D-hi )
282
283      ! Wave radiation
284      uwave = ua - uo   ;   vwave = va - vo     ! Use wind speed rel. to ocean for wave model
285      wmod  = uwave*uwave + vwave*vwave         ! The wave amplitude and length depend on the  current;
286      !                                         ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
287      ampl        = 0.5 * 0.02025 * wmod        ! This is "a", the wave amplitude
288      Lwavelength =       0.32    * wmod        ! Surface wave length fitted to data in table at
289      !                                         ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
290      Lcutoff     = 0.125 * Lwavelength
291      Ltop        = 0.25  * Lwavelength
292      Cr          = Cr0 * MIN(  MAX( 0., (L-Lcutoff) / ((Ltop-Lcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
293      !                                         ! fitted to graph from Carrieres et al.,  POAC Drift Model.
294      wave_rad    = 0.5 * rho_seawater / M * Cr * grav * ampl * MIN( ampl,F ) * (2.*W*L) / (W+L)
295      wmod        = SQRT( ua*ua + va*va )       ! Wind speed
296      IF( wmod /= 0._wp ) THEN
297         uwave = ua/wmod ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
298         vwave = va/wmod
299      ELSE
300         uwave = 0.   ;    vwave=0.   ;    wave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
301      ENDIF
302
303      ! Weighted drag coefficients
304      c_ocn = rho_seawater / M * (0.5*Cd_wv*W*(D_hi)+Cd_wh*W*L)
305      c_atm = rho_air      / M * (0.5*Cd_av*W*F     +Cd_ah*W*L)
306      c_ice = rho_ice      / M * (0.5*Cd_iv*W*hi              )
307      IF( abs(ui) + abs(vi) == 0._wp )   c_ice = 0._wp
308
309      uveln = uvel   ;   vveln = vvel ! Copy starting uvel, vvel
310      !
311      DO itloop = 1, 2  ! Iterate on drag coefficients
312         !
313         us = 0.5 * ( uveln + uvel )
314         vs = 0.5 * ( vveln + vvel )
315         drag_ocn = c_ocn * SQRT( (us-uo)*(us-uo) + (vs-vo)*(vs-vo) )
316         drag_atm = c_atm * SQRT( (us-ua)*(us-ua) + (vs-va)*(vs-va) )
317         drag_ice = c_ice * SQRT( (us-ui)*(us-ui) + (vs-vi)*(vs-vi) )
318         !
319         ! Explicit accelerations
320         !axe= zff*vvel -grav*ssh_x +wave_rad*uwave &
321         !    -drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)
322         !aye=-zff*uvel -grav*ssh_y +wave_rad*vwave &
323         !    -drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)
324         axe = -grav * ssh_x + wave_rad * uwave
325         aye = -grav * ssh_y + wave_rad * vwave
326         IF( alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
327            axe = axe + zff*vvel0
328            aye = aye - zff*uvel0
329         else
330            axe=axe+zff*vvel
331            aye=aye-zff*uvel
332         endif
333         if (beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest
334            axe=axe-drag_ocn*(uvel0-uo) -drag_atm*(uvel0-ua) -drag_ice*(uvel0-ui)
335            aye=aye-drag_ocn*(vvel0-vo) -drag_atm*(vvel0-va) -drag_ice*(vvel0-vi)
336         else
337            axe=axe-drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)
338            aye=aye-drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)
339         endif
340
341         ! Solve for implicit accelerations
342         IF( alpha + beta > 0._wp ) THEN
343            lambda = drag_ocn + drag_atm + drag_ice
344            A11 = 1.+beta*dt*lambda
345            A12 = alpha*dt*zff
346            detA = 1._wp / ( A11*A11 + A12*A12 )
347            ax = detA * ( A11*axe+A12*aye )
348            ay = detA * ( A11*aye-A12*axe )
349         else
350            ax=axe   ;   ay=aye
351         endif
352
353         uveln = uvel0 + dt*ax
354         vveln = vvel0 + dt*ay
355         !
356      END DO      ! itloop
357
358      IF( rn_speed_limit > 0.) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
359         speed = SQRT( uveln*uveln + vveln*vveln )    ! Speed of berg
360         IF( speed > 0.) THEN
361            loc_dx = MIN( e1, e2 )                          ! minimum grid spacing
362            new_speed = loc_dx / dt * rn_speed_limit     ! Speed limit as a factor of dx / dt
363            IF( new_speed < speed ) THEN
364               uveln = uveln * ( new_speed / speed )        ! Scale velocity to reduce speed
365               vveln = vveln * ( new_speed / speed )        ! without changing the direction
366               CALL speed_budget()
367            ENDIF
368         ENDIF
369      ENDIF
370      !                                      ! check the speed and acceleration limits
371      IF( ABS( uveln ) > vel_lim   .OR. ABS( vveln ) > vel_lim   )   &
372         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
373      IF( ABS( ax    ) > accel_lim .OR. ABS( ay    ) > accel_lim )   &
374         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
375      !
376   END SUBROUTINE accel
377
378END MODULE icbdyn
Note: See TracBrowser for help on using the repository browser.