Ignore:
Timestamp:
2013-12-11T18:34:00+01:00 (7 years ago)
Author:
clem
Message:

remove remaining bugs in LIM3, so that it can run in both regional and global config

Location:
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r4292 r4333  
    104104   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 
    105105 
    106 #if defined key_lim2 
    107    CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim2      ! Choice of boundary condition for sea ice variables  
    108    INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim2_dta  !: = 0 use the initial state as bdy dta ;  
     106#if ( defined key_lim2 || defined key_lim3 ) 
     107   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim       ! Choice of boundary condition for sea ice variables  
     108   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim_dta   !: = 0 use the initial state as bdy dta ;  
    109109                                                              !: = 1 read it in a NetCDF file 
    110110#endif 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4292 r4333  
    197197 
    198198#if defined key_lim2 
    199             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     199            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN  
    200200               ilen1(:) = nblen(:) 
    201201               IF( dta%ll_frld ) THEN 
     
    649649#if defined key_lim2 
    650650            ! sea ice 
    651             IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     651            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
    652652 
    653653               IF( dta%ll_frld ) THEN 
     
    724724               ENDIF 
    725725 
     726            ENDIF 
    726727#endif 
    727728            ! Recalculate field counts 
     
    850851#if defined key_lim2 
    851852         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
    852             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
     853            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
    853854               ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) ) 
    854855               ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) ) 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r4292 r4333  
    4545 
    4646   REAL(wp) ::   epsi20 = 1.e-20_wp  ! module constants 
     47   REAL(wp) ::   epsi10 = 1.e-10_wp  ! min area allowed by ice model 
    4748   !!---------------------------------------------------------------------- 
    4849   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    9192      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !! 
    9293 
     94      INTEGER  ::   jpbound            ! 0 = incoming ice 
     95                                       ! 1 = outgoing ice 
    9396      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices 
    94       INTEGER  ::   ji, jj             ! local scalar 
     97      INTEGER  ::   ji, jj, ii, ij     ! local scalar 
    9598      REAL(wp) ::   zwgt, zwgt1        ! local scalar 
    96       REAL(wp) ::   zinda, ztmelts 
     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 
    97104      !!------------------------------------------------------------------------------ 
    98105      ! 
    99106      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs') 
    100107      ! 
    101       !IF( kt /= nit000 ) THEN 
    102108      jgrd = 1      ! Everything is at T-points here 
    103109      ! 
     
    170176            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth  
    171177            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 
    172194         ENDDO 
    173195         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy ) 
    174196         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
    175197         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) 
    176  
     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 
    177206         DO jb = 1, idx%nblen(jgrd) 
    178207            ji    = idx%nbi(jb,jgrd) 
    179208            jj    = idx%nbj(jb,jgrd) 
    180             ! clem : condition on ice thickness depends on the ice velocity 
    181             ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values  
    182             IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) THEN 
    183                ht_i(ji,jj,jl) = ht_i(ji+1,jj  ,jl) 
    184                a_i (ji,jj,jl) = a_i (ji+1,jj  ,jl) 
    185                ht_s(ji,jj,jl) = ht_s(ji+1,jj  ,jl) 
    186             ENDIF 
    187             IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) THEN 
    188                ht_i(ji,jj,jl) = ht_i(ji-1,jj  ,jl) 
    189                a_i (ji,jj,jl) = a_i (ji-1,jj  ,jl) 
    190                ht_s(ji,jj,jl) = ht_s(ji-1,jj  ,jl) 
    191             ENDIF 
    192             IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) THEN 
    193                ht_i(ji,jj,jl) = ht_i(ji  ,jj+1,jl) 
    194                a_i (ji,jj,jl) = a_i (ji  ,jj+1,jl) 
    195                ht_s(ji,jj,jl) = ht_s(ji  ,jj+1,jl) 
    196             ENDIF 
    197             IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) THEN 
    198                ht_i(ji,jj,jl) = ht_i(ji  ,jj-1,jl) 
    199                a_i (ji,jj,jl) = a_i (ji  ,jj-1,jl) 
    200                ht_s(ji,jj,jl) = ht_s(ji  ,jj-1,jl) 
    201             ENDIF 
    202  
    203             zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - a_i(ji,jj,jl) ) ) ! 0 if no ice 
    204             ! -------------------- 
     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 
    205226            ! Ice and snow volumes 
    206             ! -------------------- 
    207227            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    208228            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) 
    209229 
    210             ! ------------- 
    211             ! Ice salinity 
    212             !--------------- 
    213             sm_i(ji,jj,jl)   = zinda * bulk_sal  + ( 1.0 - zinda ) * s_i_min 
     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 
    214263            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    215              
    216             !---------- 
    217             ! Ice age 
    218             !---------- 
    219             o_i(ji,jj,jl)    = zinda * 1.0   + ( 1.0 - zinda ) 
    220264            oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    221              
    222             !------------------------------ 
    223             ! Sea ice surface temperature 
    224             !------------------------------ 
    225             t_su(ji,jj,jl)   = zinda * 270.0 + ( 1.0 - zinda ) * t_bo(ji,jj) 
    226              
    227             !------------------------------------ 
    228             ! Snow temperature and heat content 
    229             !------------------------------------ 
    230265            DO jk = 1, nlay_s 
    231                t_s(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt 
    232266               ! Snow energy of melting 
    233267               e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
     
    235269               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    236270               ! Multiply by volume, so that heat content in 10^9 Joules 
    237                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    238                     v_s(ji,jj,jl)  / nlay_s 
    239             END DO !jk 
    240              
    241             !----------------------------------------------- 
    242             ! Ice salinities, temperature and heat content  
    243             !----------------------------------------------- 
     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 
    244273            DO jk = 1, nlay_i 
    245                t_i(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt  
    246                s_i(ji,jj,jk,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min 
    247                ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    248                 
     274               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K                   
    249275               ! heat content per unit volume 
    250276               e_i(ji,jj,jk,jl) = zinda * rhoic * & 
    251                     (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    252                     +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
    253                     - rcp      * ( ztmelts - rtt ) ) 
    254                 
     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 ) ) 
    255280               ! Correct dimensions to avoid big values 
    256281               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    257                 
    258282               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    259                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    260                     area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
    261             END DO ! jk 
     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 
    262286 
    263287         END DO !jb 
    264288  
    265  
    266289         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions 
    267290         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
     
    288311#endif 
    289312      !       
    290       !ENDIF !nit000/=0 
    291313      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs') 
    292314      ! 
     
    294316 
    295317 
    296    SUBROUTINE bdy_ice_lim_dyn( kn, pvar1, pvar2, pvar12 ) 
     318   SUBROUTINE bdy_ice_lim_dyn( cd_type ) 
    297319      !!------------------------------------------------------------------------------ 
    298320      !!                 ***  SUBROUTINE bdy_ice_lim_dyn  *** 
    299321      !!                     
    300322      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries. 
    301       !!      kn = 1: u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free 
     323      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free 
    302324      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities 
    303       !!      kn = 2: the stress tensor is set to 0 (i.e. pvar1, pvar2, pvar12) 
    304325      !! 
    305326      !! 2013-06 : C. Rousset 
    306327      !!------------------------------------------------------------------------------ 
    307328      !! 
    308       INTEGER,  INTENT( in    )                           ::   kn     ! set up of ice vel (kn=1) or stress tensor (kn=2) 
    309       REAL(wp), INTENT( inout ), DIMENSION(:,:), OPTIONAL ::   pvar1, pvar2, pvar12  ! stress tensor components (zs1,zs2,zs12)  
     329      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
    310330      INTEGER  ::   jb, jgrd   ! dummy loop indices 
    311331      INTEGER  ::   ji, jj             ! local scalar 
     
    317337      ! 
    318338      DO ib_bdy=1, nb_bdy 
    319       ! 
    320       SELECT CASE( nn_ice_lim(ib_bdy) ) 
    321  
    322       CASE(jp_none) 
    323          CYCLE 
    324  
    325       CASE(jp_frs) 
    326  
    327          IF ( kn == 1 ) THEN ! set up of u_ice and v_ice 
    328             
    329             jgrd = 2      ! u velocity 
    330             DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
    331                ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    332                jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    333                zflag = idx_bdy(ib_bdy)%flagu(jb) 
    334  
    335                IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
    336                   ! one of the two zmsk is always 0 (because of zflag) 
    337                   zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
    338                   zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
     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) 
    339358                   
    340                   ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    341                   u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    342                      &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
    343                      &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    344                ELSE                             ! everywhere else 
    345                   u_ice(ji,jj) = u_oce(ji,jj) 
    346                ENDIF 
    347                ! mask ice velocities 
    348                zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice 
    349                u_ice(ji,jj) = zinda * u_ice(ji,jj) 
    350  
    351             ENDDO 
     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 
    352410             
    353             jgrd = 3      ! v 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)%flagv(jb) 
    358                
    359                IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern 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,jj+1) ) ) ! 0 if no ice 
    362                   zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 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                   v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    366                      &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
    367                      &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    368                ELSE                             ! everywhere else 
    369                   v_ice(ji,jj) = v_oce(ji,jj) 
    370                ENDIF 
    371                ! mask ice velocities 
    372                zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice 
    373                v_ice(ji,jj) = zinda * v_ice(ji,jj) 
    374  
    375             ENDDO 
    376              
    377          CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) 
    378          CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) 
    379   
    380          ELSEIF ( kn == 2 ) THEN ! set up of stress tensor (not sure it works) 
     411         CASE DEFAULT 
     412            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) 
     413         END SELECT 
    381414          
    382             jgrd = 1      ! T grid 
    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  
    387                pvar1 (ji,jj) = 0._wp 
    388                pvar2 (ji,jj) = 0._wp 
    389                pvar12(ji,jj) = 0._wp 
    390  
    391             ENDDO 
    392  
    393          ENDIF 
    394  
    395       CASE DEFAULT 
    396          CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) 
    397       END SELECT 
    398  
    399415      ENDDO 
    400416 
    401417      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn') 
    402   
     418       
    403419    END SUBROUTINE bdy_ice_lim_dyn 
    404420 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4294 r4333  
    488488            DO igrd = 1, jpbgrd 
    489489               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    490                nblendta(igrd,ib_bdy) = kdimsz(1) 
    491                jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     490               !clem nblendta(igrd,ib_bdy) = kdimsz(1) 
     491               !clem jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     492               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 
     493               jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 
    492494            ENDDO 
    493495            CALL iom_close( inum ) 
Note: See TracChangeset for help on using the changeset viewer.