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.
Changeset 12143 for NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyice.F90 – NEMO

Ignore:
Timestamp:
2019-12-10T12:57:49+01:00 (4 years ago)
Author:
mathiot
Message:

update ENHANCE-02_ISF_nemo to 12072 (sette in progress)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyice.F90

    r11041 r12143  
    5555      INTEGER, INTENT(in) ::   kt   ! Main time step counter 
    5656      ! 
    57       INTEGER ::   jbdy   ! BDY set index 
     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 
    5861      !!---------------------------------------------------------------------- 
    5962      ! controls 
    6063      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing 
    6164      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 
    6266      ! 
    6367      CALL ice_var_glo2eqv 
    6468      ! 
    65       DO jbdy = 1, nb_bdy 
     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 
    6684         ! 
    67          SELECT CASE( cn_ice(jbdy) ) 
    68          CASE('none')   ;   CYCLE 
    69          CASE('frs' )   ;   CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy ) 
    70          CASE DEFAULT 
    71             CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' ) 
    72          END SELECT 
    73          ! 
    74       END DO 
     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 
    75105      ! 
    76106      CALL ice_cor( kt , 0 )      ! -- In case categories are out of bounds, do a remapping 
     
    80110      ! 
    81111      ! controls 
     112      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
    82113      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    83       IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
     114      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    84115      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
    85116      ! 
     
    87118 
    88119 
    89    SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy ) 
     120   SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy, llrim0 ) 
    90121      !!------------------------------------------------------------------------------ 
    91122      !!                 ***  SUBROUTINE bdy_ice_frs  *** 
     
    96127      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382. 
    97128      !!------------------------------------------------------------------------------ 
    98       TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices 
    99       TYPE(OBC_DATA),  INTENT(in) ::   dta     ! OBC external data 
    100       INTEGER,         INTENT(in) ::   kt      ! main time-step counter 
    101       INTEGER,         INTENT(in) ::   jbdy    ! BDY set index 
     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 
    102134      ! 
    103135      INTEGER  ::   jpbound            ! 0 = incoming ice 
    104136      !                                ! 1 = outgoing ice 
     137      INTEGER  ::   ibeg, iend         ! length of rim to be treated (rim 0 or rim 1) 
    105138      INTEGER  ::   i_bdy, jgrd        ! dummy loop indices 
    106139      INTEGER  ::   ji, jj, jk, jl, ib, jb 
    107140      REAL(wp) ::   zwgt, zwgt1        ! local scalar 
    108141      REAL(wp) ::   ztmelts, zdh 
     142      REAL(wp), POINTER  :: flagu, flagv              ! short cuts 
    109143      !!------------------------------------------------------------------------------ 
    110144      ! 
    111145      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 
    112149      ! 
    113150      DO jl = 1, jpl 
    114          DO i_bdy = 1, idx%nblenrim(jgrd) 
     151         DO i_bdy = ibeg, iend 
    115152            ji    = idx%nbi(i_bdy,jgrd) 
    116153            jj    = idx%nbj(i_bdy,jgrd) 
    117154            zwgt  = idx%nbw(i_bdy,jgrd) 
    118155            zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd) 
    119             a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction  
    120             h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth  
    121             h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth 
    122  
     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            ! 
    123174            ! ----------------- 
    124175            ! Pathological case 
     
    135186            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    136187            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )  
    137  
     188            ! 
    138189         ENDDO 
    139190      ENDDO 
    140       CALL lbc_bdy_lnk( 'bdyice', a_i(:,:,:), 'T', 1., jbdy ) 
    141       CALL lbc_bdy_lnk( 'bdyice', h_i(:,:,:), 'T', 1., jbdy ) 
    142       CALL lbc_bdy_lnk( 'bdyice', h_s(:,:,:), 'T', 1., jbdy ) 
    143191 
    144192      DO jl = 1, jpl 
    145          DO i_bdy = 1, idx%nblenrim(jgrd) 
     193         DO i_bdy = ibeg, iend 
    146194            ji = idx%nbi(i_bdy,jgrd) 
    147195            jj = idx%nbj(i_bdy,jgrd) 
    148  
     196            flagu => idx%flagu(i_bdy,jgrd) 
     197            flagv => idx%flagv(i_bdy,jgrd) 
    149198            ! condition on ice thickness depends on the ice velocity 
    150199            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 
    151200            jpbound = 0   ;   ib = ji   ;   jb = jj 
    152201            ! 
    153             IF( u_ice(ji  ,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1 
    154             IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1 
    155             IF( v_ice(ji  ,jj  ) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; jb = jj+1 
    156             IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; jb = jj-1 
     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 
    157218            ! 
    158219            IF( nn_ice_dta(jbdy) == 0 )   jpbound = 0 ; ib = ji ; jb = jj   ! case ice boundaries = initial conditions 
     
    161222            IF( a_i(ib,jb,jl) > 0._wp ) THEN   ! there is ice at the boundary 
    162223               ! 
    163                a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration 
    164                h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice 
    165                h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw 
    166                ! 
    167                SELECT CASE( jpbound ) 
    168                   ! 
    169                CASE( 0 )   ! velocity is inward 
    170                   ! 
    171                   oa_i(ji,jj,  jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age 
    172                   a_ip(ji,jj,  jl) = 0._wp                            ! pond concentration 
    173                   v_ip(ji,jj,  jl) = 0._wp                            ! pond volume 
    174                   t_su(ji,jj,  jl) = rn_ice_tem(jbdy)                 ! temperature surface 
    175                   t_s (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature snw 
    176                   t_i (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature ice 
    177                   s_i (ji,jj,  jl) = rn_ice_sal(jbdy)                 ! salinity 
    178                   sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy)                 ! salinity profile 
    179                   ! 
    180                CASE( 1 )   ! velocity is outward 
    181                   ! 
    182                   oa_i(ji,jj,  jl) = oa_i(ib,jb,  jl) ! age 
    183                   a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) ! pond concentration 
    184                   v_ip(ji,jj,  jl) = v_ip(ib,jb,  jl) ! pond volume 
    185                   t_su(ji,jj,  jl) = t_su(ib,jb,  jl) ! temperature surface 
    186                   t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw 
    187                   t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice 
    188                   s_i (ji,jj,  jl) = s_i (ib,jb,  jl) ! salinity 
    189                   sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile 
    190                   ! 
    191                END SELECT 
     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 
    192242               ! 
    193243               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     
    214264               END DO 
    215265               ! 
     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               ! 
    216274            ELSE   ! no ice at the boundary 
    217275               ! 
     
    225283               t_s (ji,jj,:,jl) = rt0 
    226284               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 
    227290                
    228291               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     
    246309         ! 
    247310      END DO ! jl 
    248  
    249       CALL lbc_bdy_lnk( 'bdyice', a_i (:,:,:)  , 'T', 1., jbdy ) 
    250       CALL lbc_bdy_lnk( 'bdyice', h_i (:,:,:)  , 'T', 1., jbdy ) 
    251       CALL lbc_bdy_lnk( 'bdyice', h_s (:,:,:)  , 'T', 1., jbdy ) 
    252       CALL lbc_bdy_lnk( 'bdyice', oa_i(:,:,:)  , 'T', 1., jbdy ) 
    253       CALL lbc_bdy_lnk( 'bdyice', a_ip(:,:,:)  , 'T', 1., jbdy ) 
    254       CALL lbc_bdy_lnk( 'bdyice', v_ip(:,:,:)  , 'T', 1., jbdy ) 
    255       CALL lbc_bdy_lnk( 'bdyice', s_i (:,:,:)  , 'T', 1., jbdy ) 
    256       CALL lbc_bdy_lnk( 'bdyice', t_su(:,:,:)  , 'T', 1., jbdy ) 
    257       CALL lbc_bdy_lnk( 'bdyice', v_i (:,:,:)  , 'T', 1., jbdy ) 
    258       CALL lbc_bdy_lnk( 'bdyice', v_s (:,:,:)  , 'T', 1., jbdy ) 
    259       CALL lbc_bdy_lnk( 'bdyice', sv_i(:,:,:)  , 'T', 1., jbdy ) 
    260       CALL lbc_bdy_lnk( 'bdyice', t_s (:,:,:,:), 'T', 1., jbdy ) 
    261       CALL lbc_bdy_lnk( 'bdyice', e_s (:,:,:,:), 'T', 1., jbdy ) 
    262       CALL lbc_bdy_lnk( 'bdyice', t_i (:,:,:,:), 'T', 1., jbdy ) 
    263       CALL lbc_bdy_lnk( 'bdyice', e_i (:,:,:,:), 'T', 1., jbdy ) 
    264311      !       
    265312   END SUBROUTINE bdy_ice_frs 
     
    279326      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
    280327      ! 
    281       INTEGER  ::   i_bdy, jgrd      ! dummy loop indices 
    282       INTEGER  ::   ji, jj           ! local scalar 
    283       INTEGER  ::   jbdy             ! BDY set index 
     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) 
    284332      REAL(wp) ::   zmsk1, zmsk2, zflag 
     333      LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
    285334      !!------------------------------------------------------------------------------ 
    286335      IF( ln_timing )   CALL timing_start('bdy_ice_dyn') 
    287336      ! 
    288       DO jbdy=1, nb_bdy 
     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 
    289424         ! 
    290          SELECT CASE( cn_ice(jbdy) ) 
    291          ! 
    292          CASE('none') 
    293             CYCLE 
    294             ! 
    295          CASE('frs') 
    296             ! 
    297             IF( nn_ice_dta(jbdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions  
    298             !                                            !      do not change ice velocity (it is only computed by rheology) 
    299             SELECT CASE ( cd_type ) 
    300             !      
    301             CASE ( 'U' )   
    302                jgrd = 2      ! u velocity 
    303                DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 
    304                   ji    = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 
    305                   jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 
    306                   zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd) 
    307                   ! 
    308                   IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
    309                      ! one of the two zmsk is always 0 (because of zflag) 
    310                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
    311                      zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice 
    312                      !   
    313                      ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice) 
    314                      u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
    315                         &            u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    316                         &            u_ice(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    317                   ELSE                             ! everywhere else 
    318                      u_ice(ji,jj) = 0._wp 
    319                   ENDIF 
    320                   ! 
    321                END DO 
    322                CALL lbc_bdy_lnk( 'bdyice', u_ice(:,:), 'U', -1., jbdy ) 
    323                ! 
    324             CASE ( 'V' ) 
    325                jgrd = 3      ! v velocity 
    326                DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd) 
    327                   ji    = idx_bdy(jbdy)%nbi(i_bdy,jgrd) 
    328                   jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 
    329                   zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd) 
    330                   ! 
    331                   IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries 
    332                      ! one of the two zmsk is always 0 (because of zflag) 
    333                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 
    334                      zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice 
    335                      !   
    336                      ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice) 
    337                      v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
    338                         &            v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    339                         &            v_ice(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    340                   ELSE                             ! everywhere else 
    341                      v_ice(ji,jj) = 0._wp 
    342                   ENDIF 
    343                   ! 
    344                END DO 
    345                CALL lbc_bdy_lnk( 'bdyice', v_ice(:,:), 'V', -1., jbdy ) 
    346                ! 
    347             END SELECT 
    348             ! 
    349          CASE DEFAULT 
    350             CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' ) 
     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 
    351454         END SELECT 
    352          ! 
    353       END DO 
     455      END DO   ! ir 
    354456      ! 
    355457      IF( ln_timing )   CALL timing_stop('bdy_ice_dyn') 
Note: See TracChangeset for help on using the changeset viewer.