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/2020/tickets_icb_1900/src/OCE/ICB – NEMO

source: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90 @ 13374

Last change on this file since 13374 was 13374, checked in by mathiot, 4 years ago

ticket #1900: fix issue about mask management in icb_utl_bilin_3d_h; add option to ground icb if icb bottom lvl hit oce bottom lvl

  • Property svn:keywords set to Id
File size: 19.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   USE in_out_manager                      ! IO parameters
17   !
18   USE icb_oce        ! define iceberg arrays
19   USE icbutl         ! iceberg utility routines
20   USE icbdia         ! iceberg budget routines
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module
26
27   !!----------------------------------------------------------------------
28   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
29   !! $Id$
30   !! Software governed by the CeCILL license (see ./LICENSE)
31   !!----------------------------------------------------------------------
32CONTAINS
33
34   SUBROUTINE icb_dyn( kt )
35      !!----------------------------------------------------------------------
36      !!                  ***  ROUTINE icb_dyn  ***
37      !!
38      !! ** Purpose :   iceberg evolution.
39      !!
40      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
41      !!----------------------------------------------------------------------
42      INTEGER, INTENT(in) ::   kt   !
43      !
44      LOGICAL  ::   ll_bounced
45      REAL(wp) ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
46      REAL(wp) ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
47      REAL(wp) ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
48      REAL(wp) ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
49      REAL(wp) ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
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( berg, 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( berg, 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( berg, 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( berg, zxi_n, zxi1, zuvel_n,   &
152            &                   zyj_n, zyj1, zvvel_n, ll_bounced )
153
154         pt%uvel = zuvel_n                        !** save in berg structure
155         pt%vvel = zvvel_n
156         pt%xi   = zxi_n
157         pt%yj   = zyj_n
158
159         berg => berg%next                         ! switch to the next berg
160         !
161      END DO                                  !==  end loop over all bergs  ==!
162      !
163   END SUBROUTINE icb_dyn
164
165
166   SUBROUTINE icb_ground( berg, pi, pi0, pu,   &
167      &                         pj, pj0, pv, ld_bounced )
168      !!----------------------------------------------------------------------
169      !!                  ***  ROUTINE icb_ground  ***
170      !!
171      !! ** Purpose :   iceberg grounding.
172      !!
173      !! ** Method  : - adjust velocity and then put iceberg back to start position
174      !!                NB two possibilities available one of which is hard-coded here
175      !!----------------------------------------------------------------------
176      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
177      !
178      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
179      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
180      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
181      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
182      !
183      INTEGER  ::   ii, ii0
184      INTEGER  ::   ij, ij0
185      INTEGER  ::   ikb
186      INTEGER  ::   ibounce_method
187      !
188      REAL(wp) :: zD 
189      REAL(wp), DIMENSION(jpk) :: ze3t
190      !!----------------------------------------------------------------------
191      !
192      ld_bounced = .FALSE.
193      !
194      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
195      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
196      !
197      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
198      !
199      ! map into current processor
200      ii0 = mi1( ii0 )
201      ij0 = mj1( ij0 )
202      ii  = mi1( ii  )
203      ij  = mj1( ij  )
204      !
205      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0
206      IF ( ln_M2016 .AND. ln_icb_grd ) THEN
207         !
208         ! draught (keel depth)
209         zD = rho_berg_1_oce * berg%current_point%thickness
210         !
211         ! interpol needed data
212         CALL icb_utl_interp( pi, pj, pe3t=ze3t )
213         !
214         !compute bottom level
215         CALL icb_utl_getkb( ikb, ze3t, zD )
216         !
217         ! berg reach a new t-cell, but an ocean one
218         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN
219         !
220      ELSE
221         IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
222      END IF
223      !
224      ! From here, berg have reach land: treat grounding/bouncing
225      ! -------------------------------
226      ld_bounced = .TRUE.
227
228      !! not obvious what should happen now
229      !! if berg tries to enter a land box, the only location we can return it to is the start
230      !! position (pi0,pj0), since it has to be in a wet box to do any melting;
231      !! first option is simply to set whole velocity to zero and move back to start point
232      !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
233      !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast
234
235      ibounce_method = 2
236      SELECT CASE ( ibounce_method )
237      CASE ( 1 )
238         pi = pi0
239         pj = pj0
240         pu = 0._wp
241         pv = 0._wp
242      CASE ( 2 )
243         IF( ii0 /= ii ) THEN
244            pi = pi0                   ! return back to the initial position
245            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
246         ENDIF
247         IF( ij0 /= ij ) THEN
248            pj = pj0                   ! return back to the initial position
249            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
250         ENDIF
251      END SELECT
252      !
253   END SUBROUTINE icb_ground
254
255
256   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
257      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
258      !!----------------------------------------------------------------------
259      !!                  ***  ROUTINE icb_accel  ***
260      !!
261      !! ** Purpose :   compute the iceberg acceleration.
262      !!
263      !! ** Method  : - sum the terms in the momentum budget
264      !!----------------------------------------------------------------------
265      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
266      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
267      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
268      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
269      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
270      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
271      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
272      !
273      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
274      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
275      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
276      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
277      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
278      !
279      INTEGER  ::   itloop, ikb, jk
280      REAL(wp) ::   zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi
281      REAL(wp) ::   zvo, zssv, zvi, zva, zvwave, zssh_y
282      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
283      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
284      REAL(wp) ::   z_ocn, z_atm, z_ice, zdep
285      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
286      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
287      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
288      REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw
289      !!----------------------------------------------------------------------
290
291      ! Interpolate gridded fields to berg
292      nknberg = berg%number(1)
293      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,     &   ! scale factor
294         &                 pssu=zssu, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
295         &                 pssv=zssv, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
296         &                 pssh_i=zssh_x, pssh_j=zssh_y,   &   ! ssh gradient
297         &                 phi=zhi, pff=zff)                   ! ice thickness and coriolis
298
299      zM = berg%current_point%mass
300      zT = berg%current_point%thickness               ! total thickness
301      zD = rho_berg_1_oce * zT                        ! draught (keel depth)
302      zF = zT - zD                                    ! freeboard
303      zW = berg%current_point%width
304      zL = berg%current_point%length
305
306      zhi   = MIN( zhi   , zD    )
307      zD_hi = MAX( 0._wp, zD-zhi )
308 
309     ! Wave radiation
310      zuwave = zua - zssu   ;   zvwave = zva - zssv   ! Use wind speed rel. to ocean for wave model
311      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
312      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
313      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
314      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
315      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
316      zLcutoff     = 0.125 * zLwavelength
317      zLtop        = 0.25  * zLwavelength
318      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
319      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
320      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
321      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
322      IF( zwmod /= 0._wp ) THEN
323         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
324         zvwave = zva/zwmod
325      ELSE
326         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
327      ENDIF
328
329      ! Weighted drag coefficients
330      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
331      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
332      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
333      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
334
335      ! lateral velocities
336      ! default ssu and ssv
337      ! ln_M2016: mean velocity along the profile
338      IF ( ln_M2016 ) THEN
339         ! interpol needed data
340         CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )   ! 3d velocities
341       
342         !compute bottom level
343         CALL icb_utl_getkb( ikb, ze3t, zD )
344         
345         ! compute mean velocity
346         CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb)
347         CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb)
348      ELSE
349         zuo = zssu
350         zvo = zssv
351      END IF
352
353      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
354      !
355      DO itloop = 1, 2  ! Iterate on drag coefficients
356         !
357         zus = 0.5 * ( zuveln + puvel )
358         zvs = 0.5 * ( zvveln + pvvel )
359         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
360         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
361         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
362         !
363         ! Explicit accelerations
364         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
365         !    -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
366         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
367         !    -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
368         zaxe = -grav * zssh_x + zwave_rad * zuwave
369         zaye = -grav * zssh_y + zwave_rad * zvwave
370         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
371            zaxe = zaxe + zff*pvvel0
372            zaye = zaye - zff*puvel0
373         ELSE
374            zaxe = zaxe + zff*pvvel
375            zaye = zaye - zff*puvel
376         ENDIF
377         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
378            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
379            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
380         ELSE
381            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
382            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
383         ENDIF
384
385         ! Solve for implicit accelerations
386         IF( pp_alpha + pp_beta > 0._wp ) THEN
387            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
388            zA11    = 1._wp + pp_beta *pdt*zlambda
389            zA12    =         pp_alpha*pdt*zff
390            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
391            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
392            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
393         ELSE
394            pax = zaxe   ;   pay = zaye
395         ENDIF
396
397         zuveln = puvel0 + pdt*pax
398         zvveln = pvvel0 + pdt*pay
399         !
400      END DO      ! itloop
401
402      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
403         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
404         IF( zspeed > 0._wp ) THEN
405            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
406            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
407            IF( zspeed_new < zspeed ) THEN
408               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
409               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
410               CALL icb_dia_speed()
411            ENDIF
412         ENDIF
413      ENDIF
414      !                                      ! check the speed and acceleration limits
415      IF (nn_verbose_level > 0) THEN
416         IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
417            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
418         IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
419            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
420      ENDIF
421      !
422   END SUBROUTINE icb_accel
423
424   !!======================================================================
425END MODULE icbdyn
Note: See TracBrowser for help on using the repository browser.