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 trunk/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90 @ 4624

Last change on this file since 4624 was 3614, checked in by acc, 11 years ago

Branch dev_NOC_2012_r3555. #1006. Step 6: Minor code changes and updated namelists to enable successful SETTE testing
File size: 17.9 KB
Line 
1MODULE icbdyn
2
3   !!======================================================================
4   !!                       ***  MODULE  icbdyn  ***
5   !! Iceberg:  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   USE par_oce        ! NEMO parameters
16   USE dom_oce        ! NEMO ocean domain
17   USE phycst         ! NEMO physical constants
18   !
19   USE icb_oce        ! define iceberg arrays
20   USE icbutl         ! iceberg utility routines
21   USE icbdia         ! iceberg budget routines
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   icb_dyn  ! routine called in icbstp.F90 module
27
28   !!----------------------------------------------------------------------
29   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
30   !! $Id:$
31   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   SUBROUTINE icb_dyn()
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE icb_dyn  ***
38      !!
39      !! ** Purpose :   iceberg evolution.
40      !!
41      !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
42      !!----------------------------------------------------------------------
43      REAL(wp)                        ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
44      REAL(wp)                        ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
45      REAL(wp)                        ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
46      REAL(wp)                        ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
47      REAL(wp)                        ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
48      REAL(wp)                        ::   zdt, zdt_2, zdt_6, ze1, ze2
49      LOGICAL                         ::   ll_bounced
50      TYPE(iceberg), POINTER          ::   berg
51      TYPE(point)  , POINTER          ::   pt
52      !!----------------------------------------------------------------------
53
54      ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
55      !                    with I.C.'s:   X=X1 and V=V1
56      !
57      !                                    ; A1=A(X1,V1)
58      !  X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2)
59      !  X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3)
60      !  X4 = X1+  dt*V3 ; V4 = V1+  dt*A3 ; A4=A(X4,V4)
61      !
62      !  Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
63      !  Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
64
65      ! time steps
66      zdt   = berg_dt
67      zdt_2 = zdt * 0.5_wp
68      zdt_6 = zdt / 6._wp
69
70      berg => first_berg                    ! start from the first berg
71      !
72      DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
73         !
74         pt => berg%current_point
75
76         ll_bounced = .false.
77
78
79         ! STEP 1 !
80         ! ====== !
81         zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
82         zyj1 = pt%yj   ;   zvvel1 = pt%vvel
83
84
85         !                                         !**   A1 = A(X1,V1)
86         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
87            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
88         !
89         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
90         zv1 = zvvel1 / ze2
91
92         ! STEP 2 !
93         ! ====== !
94         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
95         ! position using di/dt & djdt   !   V2  in m/s
96         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
97         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
98         !
99         CALL icb_ground( zxi2, zxi1, zu1,   &
100         &                zyj2, zyj1, zv1, ll_bounced )
101
102         !                                         !**   A2 = A(X2,V2)
103         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
104            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
105         !
106         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
107         zv2 = zvvel2 / ze2
108         !
109         ! STEP 3 !
110         ! ====== !
111         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
112         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
113         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
114         !
115         CALL icb_ground( zxi3, zxi1, zu3,   &
116         &                zyj3, zyj1, zv3, ll_bounced )
117
118         !                                         !**   A3 = A(X3,V3)
119         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
120            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
121         !
122         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
123         zv3 = zvvel3 / ze2
124
125         ! STEP 4 !
126         ! ====== !
127         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
128         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
129         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
130
131         CALL icb_ground( zxi4, zxi1, zu4,   &
132         &                zyj4, zyj1, zv4, ll_bounced )
133
134         !                                         !**   A4 = A(X4,V4)
135         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
136            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
137
138         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
139         zv4 = zvvel4 / ze2
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         zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
146         zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
147         zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 )
148         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
149
150         CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
151         &                      zyj_n, zyj1, zvvel_n, ll_bounced )
152
153         pt%uvel = zuvel_n                        !** save in berg structure
154         pt%vvel = zvvel_n
155         pt%xi   = zxi_n
156         pt%yj   = zyj_n
157
158         ! update actual position
159         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
160         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
161
162         berg => berg%next                         ! switch to the next berg
163         !
164      END DO                                  !==  end loop over all bergs  ==!
165      !
166   END SUBROUTINE icb_dyn
167
168
169   SUBROUTINE icb_ground( pi, pi0, pu,   &
170      &                         pj, pj0, pv, ld_bounced )
171      !!----------------------------------------------------------------------
172      !!                  ***  ROUTINE icb_ground  ***
173      !!
174      !! ** Purpose :   iceberg grounding.
175      !!
176      !! ** Method  : - adjust velocity and then put iceberg back to start position
177      !!                NB two possibilities available one of which is hard-coded here
178      !!----------------------------------------------------------------------
179      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position
180      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position
181      REAL(wp), INTENT(inout) ::   pu  , pv     ! current iceberg velocities
182      LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
183      !
184      INTEGER  ::   ii, ii0
185      INTEGER  ::   ij, ij0
186      INTEGER  ::   ibounce_method
187      !!----------------------------------------------------------------------
188      !
189      ld_bounced = .FALSE.
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 = mi1( ii0 )
198      ij0 = mj1( ij0 )
199      ii  = mi1( ii  )
200      ij  = mj1( ij  )
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      ld_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      ibounce_method = 2
216      SELECT CASE ( ibounce_method )
217         CASE ( 1 )
218            pi = pi0
219            pj = pj0
220            pu = 0._wp
221            pv = 0._wp
222         CASE ( 2 )
223            IF( ii0 /= ii ) THEN
224               pi = pi0                   ! return back to the initial position
225               pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
226            ENDIF
227            IF( ij0 /= ij ) THEN
228               pj = pj0                   ! return back to the initial position
229               pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
230            ENDIF
231      END SELECT
232      !
233   END SUBROUTINE icb_ground
234
235
236   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                &
237      &                         pyj, pe2, pvvel, pvvel0, pay, pdt )
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE icb_accel  ***
240      !!
241      !! ** Purpose :   compute the iceberg acceleration.
242      !!
243      !! ** Method  : - sum the terms in the momentum budget
244      !!----------------------------------------------------------------------
245      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg
246      REAL(wp)               , INTENT(in   ) ::   pxi   , pyj      ! berg position in (i,j) referential
247      REAL(wp)               , INTENT(in   ) ::   puvel , pvvel    ! berg velocity [m/s]
248      REAL(wp)               , INTENT(in   ) ::   puvel0, pvvel0   ! initial berg velocity [m/s]
249      REAL(wp)               , INTENT(  out) ::   pe1, pe2         ! horizontal scale factor at (xi,yj)
250      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration
251      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step
252      !
253      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      !
254      REAL(wp), PARAMETER ::   pp_beta      = 1._wp      !
255      REAL(wp), PARAMETER ::   pp_vel_lim   =15._wp      ! max allowed berg speed
256      REAL(wp), PARAMETER ::   pp_accel_lim = 1.e-2_wp   ! max allowed berg acceleration
257      REAL(wp), PARAMETER ::   pp_Cr0       = 0.06_wp    !
258      !
259      INTEGER  ::   itloop
260      REAL(wp) ::   zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi
261      REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
262      REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
263      REAL(wp) ::   z_ocn, z_atm, z_ice
264      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop
265      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi
266      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new
267      !!----------------------------------------------------------------------
268
269      ! Interpolate gridded fields to berg
270      nknberg = berg%number(1)
271      CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     &
272         &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff )
273
274      zM = berg%current_point%mass
275      zT = berg%current_point%thickness               ! total thickness
276      zD = ( rn_rho_bergs / pp_rho_seawater ) * zT    ! draught (keel depth)
277      zF = zT - zD                                    ! freeboard
278      zW = berg%current_point%width
279      zL = berg%current_point%length
280
281      zhi   = MIN( zhi   , zD    )
282      zD_hi = MAX( 0._wp, zD-zhi )
283
284      ! Wave radiation
285      zuwave = zua - zuo   ;   zvwave = zva - zvo     ! Use wind speed rel. to ocean for wave model
286      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
287      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
288      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
289      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
290      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
291      zLcutoff     = 0.125 * zLwavelength
292      zLtop        = 0.25  * zLwavelength
293      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
294      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
295      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
296      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
297      IF( zwmod /= 0._wp ) THEN
298         zuwave = zua/zwmod   ! Wave radiation force acts in wind direction ...       !!gm  this should be the wind rel. to ocean ?
299         zvwave = zva/zwmod
300      ELSE
301         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
302      ENDIF
303
304      ! Weighted drag coefficients
305      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
306      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
307      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
308      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
309
310      zuveln = puvel   ;   zvveln = pvvel ! Copy starting uvel, vvel
311      !
312      DO itloop = 1, 2  ! Iterate on drag coefficients
313         !
314         zus = 0.5 * ( zuveln + puvel )
315         zvs = 0.5 * ( zvveln + pvvel )
316         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
317         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )
318         zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) )
319         !
320         ! Explicit accelerations
321         !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave &
322         !    -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)
323         !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave &
324         !    -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)
325         zaxe = -grav * zssh_x + zwave_rad * zuwave
326         zaye = -grav * zssh_y + zwave_rad * zvwave
327         IF( pp_alpha > 0._wp ) THEN   ! If implicit, use time-level (n) rather than RK4 latest
328            zaxe = zaxe + zff*pvvel0
329            zaye = zaye - zff*puvel0
330         ELSE
331            zaxe = zaxe + zff*pvvel
332            zaye = zaye - zff*puvel
333         ENDIF
334         IF( pp_beta > 0._wp ) THEN    ! If implicit, use time-level (n) rather than RK4 latest
335            zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui)
336            zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi)
337         ELSE
338            zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
339            zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
340         endif
341
342         ! Solve for implicit accelerations
343         IF( pp_alpha + pp_beta > 0._wp ) THEN
344            zlambda = zdrag_ocn + zdrag_atm + zdrag_ice
345            zA11    = 1._wp + pp_beta *pdt*zlambda
346            zA12    =         pp_alpha*pdt*zff
347            zdetA   = 1._wp / ( zA11*zA11 + zA12*zA12 )
348            pax     = zdetA * ( zA11*zaxe + zA12*zaye )
349            pay     = zdetA * ( zA11*zaye - zA12*zaxe )
350         else
351            pax = zaxe   ;   pay = zaye
352         endif
353
354         zuveln = puvel0 + pdt*pax
355         zvveln = pvvel0 + pdt*pay
356         !
357      END DO      ! itloop
358
359      IF( rn_speed_limit > 0._wp ) THEN       ! Limit speed of bergs based on a CFL criteria (if asked)
360         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg
361         IF( zspeed > 0._wp ) THEN
362            zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
363            zspeed_new = zloc_dx / pdt * rn_speed_limit     ! Speed limit as a factor of dx / dt
364            IF( zspeed_new < zspeed ) THEN
365               zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
366               zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction
367               CALL icb_dia_speed()
368            ENDIF
369         ENDIF
370      ENDIF
371      !                                      ! check the speed and acceleration limits
372      IF( ABS( zuveln ) > pp_vel_lim   .OR. ABS( zvveln ) > pp_vel_lim   )   &
373         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity'
374      IF( ABS( pax    ) > pp_accel_lim .OR. ABS( pay    ) > pp_accel_lim )   &
375         WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration'
376      !
377   END SUBROUTINE icb_accel
378
379   !!======================================================================
380END MODULE icbdyn
Note: See TracBrowser for help on using the repository browser.