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.
bdyice_lim.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90 @ 4528

Last change on this file since 4528 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config
File size: 20.7 KB
Line 
1MODULE bdyice_lim
2   !!======================================================================
3   !!                       ***  MODULE  bdyice_lim  ***
4   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (LIM2 and LIM3)
5   !!======================================================================
6   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
7   !!             3.4  !  2011    (D. Storkey)  rewrite in preparation for OBC-BDY merge
8   !!              -   !  2012-01 (C. Rousset)  add lim3 and remove useless jk loop
9   !!----------------------------------------------------------------------
10#if defined   key_bdy   &&  ( defined key_lim2 || defined key_lim3 )
11   !!----------------------------------------------------------------------
12   !!   'key_bdy'            and                 Unstructured Open Boundary Conditions
13   !!   'key_lim2'                                                 LIM-2 sea ice model
14   !!   'key_lim3'                                                 LIM-3 sea ice model
15   !!----------------------------------------------------------------------
16   !!   bdy_ice_lim        : Application of open boundaries to ice
17   !!   bdy_ice_frs        : Application of Flow Relaxation Scheme
18   !!----------------------------------------------------------------------
19   USE timing          ! Timing
20   USE phycst          ! physical constant
21   USE eosbn2          ! equation of state
22   USE oce             ! ocean dynamics and tracers variables
23#if defined key_lim2
24   USE par_ice_2
25   USE ice_2           ! LIM_2 ice variables
26#elif defined key_lim3
27   USE par_ice
28   USE ice             ! LIM_3 ice variables
29#endif
30   USE par_oce         ! ocean parameters
31   USE dom_oce         ! ocean space and time domain variables
32   USE dom_ice          ! sea-ice domain
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE bdy_oce         ! ocean open boundary conditions
35   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
36   USE in_out_manager  ! write to numout file
37   USE lib_mpp         ! distributed memory computing
38   USE lib_fortran     ! to use key_nosignedzero
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   bdy_ice_lim    ! routine called in sbcmod
44   PUBLIC   bdy_ice_lim_dyn ! routine called in limrhg
45
46   REAL(wp) ::   epsi20 = 1.e-20_wp  ! module constants
47   REAL(wp) ::   epsi10 = 1.e-10_wp  ! min area allowed by ice model
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id: bdyice.F90 2715 2011-03-30 15:58:35Z rblod $
51   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE bdy_ice_lim( kt )
56      !!----------------------------------------------------------------------
57      !!                  ***  SUBROUTINE bdy_ice_lim  ***
58      !!
59      !! ** Purpose : - Apply open boundary conditions for ice (LIM2 and LIM3)
60      !!
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT( in ) :: kt     ! Main time step counter
63      !!
64      INTEGER               :: ib_bdy ! Loop index
65      DO ib_bdy=1, nb_bdy
66
67         SELECT CASE( cn_ice_lim(ib_bdy) )
68         CASE('none')
69            CYCLE
70         CASE('frs')
71            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
72         CASE DEFAULT
73            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' )
74         END SELECT
75      ENDDO
76
77   END SUBROUTINE bdy_ice_lim
78
79   SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy )
80      !!------------------------------------------------------------------------------
81      !!                 ***  SUBROUTINE bdy_ice_frs  ***
82      !!                   
83      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case
84      !!              of unstructured open boundaries.
85      !!
86      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
87      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
88      !!------------------------------------------------------------------------------
89      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
90      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
91      INTEGER,         INTENT(in) ::   kt   ! main time-step counter
92      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !!
93
94      INTEGER  ::   jpbound            ! 0 = incoming ice
95                                       ! 1 = outgoing ice
96      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices
97      INTEGER  ::   ji, jj, ii, ij     ! local scalar
98      REAL(wp) ::   zwgt, zwgt1        ! local scalar
99      REAL(wp) ::   zinda, ztmelts, zdh
100
101      REAL(wp), PARAMETER  ::   zsal = 6.3    ! arbitrary salinity    for incoming ice
102      REAL(wp), PARAMETER  ::   ztem = 270.0  ! arbitrary temperature for incoming ice
103      REAL(wp), PARAMETER  ::   zage = 30.0   ! arbitrary age         for incoming ice
104      !!------------------------------------------------------------------------------
105      !
106      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs')
107      !
108      jgrd = 1      ! Everything is at T-points here
109      !
110#if defined key_lim2
111      DO jb = 1, idx%nblen(jgrd)
112         ji    = idx%nbi(jb,jgrd)
113         jj    = idx%nbj(jb,jgrd)
114         zwgt  = idx%nbw(jb,jgrd)
115         zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
116         frld (ji,jj) = ( frld (ji,jj) * zwgt1 + dta%frld (jb) * zwgt ) * tmask(ji,jj,1)     ! Leads fraction
117         hicif(ji,jj) = ( hicif(ji,jj) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ji,jj,1)     ! Ice depth
118         hsnif(ji,jj) = ( hsnif(ji,jj) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ji,jj,1)     ! Snow depth
119
120!         zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - frld(ji,jj) ) ) ! 0 if no ice
121!         !------------------------------
122!         ! Sea ice surface temperature
123!         !------------------------------
124!         sist(ji,jj)   = zinda * 270.0 + ( 1.0 - zinda ) * tfu(ji,jj)
125!         !-----------------------------------------------
126!         ! Ice/snow temperatures and energy stored in brines
127!         !-----------------------------------------------
128!         !!! TO BE CONTIUNED (as LIM3 below) !!!
129!            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
130!
131!               ! Recover in situ values.
132!               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
133!               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
134!               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
135!               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
136!               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
137!               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
138!               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
139!               zindb         = MAX( zindsn, zindic )
140!               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
141!               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
142!               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
143!               hicif(ji,jj)  = zindic * hicif(ji,jj)
144!               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
145!               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
146!               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
147!               zrtt          = 173.15 * rone
148!               ztsn          =          zignm   * tbif(ji,jj,1)  &
149!                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) )
150!               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
151!               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
152!
153!               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
154!               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
155!               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
156!               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
157
158      END DO
159
160      CALL lbc_bdy_lnk( frld,  'T', 1., ib_bdy )                                         ! lateral boundary conditions
161      CALL lbc_bdy_lnk( hicif, 'T', 1., ib_bdy )
162      CALL lbc_bdy_lnk( hsnif, 'T', 1., ib_bdy )
163
164      vt_i(:,:) = hicif(:,:) * frld(:,:)
165      vt_s(:,:) = hsnif(:,:) * frld(:,:)
166      !
167#elif defined key_lim3
168
169      DO jl = 1, jpl
170         DO jb = 1, idx%nblen(jgrd)
171            ji    = idx%nbi(jb,jgrd)
172            jj    = idx%nbj(jb,jgrd)
173            zwgt  = idx%nbw(jb,jgrd)
174            zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
175            a_i (ji,jj,jl) = ( a_i (ji,jj,jl) * zwgt1 + dta%a_i (jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction
176            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth
177            ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
178
179            ! -----------------
180            ! Pathological case
181            ! -----------------
182            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to
183            ! very large transformation from snow to ice (see limthd_dh.F90)
184
185            ! Then, a) transfer the snow excess into the ice (different from limthd_dh)
186            zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
187            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
188            !zdh = MAX( 0._wp, ht_s(ji,jj,jl) * rhosn / rhoic )
189
190            ! recompute ht_i, ht_s
191            ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
192            ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
193
194         ENDDO
195         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
196         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
197         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
198      ENDDO
199      ! retrieve at_i
200      at_i(:,:) = 0._wp
201      DO jl = 1, jpl
202         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
203      END DO
204
205      DO jl = 1, jpl
206         DO jb = 1, idx%nblen(jgrd)
207            ji    = idx%nbi(jb,jgrd)
208            jj    = idx%nbj(jb,jgrd)
209
210            ! condition on ice thickness depends on the ice velocity
211            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
212            jpbound = 0; ii = ji; ij = jj;
213
214            IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj
215            IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj
216            IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1
217            IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1
218
219            zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice
220
221            ! concentration and thickness
222            a_i (ji,jj,jl) = a_i (ii,ij,jl) * zinda
223            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * zinda
224            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * zinda
225
226            ! Ice and snow volumes
227            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
228            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl)
229
230            SELECT CASE( jpbound )
231
232            CASE( 0 ) ! velocity is inward
233
234               ! Ice salinity, age, temperature
235               sm_i(ji,jj,jl)   = zinda * zsal  + ( 1.0 - zinda ) * s_i_min
236               o_i(ji,jj,jl)    = zinda * zage  + ( 1.0 - zinda )
237               t_su(ji,jj,jl)   = zinda * ztem  + ( 1.0 - zinda ) * ztem
238               DO jk = 1, nlay_s
239                  t_s(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt
240               END DO
241               DO jk = 1, nlay_i
242                  t_i(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt 
243                  s_i(ji,jj,jk,jl) = zinda * zsal + ( 1.0 - zinda ) * s_i_min
244               END DO
245               
246            CASE( 1 ) ! velocity is outward
247 
248               ! Ice salinity, age, temperature
249               sm_i(ji,jj,jl)   = zinda * sm_i(ii,ij,jl)  + ( 1.0 - zinda ) * s_i_min
250               o_i(ji,jj,jl)    = zinda * o_i(ii,ij,jl)   + ( 1.0 - zinda )
251               t_su(ji,jj,jl)   = zinda * t_su(ii,ij,jl)  + ( 1.0 - zinda ) * rtt
252               DO jk = 1, nlay_s
253                  t_s(ji,jj,jk,jl) = zinda * t_s(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt
254               END DO
255               DO jk = 1, nlay_i
256                  t_i(ji,jj,jk,jl) = zinda * t_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt
257                  s_i(ji,jj,jk,jl) = zinda * s_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * s_i_min
258               END DO
259
260            END SELECT
261
262            ! contents
263            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
264            oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
265            DO jk = 1, nlay_s
266               ! Snow energy of melting
267               e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
268               ! Change dimensions
269               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
270               ! Multiply by volume, so that heat content in 10^9 Joules
271               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
272            END DO
273            DO jk = 1, nlay_i
274               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K                 
275               ! heat content per unit volume
276               e_i(ji,jj,jk,jl) = zinda * rhoic * &
277                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
278                  +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
279                  - rcp      * ( ztmelts - rtt ) )
280               ! Correct dimensions to avoid big values
281               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
282               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
283               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
284            END DO
285
286
287         END DO !jb
288 
289         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions
290         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
291         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
292         CALL lbc_bdy_lnk(  v_i(:,:,jl), 'T', 1., ib_bdy )
293         CALL lbc_bdy_lnk(  v_s(:,:,jl), 'T', 1., ib_bdy )
294 
295         CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy )
296         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy )
297         CALL lbc_bdy_lnk(  oa_i(:,:,jl), 'T', 1., ib_bdy )
298         CALL lbc_bdy_lnk(   o_i(:,:,jl), 'T', 1., ib_bdy )
299         CALL lbc_bdy_lnk(  t_su(:,:,jl), 'T', 1., ib_bdy )
300         DO jk = 1, nlay_s
301            CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy )
302            CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy )
303         END DO
304         DO jk = 1, nlay_i
305            CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy )
306            CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy )
307         END DO
308
309      END DO !jl
310   
311#endif
312      !     
313      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs')
314      !
315   END SUBROUTINE bdy_ice_frs
316
317
318   SUBROUTINE bdy_ice_lim_dyn( cd_type )
319      !!------------------------------------------------------------------------------
320      !!                 ***  SUBROUTINE bdy_ice_lim_dyn  ***
321      !!                   
322      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries.
323      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
324      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities
325      !!
326      !! 2013-06 : C. Rousset
327      !!------------------------------------------------------------------------------
328      !!
329      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points
330      INTEGER  ::   jb, jgrd   ! dummy loop indices
331      INTEGER  ::   ji, jj             ! local scalar
332      INTEGER  ::   ib_bdy ! Loop index
333      REAL(wp) ::   zmsk1, zmsk2, zflag, zinda 
334     !!------------------------------------------------------------------------------
335      !
336      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_lim_dyn')
337      !
338      DO ib_bdy=1, nb_bdy
339         !
340         SELECT CASE( nn_ice_lim(ib_bdy) )
341
342         CASE('none')
343
344            CYCLE
345           
346         CASE('frs')
347           
348
349            SELECT CASE ( cd_type )
350
351            CASE ( 'U' )
352               
353               jgrd = 2      ! u velocity
354               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
355                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
356                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
357                  zflag = idx_bdy(ib_bdy)%flagu(jb)
358                 
359                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
360                     ! one of the two zmsk is always 0 (because of zflag)
361                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
362                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
363                     
364                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
365                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
366                        &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
367                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
368                  ELSE                             ! everywhere else
369                     !u_ice(ji,jj) = u_oce(ji,jj)
370                     u_ice(ji,jj) = 0._wp
371                  ENDIF
372                  ! mask ice velocities
373                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice
374                  u_ice(ji,jj) = zinda * u_ice(ji,jj)
375                 
376               ENDDO
377
378               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
379               
380            CASE ( 'V' )
381               
382               jgrd = 3      ! v velocity
383               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
384                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
385                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
386                  zflag = idx_bdy(ib_bdy)%flagv(jb)
387                 
388                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
389                     ! one of the two zmsk is always 0 (because of zflag)
390                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
391                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
392                     
393                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
394                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
395                        &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
396                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
397                  ELSE                             ! everywhere else
398                     !v_ice(ji,jj) = v_oce(ji,jj)
399                     v_ice(ji,jj) = 0._wp
400                  ENDIF
401                  ! mask ice velocities
402                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice
403                  v_ice(ji,jj) = zinda * v_ice(ji,jj)
404                 
405               ENDDO
406               
407               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
408               
409            END SELECT
410           
411         CASE DEFAULT
412            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' )
413         END SELECT
414         
415      ENDDO
416
417      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn')
418     
419    END SUBROUTINE bdy_ice_lim_dyn
420
421#else
422   !!---------------------------------------------------------------------------------
423   !!   Default option                                                    Empty module
424   !!---------------------------------------------------------------------------------
425CONTAINS
426   SUBROUTINE bdy_ice_lim( kt )      ! Empty routine
427      WRITE(*,*) 'bdy_ice_lim: You should not have seen this print! error?', kt
428   END SUBROUTINE bdy_ice_lim
429#endif
430
431   !!=================================================================================
432END MODULE bdyice_lim
Note: See TracBrowser for help on using the repository browser.