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/2019/fix_ticket2238_solution2/src/OCE/ICB – NEMO

source: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbdyn.F90 @ 10693

Last change on this file since 10693 was 10679, checked in by mathiot, 5 years ago

branch for solution 2 of ticket #2238
  • Property svn:keywords set to Id
File size: 23.7 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   USE icblbc         ! iceberg: lateral boundary routines (including mpp)
21   !
22   USE in_out_manager
23   USE lib_mpp        ! massively parallel library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module
29
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
33   !! Software governed by the CeCILL license (see ./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE icb_dyn( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE icb_dyn  ***
40      !!
41      !! ** Purpose :   iceberg evolution.
42      !!
43      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
44      !!----------------------------------------------------------------------
45      INTEGER, INTENT(in) ::   kt   !
46      !
47      LOGICAL  ::   ll_bounced
48      REAL(wp) ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
49      REAL(wp) ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
50      REAL(wp) ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
51      REAL(wp) ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
52      REAL(wp) ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
53      REAL(wp) ::   zdt, zdt_2, zdt_6, ze1, ze2
54      TYPE(iceberg), POINTER ::   berg
55      TYPE(point)  , POINTER ::   pt
56      !!----------------------------------------------------------------------
57      !
58      ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
59      !                    with I.C.'s:   X=X1 and V=V1
60      !
61      !                                    ; A1=A(X1,V1)
62      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
63      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
64      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
65      !
66      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
67      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
68
69      ! time steps
70      zdt   = berg_dt
71      zdt_2 = zdt * 0.5_wp
72      zdt_6 = zdt / 6._wp
73
74      IF( ASSOCIATED(first_berg) ) THEN
75         berg => first_berg                    ! start from the first berg
76         !
77         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
78            !
79            pt => berg%current_point
80            !
81            ! STEP 1 !
82            ! ====== !
83            zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
84            zyj1 = pt%yj   ;   zvvel1 = pt%vvel
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            ! check ground iceberg
101            CALL icb_ground( zxi2, zxi1, zu1,   &
102               &             zyj2, zyj1, zv1, ll_bounced )
103            !
104            ! SAVE intermediaite properties
105            pt%xRK1 = (/ zxi1, zuvel1, zu1, zax1 /); pt%yRK1 = (/ zyj1, zvvel1, zv1, zay1 /); 
106            pt%xRK2 = (/ zxi2, zuvel2, 0.0, 0.0  /); pt%yRK2 = (/ zyj2, zvvel2, 0.0, 0.0  /); 
107            pt%xi = zxi2 ; pt%yj = zyj2 ;
108            !
109            berg => berg%next                         ! switch to the next berg
110            !
111         END DO
112      END IF
113      !
114      ! EXCHANGE BERG (in case first step fall out of the inner domain)
115      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
116      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
117      ENDIF
118      !
119      IF( ASSOCIATED(first_berg) ) THEN
120         berg => first_berg                    ! start from the first berg
121         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
122            !
123            pt => berg%current_point
124            !
125            ! LOAD intermediaite properties
126            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
127            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 
128            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
129            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 
130            !
131            !                                         !**   A2 = A(X2,V2)
132            CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
133               &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
134            !
135            zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
136            zv2 = zvvel2 / ze2
137            !
138            ! STEP 3 !
139            ! ====== !
140            !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
141            zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
142            zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
143            !
144            CALL icb_ground( zxi3, zxi1, zu3,   &
145               &             zyj3, zyj1, zv3, ll_bounced )
146            !
147            ! SAVE intermediaite properties
148            pt%xRK2(3:4) =               (/ zu2, zax2 /); pt%yRK2(3:4) =               (/ zv2, zay2 /);
149            pt%xRK3      = (/ zxi3, zuvel3, 0.0, 0.0  /); pt%yRK3      = (/ zyj3, zvvel3, 0.0, 0.0  /); 
150            pt%xi = zxi3 ; pt%yj = zyj3 ;
151            !
152            berg => berg%next                         ! switch to the next berg
153            !
154         END DO
155      END IF
156      !
157      ! EXCHANGE BERG (in case first step fall out of the inner domain)
158      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
159      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
160      ENDIF
161      !
162      IF( ASSOCIATED(first_berg) ) THEN
163         berg => first_berg                    ! start from the first berg
164         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
165            !
166            pt => berg%current_point
167            !
168            ! LOAD intermediaite properties
169            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
170            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 
171            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
172            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 
173            !
174            !                                         !**   A3 = A(X3,V3)
175            CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
176               &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
177            !
178            zu3 = zuvel3 / ze1                        !**   V3 in d(i,j)/dt
179            zv3 = zvvel3 / ze2
180            !
181            ! STEP 4 !
182            ! ====== !
183            !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
184            zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
185            zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
186            !
187            CALL icb_ground( zxi4, zxi1, zu4,   &
188               &             zyj4, zyj1, zv4, ll_bounced )
189            !
190            ! SAVE intermediaite properties
191            pt%xRK3(3:4) =               (/ zu3, zax3 /); pt%yRK3(3:4) =               (/ zv3, zay3 /); 
192            pt%xRK4      = (/ zxi4, zuvel4, 0.0, 0.0  /); pt%yRK4      = (/ zyj4, zvvel4, 0.0, 0.0  /); 
193            pt%xi = zxi4 ; pt%yj = zyj4 ;
194            !
195            berg => berg%next                         ! switch to the next berg
196            !
197         END DO
198      END IF
199      !
200      ! EXCHANGE BERG (in case first step fall out of the inner domain)
201      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
202      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
203      ENDIF
204      !
205      IF( ASSOCIATED(first_berg) ) THEN
206         berg => first_berg                    ! start from the first berg
207         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
208            !
209            pt => berg%current_point
210            !
211            ! LOAD intermediaite properties
212            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
213            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 
214            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 
215            zxi4 = pt%xRK4(1) ; zuvel4 = pt%xRK4(2) ; zu4 = pt%xRK4(3) ; zax4 = pt%xRK4(4) 
216            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
217            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 
218            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 
219            zyj4 = pt%yRK4(1) ; zvvel4 = pt%yRK4(2) ; zv4 = pt%yRK4(3) ; zay4 = pt%yRK4(4) 
220            !
221            !                                         !**   A4 = A(X4,V4)
222            CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
223               &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
224            !
225            zu4 = zuvel4 / ze1                        !**   V4 in d(i,j)/dt
226            zv4 = zvvel4 / ze2
227            !
228            ! FINAL STEP !
229            ! ========== !
230            !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
231            !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
232            zxi_n   = zxi1   + zdt_6 * (  zu1  + 2._wp*(zu2  + zu3 ) + zu4  )
233            zyj_n   = zyj1   + zdt_6 * (  zv1  + 2._wp*(zv2  + zv3 ) + zv4  )
234            zuvel_n = zuvel1 + zdt_6 * (  zax1 + 2._wp*(zax2 + zax3) + zax4 )
235            zvvel_n = zvvel1 + zdt_6 * (  zay1 + 2._wp*(zay2 + zay3) + zay4 )
236            !
237            CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
238               &             zyj_n, zyj1, zvvel_n, ll_bounced )
239            !
240            !                                        !** save in berg structure
241            pt%xi   = zxi_n   ; pt%yj   = zyj_n ;
242            pt%uvel = zuvel_n ; pt%vvel = zvvel_n
243            !
244            berg => berg%next                        ! switch to the next berg
245            !
246         END DO
247      END IF
248      !
249      ! EXCHANGE BERG (in case first step fall out of the inner domain)
250      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
251      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
252      ENDIF
253      !
254      IF( ASSOCIATED(first_berg) ) THEN
255         berg => first_berg                    ! start from the first berg
256         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
257            !
258            pt => berg%current_point
259            !
260            ! update actual position
261            pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
262            pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
263            !
264            berg => berg%next                         ! switch to the next berg
265            !
266         END DO                                  !==  end loop over all bergs  ==!
267      END IF
268      !
269   END SUBROUTINE icb_dyn
270
271
272   SUBROUTINE icb_ground( pi, pi0, pu,   &
273      &                   pj, pj0, pv, ld_bounced )
274      !!----------------------------------------------------------------------
275      !!                  ***  ROUTINE icb_ground  ***
276      !!
277      !! ** Purpose :   iceberg grounding.
278      !!
279      !! ** Method  : - adjust velocity and then put iceberg back to start position
280      !!                NB two possibilities available one of which is hard-coded here
281      !!----------------------------------------------------------------------
282      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
283      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
284      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
285      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
286      !
287      INTEGER  ::   ii, iig, iig0
288      INTEGER  ::   ij, ijg, ijg0
289      INTEGER  ::   ibounce_method
290      !!----------------------------------------------------------------------
291      !
292      ! (PM) is this variable used somewhere ?
293      ld_bounced = .FALSE.
294      !
295      ! work with global position in case pi and pi0 are on different domains
296      iig0 = INT( pi0 + 0.5_wp )   ;   ijg0 = INT( pj0 + 0.5_wp )       ! initial gridpoint position (T-cell)
297      iig  = INT( pi  + 0.5_wp )   ;   ijg  = INT( pj  + 0.5_wp )       ! current     -         -
298      !
299      IF( iig == iig0  .AND.  ijg == ijg0  )   RETURN           ! berg remains in the same cell
300      !
301      ! map into current processor
302      ii  = mi1( iig  )
303      ij  = mj1( ijg  )
304      !
305      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
306      !
307      ! From here, berg have reach land: treat grounding/bouncing
308      ! -------------------------------
309      ld_bounced = .TRUE.
310
311      !! not obvious what should happen now
312      !! if berg tries to enter a land box, the only location we can return it to is the start
313      !! position (pi0,pj0), since it has to be in a wet box to do any melting;
314      !! first option is simply to set whole velocity to zero and move back to start point
315      !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
316      !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast
317
318      ibounce_method = 2
319      SELECT CASE ( ibounce_method )
320      CASE ( 1 )
321         pi = pi0
322         pj = pj0
323         pu = 0._wp
324         pv = 0._wp
325      CASE ( 2 )
326         IF( iig0 /= iig ) THEN
327            pi = pi0                   ! return back to the initial position
328            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
329         ENDIF
330         IF( ijg0 /= ijg ) THEN
331            pj = pj0                   ! return back to the initial position
332            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
333         ENDIF
334      END SELECT
335      !
336   END SUBROUTINE icb_ground
337
338
339   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
340      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
341      !!----------------------------------------------------------------------
342      !!                  ***  ROUTINE icb_accel  ***
343      !!
344      !! ** Purpose :   compute the iceberg acceleration.
345      !!
346      !! ** Method  : - sum the terms in the momentum budget
347      !!----------------------------------------------------------------------
348      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
349      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
350      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
351      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
352      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
353      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
354      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
355      !
356      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
357      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
358      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
359      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
360      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
361      !
362      INTEGER  ::   itloop
363      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
364      REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
365      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
366      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
367      REAL(wp) ::   z_ocn, z_atm, z_ice
368      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
369      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
370      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
371      !!----------------------------------------------------------------------
372
373      ! Interpolate gridded fields to berg
374      nknberg = berg%number(1)
375      CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     &
376         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff )
377
378      zM = berg%current_point%mass
379      zT = berg%current_point%thickness               ! total thickness
380      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
381      zF = zT - zD                                    ! freeboard
382      zW = berg%current_point%width
383      zL = berg%current_point%length
384
385      zhi   = MIN( zhi   , zD    )
386      zD_hi = MAX( 0._wp, zD-zhi )
387
388      ! Wave radiation
389      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
390      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
391      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
392      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
393      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
394      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
395      zLcutoff     = 0.125 * zLwavelength
396      zLtop        = 0.25  * zLwavelength
397      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
398      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
399      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
400      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
401      IF( zwmod /= 0._wp ) THEN
402         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
403         zvwave = zva/zwmod
404      ELSE
405         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
406      ENDIF
407
408      ! Weighted drag coefficients
409      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
410      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
411      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
412      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
413
414      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
415      !
416      DO itloop = 1, 2  ! Iterate on drag coefficients
417         !
418         zus = 0.5 * ( zuveln + puvel )
419         zvs = 0.5 * ( zvveln + pvvel )
420         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
421         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
422         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
423         !
424         ! Explicit accelerations
425         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
426         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
427         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
428         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
429         zaxe = -grav * zssh_x + zwave_rad * zuwave
430         zaye = -grav * zssh_y + zwave_rad * zvwave
431         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
432            zaxe = zaxe + zff*pvvel0
433            zaye = zaye - zff*puvel0
434         ELSE
435            zaxe = zaxe + zff*pvvel
436            zaye = zaye - zff*puvel
437         ENDIF
438         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
439            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
440            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
441         ELSE
442            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
443            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
444         ENDIF
445
446         ! Solve for implicit accelerations
447         IF( pp_alpha + pp_beta > 0._wp ) THEN
448            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
449            zA11    = 1._wp + pp_beta *pdt*zlambda
450            zA12    =         pp_alpha*pdt*zff
451            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
452            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
453            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
454         ELSE
455            pax = zaxe   ;   pay = zaye
456         ENDIF
457
458         zuveln = puvel0 + pdt*pax
459         zvveln = pvvel0 + pdt*pay
460         !
461      END DO      ! itloop
462
463      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
464         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
465         IF( zspeed > 0._wp ) THEN
466            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
467            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
468            IF( zspeed_new < zspeed ) THEN
469               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
470               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
471               CALL icb_dia_speed()
472            ENDIF
473         ENDIF
474      ENDIF
475      !                                      ! check the speed and acceleration limits
476      IF (nn_verbose_level > 0) THEN
477         IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
478            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
479         IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
480            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
481      ENDIF
482      !
483   END SUBROUTINE icb_accel
484
485   !!======================================================================
486END MODULE icbdyn
Note: See TracBrowser for help on using the repository browser.