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 branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90 @ 3372

Last change on this file since 3372 was 3372, checked in by sga, 12 years ago

NEMO branch dev_r3337_NOCS10_ICB: change all routine names and create more Gurvanistic havoc

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