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

Last change on this file since 13357 was 13357, checked in by mathiot, 2 months ago

ticket #1900: add changes for Merino 2016 works. Results unchanged if flag ln_M2016 is set to .false. . Remaining step: test ln_M2016 = .true.

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