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

source: NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbdyn.F90 @ 13265

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

ticket #2494 and #2375: ticket #2494 changes and part of ticket #2375 (lbc_icb_lnk on extended variable at T point removed), ff_e initialisation and lbc_lnk move in the icbini.F90; tests on possible useless lbc_lnk in icbclv not yet done
  • Property svn:keywords set to Id
File size: 18.2 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( 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         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( 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      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) ::   ld_bounced   ! bounced indicator
180      !
181      INTEGER  ::   ii, ii0
182      INTEGER  ::   ij, ij0
183      INTEGER  ::   ibounce_method
184      !!----------------------------------------------------------------------
185      !
186      ld_bounced = .FALSE.
187      !
188      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
189      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
190      !
191      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
192      !
193      ! map into current processor
194      ii0 = mi1( ii0 )
195      ij0 = mj1( ij0 )
196      ii  = mi1( ii  )
197      ij  = mj1( ij  )
198      !
199      IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
200      !
201      ! From here, berg have reach land: treat grounding/bouncing
202      ! -------------------------------
203      ld_bounced = .TRUE.
204
205      !! not obvious what should happen now
206      !! if berg tries to enter a land box, the only location we can return it to is the start
207      !! position (pi0,pj0), since it has to be in a wet box to do any melting;
208      !! first option is simply to set whole velocity to zero and move back to start point
209      !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction
210      !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast
211
212      ibounce_method = 2
213      SELECT CASE ( ibounce_method )
214      CASE ( 1 )
215         pi = pi0
216         pj = pj0
217         pu = 0._wp
218         pv = 0._wp
219      CASE ( 2 )
220         IF( ii0 /= ii ) THEN
221            pi = pi0                   ! return back to the initial position
222            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
223         ENDIF
224         IF( ij0 /= ij ) THEN
225            pj = pj0                   ! return back to the initial position
226            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
227         ENDIF
228      END SELECT
229      !
230   END SUBROUTINE icb_ground
231
232
233   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
234      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
235      !!----------------------------------------------------------------------
236      !!                  ***  ROUTINE icb_accel  ***
237      !!
238      !! ** Purpose :   compute the iceberg acceleration.
239      !!
240      !! ** Method  : - sum the terms in the momentum budget
241      !!----------------------------------------------------------------------
242      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
243      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
244      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
245      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
246      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
247      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
248      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
249      !
250      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
251      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
252      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
253      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
254      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
255      !
256      INTEGER  ::   itloop
257      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
258      REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
259      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
260      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
261      REAL(wp) ::   z_ocn, z_atm, z_ice
262      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
263      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
264      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
265      !!----------------------------------------------------------------------
266
267      ! Interpolate gridded fields to berg
268      nknberg = berg%number(1)
269      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,   &   ! scale factor
270         &                 puo=zuo, pui=zui, pua=zua,    &   ! oce/ice/atm velocities
271         &                 pvo=zvo, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities
272         &                 pssh_i=zssh_x, pssh_j=zssh_y, &   ! ssh gradient
273         &                 phi=zhi, pff=zff )               ! ice thickness and coriolis
274
275      IF (lwp) THEN
276         WRITE(numout,*) 'icbdyn ', pxi, pyj 
277         WRITE(numout,*) pe1, zuo, zui, zua, zssh_x
278         WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y
279         WRITE(numout,*) zhi, zff 
280         WRITE(numout,*) ''
281      END IF
282
283
284      zM = berg%current_point%mass
285      zT = berg%current_point%thickness               ! total thickness
286      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
287      zF = zT - zD                                    ! freeboard
288      zW = berg%current_point%width
289      zL = berg%current_point%length
290
291      zhi   = MIN( zhi   , zD    )
292      zD_hi = MAX( 0._wp, zD-zhi )
293
294      ! Wave radiation
295      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
296      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
297      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
298      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
299      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
300      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
301      zLcutoff     = 0.125 * zLwavelength
302      zLtop        = 0.25  * zLwavelength
303      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
304      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
305      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
306      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
307      IF( zwmod /= 0._wp ) THEN
308         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
309         zvwave = zva/zwmod
310      ELSE
311         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
312      ENDIF
313
314      ! Weighted drag coefficients
315      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
316      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
317      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
318      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
319
320      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
321      !
322      DO itloop = 1, 2  ! Iterate on drag coefficients
323         !
324         zus = 0.5 * ( zuveln + puvel )
325         zvs = 0.5 * ( zvveln + pvvel )
326         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
327         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
328         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
329         !
330         ! Explicit accelerations
331         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
332         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
333         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
334         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
335         zaxe = -grav * zssh_x + zwave_rad * zuwave
336         zaye = -grav * zssh_y + zwave_rad * zvwave
337         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
338            zaxe = zaxe + zff*pvvel0
339            zaye = zaye - zff*puvel0
340         ELSE
341            zaxe = zaxe + zff*pvvel
342            zaye = zaye - zff*puvel
343         ENDIF
344         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
345            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
346            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
347         ELSE
348            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
349            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
350         ENDIF
351
352         ! Solve for implicit accelerations
353         IF( pp_alpha + pp_beta > 0._wp ) THEN
354            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
355            zA11    = 1._wp + pp_beta *pdt*zlambda
356            zA12    =         pp_alpha*pdt*zff
357            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
358            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
359            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
360         ELSE
361            pax = zaxe   ;   pay = zaye
362         ENDIF
363
364         zuveln = puvel0 + pdt*pax
365         zvveln = pvvel0 + pdt*pay
366         !
367      END DO      ! itloop
368
369      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
370         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
371         IF( zspeed > 0._wp ) THEN
372            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
373            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
374            IF( zspeed_new < zspeed ) THEN
375               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
376               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
377               CALL icb_dia_speed()
378            ENDIF
379         ENDIF
380      ENDIF
381      !                                      ! check the speed and acceleration limits
382      IF (nn_verbose_level > 0) THEN
383         IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
384            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
385         IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
386            WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
387      ENDIF
388      !
389   END SUBROUTINE icb_accel
390
391   !!======================================================================
392END MODULE icbdyn
Note: See TracBrowser for help on using the repository browser.