source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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