source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90 @ 11518

Last change on this file since 11518 was 11518, checked in by clem, 14 months ago

add the final touch to the famous gaston's branch. More precisely, add the possibility to have melt ponds as input file when using bdy

  • Property svn:keywords set to Id
File size: 25.3 KB
Line 
1MODULE bdyice
2   !!======================================================================
3   !!                       ***  MODULE  bdyice  ***
4   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (SI3)
5   !!======================================================================
6   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
7   !!             3.4  !  2012-01 (C. Rousset)  add new sea ice model
8   !!             4.0  !  2018    (C. Rousset)  SI3 compatibility
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                          SI3 sea ice model
13   !!----------------------------------------------------------------------
14   !!   bdy_ice        : Application of open boundaries to ice
15   !!   bdy_ice_frs    : Application of Flow Relaxation Scheme
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE ice             ! sea-ice: variables
19   USE icevar          ! sea-ice: operations
20   USE icecor          ! sea-ice: corrections
21   USE icectl          ! sea-ice: control prints
22   USE phycst          ! physical constant
23   USE eosbn2          ! equation of state
24   USE par_oce         ! ocean parameters
25   USE dom_oce         ! ocean space and time domain variables
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE bdy_oce         ! ocean open boundary conditions
28   !
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE in_out_manager  ! write to numout file
31   USE lib_mpp         ! distributed memory computing
32   USE lib_fortran     ! to use key_nosignedzero
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   bdy_ice     ! routine called in sbcmod
39   PUBLIC   bdy_ice_dyn ! routine called in icedyn_rhg_evp
40
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE bdy_ice( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  SUBROUTINE bdy_ice  ***
51      !!
52      !! ** Purpose : Apply open boundary conditions for sea ice
53      !!
54      !!----------------------------------------------------------------------
55      INTEGER, INTENT(in) ::   kt   ! Main time step counter
56      !
57      INTEGER ::   jbdy, ir                             ! BDY set index, rim index
58      INTEGER ::   ibeg, iend                           ! length of rim to be treated (rim 0 or rim 1)
59      LOGICAL ::   llrim0                               ! indicate if rim 0 is treated
60      LOGICAL, DIMENSION(4)  :: llsend1, llrecv1        ! indicate how communications are to be carried out
61      !!----------------------------------------------------------------------
62      ! controls
63      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
64      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
65      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
66      !
67      CALL ice_var_glo2eqv
68      !
69      llsend1(:) = .false.   ;   llrecv1(:) = .false.
70      DO ir = 1, 0, -1   ! treat rim 1 before rim 0
71         IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
72         ELSE                 ;   llrim0 = .FALSE.
73         END IF
74         DO jbdy = 1, nb_bdy
75            !
76            SELECT CASE( cn_ice(jbdy) )
77            CASE('none')   ;   CYCLE
78            CASE('frs' )   ;   CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy, llrim0 )
79            CASE DEFAULT
80               CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' )
81            END SELECT
82            !
83         END DO
84         !
85         ! Update bdy points       
86         IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
87         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   END IF
88         DO jbdy = 1, nb_bdy
89            IF( cn_ice(jbdy) == 'frs' ) THEN
90               llsend1(:) = llsend1(:) .OR. lsend_bdyint(jbdy,1,:,ir)   ! possibly every direction, T points
91               llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(jbdy,1,:,ir)   ! possibly every direction, T points
92            END IF
93         END DO   ! jbdy
94         IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
95            ! exchange 3d arrays
96            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
97                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
98                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                &
99                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      )
100            ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk
101            CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
102            CALL lbc_lnk_multi( 'bdyice', t_i , 'T', 1., e_i , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
103         END IF
104      END DO   ! ir
105      !
106      CALL ice_cor( kt , 0 )      ! -- In case categories are out of bounds, do a remapping
107      !                           !    i.e. inputs have not the same ice thickness distribution (set by rn_himean)
108      !                           !         than the regional simulation
109      CALL ice_var_agg(1)
110      !
111      ! controls
112      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
113      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
114      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
115      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
116      !
117   END SUBROUTINE bdy_ice
118
119
120   SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy, llrim0 )
121      !!------------------------------------------------------------------------------
122      !!                 ***  SUBROUTINE bdy_ice_frs  ***
123      !!                   
124      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields
125      !!
126      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
127      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
128      !!------------------------------------------------------------------------------
129      TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices
130      TYPE(OBC_DATA),  INTENT(in) ::   dta      ! OBC external data
131      INTEGER,         INTENT(in) ::   kt       ! main time-step counter
132      INTEGER,         INTENT(in) ::   jbdy     ! BDY set index
133      LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
134      !
135      INTEGER  ::   jpbound            ! 0 = incoming ice
136      !                                ! 1 = outgoing ice
137      INTEGER  ::   ibeg, iend         ! length of rim to be treated (rim 0 or rim 1)
138      INTEGER  ::   i_bdy, jgrd        ! dummy loop indices
139      INTEGER  ::   ji, jj, jk, jl, ib, jb
140      REAL(wp) ::   zwgt, zwgt1        ! local scalar
141      REAL(wp) ::   ztmelts, zdh
142      REAL(wp), POINTER  :: flagu, flagv              ! short cuts
143      !!------------------------------------------------------------------------------
144      !
145      jgrd = 1      ! Everything is at T-points here
146      IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(jgrd)
147      ELSE                ;   ibeg = idx%nblenrim0(jgrd)+1   ;   iend = idx%nblenrim(jgrd)
148      END IF
149      !
150      DO jl = 1, jpl
151         DO i_bdy = ibeg, iend
152            ji    = idx%nbi(i_bdy,jgrd)
153            jj    = idx%nbj(i_bdy,jgrd)
154            zwgt  = idx%nbw(i_bdy,jgrd)
155            zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd)
156            a_i (ji,jj,  jl) = ( a_i (ji,jj,  jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  concentration
157            h_i (ji,jj,  jl) = ( h_i (ji,jj,  jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  depth
158            h_s (ji,jj,  jl) = ( h_s (ji,jj,  jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
159            t_i (ji,jj,:,jl) = ( t_i (ji,jj,:,jl) * zwgt1 + dta%t_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  temperature
160            t_s (ji,jj,:,jl) = ( t_s (ji,jj,:,jl) * zwgt1 + dta%t_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow temperature
161            t_su(ji,jj,  jl) = ( t_su(ji,jj,  jl) * zwgt1 + dta%tsu(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Surf temperature
162            s_i (ji,jj,  jl) = ( s_i (ji,jj,  jl) * zwgt1 + dta%s_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  salinity
163            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration
164            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth
165            !
166            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
167            !
168            ! make sure ponds = 0 if no ponds scheme
169            IF( .NOT.ln_pnd ) THEN
170               a_ip(ji,jj,jl) = 0._wp
171               h_ip(ji,jj,jl) = 0._wp
172            ENDIF
173            !
174            ! -----------------
175            ! Pathological case
176            ! -----------------
177            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to
178            ! very large transformation from snow to ice (see icethd_dh.F90)
179
180            ! Then, a) transfer the snow excess into the ice (different from icethd_dh)
181            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
182            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
183            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
184
185            ! recompute h_i, h_s
186            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
187            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos ) 
188            !
189         ENDDO
190      ENDDO
191
192      DO jl = 1, jpl
193         DO i_bdy = ibeg, iend
194            ji = idx%nbi(i_bdy,jgrd)
195            jj = idx%nbj(i_bdy,jgrd)
196            flagu => idx%flagu(i_bdy,jgrd)
197            flagv => idx%flagv(i_bdy,jgrd)
198            ! condition on ice thickness depends on the ice velocity
199            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
200            jpbound = 0   ;   ib = ji   ;   jb = jj
201            !
202            IF( flagu ==  1. )   THEN
203               IF( ji+1 > jpi )   CYCLE
204               IF( u_ice(ji  ,jj  ) < 0. )   jpbound = 1 ; ib = ji+1
205            END IF
206            IF( flagu == -1. )   THEN
207               IF( ji-1 < 1   )   CYCLE
208               IF( u_ice(ji-1,jj  ) < 0. )   jpbound = 1 ; ib = ji-1
209            END IF
210            IF( flagv ==  1. )   THEN
211               IF( jj+1 > jpj )   CYCLE
212               IF( v_ice(ji  ,jj  ) < 0. )   jpbound = 1 ; jb = jj+1
213            END IF
214            IF( flagv == -1. )   THEN
215               IF( jj-1 < 1   )   CYCLE
216               IF( v_ice(ji  ,jj-1) < 0. )   jpbound = 1 ; jb = jj-1
217            END IF
218            !
219            IF( nn_ice_dta(jbdy) == 0 )   jpbound = 0 ; ib = ji ; jb = jj   ! case ice boundaries = initial conditions
220            !                                                               !      do not make state variables dependent on velocity
221            !
222            IF( a_i(ib,jb,jl) > 0._wp ) THEN   ! there is ice at the boundary
223               !
224               a_i (ji,jj,  jl) = a_i (ib,jb,  jl)
225               h_i (ji,jj,  jl) = h_i (ib,jb,  jl)
226               h_s (ji,jj,  jl) = h_s (ib,jb,  jl)
227               t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl)
228               t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl)
229               t_su(ji,jj,  jl) = t_su(ib,jb,  jl)
230               s_i (ji,jj,  jl) = s_i (ib,jb,  jl)
231               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
232               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
233               !
234               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
235               !
236               ! ice age
237               IF    ( jpbound == 0 ) THEN  ! velocity is inward
238                  oa_i(ji,jj,jl) = rice_age(jbdy) * a_i(ji,jj,jl)
239               ELSEIF( jpbound == 1 ) THEN  ! velocity is outward
240                  oa_i(ji,jj,jl) = oa_i(ib,jb,jl)
241               ENDIF
242               !
243               IF( nn_icesal == 1 ) THEN     ! if constant salinity
244                  s_i (ji,jj  ,jl) = rn_icesal
245                  sz_i(ji,jj,:,jl) = rn_icesal
246               ENDIF
247               !
248               ! global fields
249               v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)                       ! volume ice
250               v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)                       ! volume snw
251               sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
252               DO jk = 1, nlay_s
253                  e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )   ! enthalpy in J/m3
254                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s           ! enthalpy in J/m2
255               END DO               
256               DO jk = 1, nlay_i
257                  ztmelts          = - rTmlt  * sz_i(ji,jj,jk,jl)             ! Melting temperature in C
258                  t_i(ji,jj,jk,jl) = MIN( t_i(ji,jj,jk,jl), ztmelts + rt0 )   ! Force t_i to be lower than melting point => likely conservation issue
259                  !
260                  e_i(ji,jj,jk,jl) = rhoi * ( rcpi  * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) )           &   ! enthalpy in J/m3
261                     &                      + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) )   &
262                     &                      - rcp   *   ztmelts )                 
263                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i                            ! enthalpy in J/m2
264               END DO
265               !
266               ! melt ponds
267               IF( a_i(ji,jj,jl) > epsi10 ) THEN
268                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl)
269               ELSE
270                  a_ip_frac(ji,jj,jl) = 0._wp
271               ENDIF
272               v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
273               !
274            ELSE   ! no ice at the boundary
275               !
276               a_i (ji,jj,  jl) = 0._wp
277               h_i (ji,jj,  jl) = 0._wp
278               h_s (ji,jj,  jl) = 0._wp
279               oa_i(ji,jj,  jl) = 0._wp
280               a_ip(ji,jj,  jl) = 0._wp
281               v_ip(ji,jj,  jl) = 0._wp
282               t_su(ji,jj,  jl) = rt0
283               t_s (ji,jj,:,jl) = rt0
284               t_i (ji,jj,:,jl) = rt0 
285
286               a_ip_frac(ji,jj,jl) = 0._wp
287               h_ip     (ji,jj,jl) = 0._wp
288               a_ip     (ji,jj,jl) = 0._wp
289               v_ip     (ji,jj,jl) = 0._wp
290               
291               IF( nn_icesal == 1 ) THEN     ! if constant salinity
292                  s_i (ji,jj  ,jl) = rn_icesal
293                  sz_i(ji,jj,:,jl) = rn_icesal
294               ELSE                          ! if variable salinity
295                  s_i (ji,jj,jl)   = rn_simin
296                  sz_i(ji,jj,:,jl) = rn_simin
297               ENDIF
298               !
299               ! global fields
300               v_i (ji,jj,  jl) = 0._wp
301               v_s (ji,jj,  jl) = 0._wp
302               sv_i(ji,jj,  jl) = 0._wp
303               e_s (ji,jj,:,jl) = 0._wp
304               e_i (ji,jj,:,jl) = 0._wp
305
306            ENDIF
307                       
308         END DO
309         !
310      END DO ! jl
311      !     
312   END SUBROUTINE bdy_ice_frs
313
314
315   SUBROUTINE bdy_ice_dyn( cd_type )
316      !!------------------------------------------------------------------------------
317      !!                 ***  SUBROUTINE bdy_ice_dyn  ***
318      !!                   
319      !! ** Purpose : Apply dynamics boundary conditions for sea-ice.
320      !!
321      !! ** Method :  if this adjacent grid point is not ice free, then u_ice and v_ice take its value
322      !!              if                          is     ice free, then u_ice and v_ice are unchanged by BDY
323      !!                                                           they keep values calculated in rheology
324      !!
325      !!------------------------------------------------------------------------------
326      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points
327      !
328      INTEGER  ::   i_bdy, jgrd       ! dummy loop indices
329      INTEGER  ::   ji, jj            ! local scalar
330      INTEGER  ::   jbdy, ir     ! BDY set index, rim index
331      INTEGER  ::   ibeg, iend   ! length of rim to be treated (rim 0 or rim 1)
332      REAL(wp) ::   zmsk1, zmsk2, zflag
333      LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out
334      !!------------------------------------------------------------------------------
335      IF( ln_timing )   CALL timing_start('bdy_ice_dyn')
336      !
337      llsend2(:) = .false.   ;   llrecv2(:) = .false.
338      llsend3(:) = .false.   ;   llrecv3(:) = .false.
339      DO ir = 1, 0, -1
340         DO jbdy = 1, nb_bdy
341            !
342            SELECT CASE( cn_ice(jbdy) )
343               !
344            CASE('none')
345               CYCLE
346               !
347            CASE('frs')
348               !
349               IF( nn_ice_dta(jbdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions
350               !                                            !      do not change ice velocity (it is only computed by rheology)
351               SELECT CASE ( cd_type )
352                  !     
353               CASE ( 'U' ) 
354                  jgrd = 2      ! u velocity
355                  IF( ir == 0 ) THEN   ;   ibeg = 1                                 ;   iend = idx_bdy(jbdy)%nblenrim0(jgrd)
356                  ELSE                 ;   ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1   ;   iend = idx_bdy(jbdy)%nblenrim(jgrd)
357                  END IF
358                  DO i_bdy = ibeg, iend
359                     ji    = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
360                     jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
361                     zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd)
362                     !     i-1  i   i    |  !        i  i i+1 |  !          i  i i+1 |
363                     !      >  ice  >    |  !        o  > ice |  !          o  >  o  |     
364                     ! => set at u_ice(i-1) !  => set to O       !  => unchanged
365                     IF( zflag == -1. .AND. ji > 1 .AND. ji < jpi )   THEN 
366                        IF    ( vt_i(ji  ,jj) > 0. )   THEN   ;   u_ice(ji,jj) = u_ice(ji-1,jj) 
367                        ELSEIF( vt_i(ji+1,jj) > 0. )   THEN   ;   u_ice(ji,jj) = 0._wp
368                        END IF
369                     END IF
370                     ! |    i  i+1 i+1        !  |  i   i i+1        !  | i  i i+1
371                     ! |    >  ice  >         !  | ice  >  o         !  | o  >  o   
372                     ! => set at u_ice(i+1)   !     => set to O      !     =>  unchanged
373                     IF( zflag ==  1. .AND. ji+1 < jpi+1 )   THEN
374                        IF    ( vt_i(ji+1,jj) > 0. )   THEN   ;   u_ice(ji,jj) = u_ice(ji+1,jj)
375                        ELSEIF( vt_i(ji  ,jj) > 0. )   THEN   ;   u_ice(ji,jj) = 0._wp
376                        END IF
377                     END IF
378                     !
379                     IF( zflag ==  0. )   u_ice(ji,jj) = 0._wp   ! u_ice = 0  if north/south bdy 
380                     !
381                  END DO
382                  !
383               CASE ( 'V' )
384                  jgrd = 3      ! v velocity
385                  IF( ir == 0 ) THEN   ;   ibeg = 1                                 ;   iend = idx_bdy(jbdy)%nblenrim0(jgrd)
386                  ELSE                 ;   ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1   ;   iend = idx_bdy(jbdy)%nblenrim(jgrd)
387                  END IF
388                  DO i_bdy = ibeg, iend
389                     ji    = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
390                     jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
391                     zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd)
392                     !                         !      ice   (jj+1)       !       o    (jj+1)
393                     !       ^    (jj  )       !       ^    (jj  )       !       ^    (jj  )       
394                     !      ice   (jj  )       !       o    (jj  )       !       o    (jj  )       
395                     !       ^    (jj-1)       !                         !
396                     ! => set to u_ice(jj-1)   !  =>   set to 0          !   => unchanged       
397                     IF( zflag == -1. .AND. jj > 1 .AND. jj < jpj )   THEN                 
398                        IF    ( vt_i(ji,jj  ) > 0. )   THEN   ;   v_ice(ji,jj) = v_ice(ji,jj-1)
399                        ELSEIF( vt_i(ji,jj+1) > 0. )   THEN   ;   v_ice(ji,jj) = 0._wp
400                        END IF
401                     END IF
402                     !       ^    (jj+1)       !                         !             
403                     !      ice   (jj+1)       !       o    (jj+1)       !       o    (jj+1)       
404                     !       ^    (jj  )       !       ^    (jj  )       !       ^    (jj  )
405                     !   ________________      !  ____ice___(jj  )_      !  _____o____(jj  )
406                     ! => set to u_ice(jj+1)   !    => set to 0          !    => unchanged 
407                     IF( zflag ==  1. .AND. jj < jpj )   THEN             
408                        IF    ( vt_i(ji,jj+1) > 0. )   THEN   ;   v_ice(ji,jj) = v_ice(ji,jj+1)
409                        ELSEIF( vt_i(ji,jj  ) > 0. )   THEN   ;   v_ice(ji,jj) = 0._wp
410                        END IF
411                     END IF
412                     !
413                     IF( zflag ==  0. )   v_ice(ji,jj) = 0._wp   ! v_ice = 0  if west/east bdy 
414                     !
415                  END DO
416                  !
417               END SELECT
418               !
419            CASE DEFAULT
420               CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' )
421            END SELECT
422            !
423         END DO    ! jbdy
424         !
425         SELECT CASE ( cd_type )       
426         CASE ( 'U' ) 
427         IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
428         IF( nn_hls == 1 ) THEN   ;   llsend2(:) = .false.   ;   llrecv2(:) = .false.   ;   END IF
429            DO jbdy = 1, nb_bdy
430               IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN
431                  llsend2(:) = llsend2(:) .OR. lsend_bdyint(jbdy,2,:,ir)   ! possibly every direction, U points
432                  llsend2(1) = llsend2(1) .OR. lsend_bdyext(jbdy,2,1,ir)   ! neighbour might search point towards its west bdy
433                  llrecv2(:) = llrecv2(:) .OR. lrecv_bdyint(jbdy,2,:,ir)   ! possibly every direction, U points
434                  llrecv2(2) = llrecv2(2) .OR. lrecv_bdyext(jbdy,2,2,ir)   ! might search point towards east bdy
435               END IF
436            END DO
437            IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
438               CALL lbc_lnk( 'bdyice', u_ice, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
439            END IF
440         CASE ( 'V' )
441         IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
442         IF( nn_hls == 1 ) THEN   ;   llsend3(:) = .false.   ;   llrecv3(:) = .false.   ;   END IF
443            DO jbdy = 1, nb_bdy
444               IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN
445                  llsend3(:) = llsend3(:) .OR. lsend_bdyint(jbdy,3,:,ir)   ! possibly every direction, V points
446                  llsend3(3) = llsend3(3) .OR. lsend_bdyext(jbdy,3,3,ir)   ! neighbour might search point towards its south bdy
447                  llrecv3(:) = llrecv3(:) .OR. lrecv_bdyint(jbdy,3,:,ir)   ! possibly every direction, V points
448                  llrecv3(4) = llrecv3(4) .OR. lrecv_bdyext(jbdy,3,4,ir)   ! might search point towards north bdy
449               END IF
450            END DO
451            IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
452               CALL lbc_lnk( 'bdyice', v_ice, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
453            END IF
454         END SELECT
455      END DO   ! ir
456      !
457      IF( ln_timing )   CALL timing_stop('bdy_ice_dyn')
458      !
459    END SUBROUTINE bdy_ice_dyn
460
461#else
462   !!---------------------------------------------------------------------------------
463   !!   Default option                                                    Empty module
464   !!---------------------------------------------------------------------------------
465CONTAINS
466   SUBROUTINE bdy_ice( kt )      ! Empty routine
467      IMPLICIT NONE
468      INTEGER, INTENT( in ) :: kt
469      WRITE(*,*) 'bdy_ice: You should not have seen this print! error?', kt
470   END SUBROUTINE bdy_ice
471#endif
472
473   !!=================================================================================
474END MODULE bdyice
Note: See TracBrowser for help on using the repository browser.