Changeset 6763


Ignore:
Timestamp:
2016-06-30T17:17:35+02:00 (4 years ago)
Author:
clem
Message:

new parameterization from Lupkes et al. (2012) of ice-atm and ocean-atm drags (dependent on ice morphology) ⇒ ln_Cd_L12 in namelist_ref

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r6747 r6763  
    9090   rn_gamma       =    0.2          !     (ln_landfast=T)  fraction of ocean depth that ice must reach to initiate landfast 
    9191                                    !                      recommended range: [0.1 ; 0.25] 
    92    rn_icebfr      =    1.e-05       !     (ln_landfast=T)  maximum bottom stress per unit area of contact                   
     92   rn_icebfr      =    5.6e-06      !     (ln_landfast=T)  maximum bottom stress per unit area of contact                   
     93                                    !                      a very large value ensures ice velocity=0 even with a small contact area 
    9394                                    !                      recommended range: ?? 
    9495/ 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6319 r6763  
    336336   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
    337337                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 
     338   ln_Cd_L12   = .false.   !  Modify the drag ice-atm and oce-atm depending on ice concentration 
     339                           !  This parameterization is from Lupkes et al. (JGR 2012) 
    338340/ 
    339341!----------------------------------------------------------------------- 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6746 r6763  
    422422   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                      !: ice temperatures 
    423423   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b           !: ice velocity 
     424   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                     !: ice concentration (total) 
    424425             
    425426   !!-------------------------------------------------------------------------- 
     
    528529         &      oa_i_b (jpi,jpj,jpl)                                                    , STAT=ierr(ii) ) 
    529530      ii = ii + 1 
    530       ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 
     531      ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) ) 
    531532       
    532533      ! * Ice thickness distribution variables 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6746 r6763  
    116116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2  ! alpha and beta from Bouillon 2009 and 2013 
    117117      REAL(wp) ::   zm1, zm2, zm3                              ! ice/snow mass 
    118       REAL(wp) ::   delta, zp_delf, zds2 
     118      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 
    119119      REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel  ! temporary scalars 
    120120 
     
    134134      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2               ! internal stresses 
    135135       
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt, zds                  ! tension and shear 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                       ! shear 
    137137      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12            ! stress tensor components 
    138138      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr     ! check convergence 
     
    149149      CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    150150      CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
    151       CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     151      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    152152      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 ) 
    153153 
     
    257257         END DO 
    258258      END DO 
    259  
    260259      ! 
    261260      !------------------------------------------------------------------------------! 
     
    266265      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    267266         !                                            !----------------------!         
    268          ! Convergence test 
    269          IF(ln_ctl) THEN 
     267         IF(ln_ctl) THEN   ! Convergence test 
    270268            DO jj = 1, jpjm1 
    271269               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     
    275273 
    276274         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    277          DO jj = 2, jpjm1 
    278             DO ji = fs_2, fs_jpim1 
    279                 
    280                ! divergence at T points 
    281                divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    282                   &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    283                   &            ) * r1_e12t(ji,jj) 
    284  
    285                ! tension at T points 
    286                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    287                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    288                   &         ) * r1_e12t(ji,jj) 
     275         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     276            DO ji = 1, fs_jpim1 
    289277 
    290278               ! shear at F points 
     
    292280                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    293281                  &         ) * r1_e12f(ji,jj) 
     282 
     283            END DO 
     284         END DO 
     285         CALL lbc_lnk( zds, 'F', 1. ) 
     286 
     287         DO jj = 2, jpjm1 
     288            DO ji = 2, jpim1 ! no vector loop 
     289 
     290               ! shear**2 at T points (doc eq. A16) 
     291               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     292                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     293                  &   ) * 0.25_wp * r1_e12t(ji,jj) 
     294               
     295               ! divergence at T points 
     296               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     297                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     298                  &    ) * r1_e12t(ji,jj) 
     299               zdiv2 = zdiv * zdiv 
     300                
     301               ! tension at T points 
     302               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     303                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     304                  &   ) * r1_e12t(ji,jj) 
     305               zdt2 = zdt * zdt 
     306                
     307               ! delta at T points 
     308               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )   
     309 
     310               ! P/delta at T points 
     311               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 
     312                
     313               ! stress at T points 
     314               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
     315               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     316              
     317            END DO 
     318         END DO 
     319         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     320 
     321         DO jj = 1, jpjm1 
     322            DO ji = 1, jpim1 
     323 
     324               ! P/delta at F points 
     325               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     326                
     327               ! stress at F points 
     328               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     329 
     330            END DO 
     331         END DO 
     332         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     333  
     334 
     335         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
     336         DO jj = 2, jpjm1 
     337            DO ji = fs_2, fs_jpim1                
     338 
     339               ! U points 
     340               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     341                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     342                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     343                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     344                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     345                  &                  ) * r1_e12u(ji,jj) 
     346 
     347               ! V points 
     348               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     349                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     350                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     351                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     352                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     353                  &                  ) * r1_e12v(ji,jj) 
    294354 
    295355               ! u_ice at V point 
     
    303363            END DO 
    304364         END DO 
    305          CALL lbc_lnk( zds, 'F', 1. ) 
    306  
    307          DO jj = 2, jpjm1 
    308             DO ji = 2, jpim1 ! no vector loop 
    309  
    310                ! shear**2 at T points (doc eq. A16) 
    311                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
    312                   &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
    313                   &   ) * 0.25_wp * r1_e12t(ji,jj) 
    314                
    315                ! delta at T points 
    316                delta          = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )   
    317                delta_i(ji,jj) = delta + rn_creepl 
    318  
    319                ! P/delta at T points 
    320                zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj) 
    321             END DO 
    322          END DO 
    323          CALL lbc_lnk( zp_delt, 'T', 1. ) 
    324  
    325          ! --- Stress tensor --- ! 
    326          DO jj = 2, jpjm1 
    327             DO ji = 2, jpim1 ! no vector loop 
    328                 
    329                ! P/delta at F points 
    330                zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
    331                 
    332                ! stress tensor at T points 
    333                zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1 
    334                zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2  )                      ) * z1_alph2 
    335                ! F points 
    336                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf        * ( zds(ji,jj) * z1_ecc2  ) * 0.5_wp             ) * z1_alph2 
    337             END DO 
    338          END DO 
    339          CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    340   
    341          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    342          DO jj = 2, jpjm1 
    343             DO ji = fs_2, fs_jpim1                
    344                ! U points 
    345                zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    346                   &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
    347                   &                    ) * r1_e2u(ji,jj)                                                                      & 
    348                   &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
    349                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    350                   &                  ) * r1_e12u(ji,jj) 
    351  
    352                ! V points 
    353                zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    354                   &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    355                   &                    ) * r1_e1v(ji,jj)                                                                      & 
    356                   &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
    357                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    358                   &                  ) * r1_e12v(ji,jj) 
    359             END DO 
    360          END DO 
    361365         ! 
    362366         ! --- Computation of ice velocity --- ! 
     
    368372               DO ji = fs_2, fs_jpim1 
    369373 
    370                   zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    371                      &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
    372  
     374                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
     375                   
    373376                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    374                   zTauO   =   za2(ji,jj) * rhoco * zvel 
    375                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     377                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     378                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     379                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    376380                  zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj) 
    377381                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     
    403407                     &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)  
    404408                   
    405                   zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    406                      &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
    407  
     409                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     410                   
    408411                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    409                   zTauO   =   za1(ji,jj) * rhoco * zvel 
    410                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     412                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     413                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     414                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    411415                  zCor    =   zcor1(ji,jj)  * zv_ice1 
    412416                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     
    435439               DO ji = fs_2, fs_jpim1 
    436440 
    437                   zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    438                      &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     441                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    439442 
    440443                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    441                   zTauO   =   za1(ji,jj) * rhoco * zvel 
    442                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     444                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     445                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     446                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    443447                  zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj) 
    444448                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     
    469473                     &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    470474 
    471                   zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    472                      &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     475                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
    473476          
    474477                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    475                   zTauO   =   za2(ji,jj) * rhoco * zvel 
    476                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     478                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     479                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     480                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    477481                  zCor    =   zcor2(ji,jj)  * zu_ice2 
    478482                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     
    485489                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    486490                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
    487  
    488491               END DO 
    489492            END DO 
     
    499502         ENDIF 
    500503          
    501          ! Convergence test 
    502          IF(ln_ctl) THEN 
     504         IF(ln_ctl) THEN   ! Convergence test 
    503505            DO jj = 2 , jpjm1 
    504506               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     
    539541      ! 5) Recompute delta, shear and div (inputs for mechanical redistribution)  
    540542      !------------------------------------------------------------------------------! 
    541  
    542       ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     543      DO jj = 1, jpjm1 
     544         DO ji = 1, fs_jpim1 
     545 
     546            ! shear at F points 
     547            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     548               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     549               &         ) * r1_e12f(ji,jj) 
     550 
     551         END DO 
     552      END DO            
     553      CALL lbc_lnk( zds, 'F', 1. ) 
     554       
    543555      DO jj = 2, jpjm1 
    544          DO ji = fs_2, fs_jpim1 
    545              
     556         DO ji = 2, jpim1 ! no vector loop 
     557             
     558            ! tension**2 at T points 
     559            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     560               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     561               &   ) * r1_e12t(ji,jj) 
     562            zdt2 = zdt * zdt 
     563             
     564            ! shear**2 at T points (doc eq. A16) 
     565            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     566               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     567               &   ) * 0.25_wp * r1_e12t(ji,jj) 
     568             
     569            ! shear at T points 
     570            shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     571 
    546572            ! divergence at T points 
    547573            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     
    549575               &            ) * r1_e12t(ji,jj) 
    550576             
    551             ! tension at T points 
    552             zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    553                &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    554                &         ) * r1_e12t(ji,jj) 
    555              
    556             ! shear at F points 
    557             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    558                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    559                &         ) * r1_e12f(ji,jj) 
    560              
     577            ! delta at T points 
     578            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
     579            delta_i(ji,jj) = zdelta + rn_creepl 
     580 
    561581         END DO 
    562582      END DO 
    563       CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. ) 
    564        
    565  
    566       DO jj = 2, jpjm1 
    567          DO ji = 2, jpim1 ! no vector loop 
    568              
    569             ! shear**2 at T points (doc eq. A16) 
    570             zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
    571                &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
    572                &   ) * 0.25_wp * r1_e12t(ji,jj) 
    573              
    574             ! delta at T points 
    575             delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 )   
    576             delta_i(ji,jj) = delta + rn_creepl 
    577              
    578             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 ) 
    579          END DO 
    580       END DO 
    581       CALL lbc_lnk_multi( delta_i, 'T', 1.,  shear_i, 'T', 1. ) 
     583      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
    582584       
    583585      ! --- Store the stress tensor for the next time step --- ! 
     
    626628      CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    627629      CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
    628       CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     630      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    629631      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 ) 
    630632 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r6515 r6763  
    137137 
    138138      zalb(:,:) = 0._wp 
    139       WHERE     ( SUM( a_i_b, dim=3 ) <= epsi06 )  ;  zalb(:,:) = 0.066_wp 
    140       ELSEWHERE                                    ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     139      WHERE     ( at_i_b <= epsi06 )  ;  zalb(:,:) = 0.066_wp 
     140      ELSEWHERE                       ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
    141141      END WHERE 
    142142      IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output 
    143143 
    144       zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) )       
     144      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - at_i_b )       
    145145      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output 
    146146 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6746 r6763  
    185185      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    186186      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    187       CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base  
     187      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base  
    188188      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    189189      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6399 r6763  
    1616   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE 
    1717   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk 
     18   !!            4.0  !  2016-06  (C. Rousset) Add new param of drags with sea-ice (Lupkes at al 2012) 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    4546   USE lib_fortran     ! to use key_nosignedzero 
    4647#if defined key_lim3 
    47    USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b 
     48   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 
    4849   USE limthd_dh       ! for CALL lim_thd_snwblow 
    4950#elif defined key_lim2 
     
    6061   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module 
    6162#endif 
    62    PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
     63   PUBLIC   turb_core_2z         ! routine called in sbcblk_mfs module 
    6364 
    6465   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
     
    7677 
    7778   !                                             !!! CORE bulk parameters 
    78    REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density 
    79    REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air 
    80    REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization 
    81    REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation 
    82    REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant 
    83    REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
    84    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant 
     79   REAL(wp), PARAMETER ::   rhoa   =    1.22        ! air density 
     80   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air 
     81   REAL(wp), PARAMETER ::   Lv     =    2.5e6       ! latent heat of vaporization 
     82   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
     83   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
     84   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
     85   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    8586 
    8687   !                                  !!* Namelist namsbc_core : CORE bulk parameters 
     
    9192   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9293   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
    93  
     94   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     95 
     96   ! 
     97   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     98    
    9499   !! * Substitutions 
    95100#  include "domzgr_substitute.h90" 
     
    102107CONTAINS 
    103108 
     109   INTEGER FUNCTION sbc_blk_core_alloc() 
     110      !!------------------------------------------------------------------- 
     111      !!             ***  ROUTINE sbc_blk_core_alloc (clem) *** 
     112      !!------------------------------------------------------------------- 
     113      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_core_alloc ) 
     114      ! 
     115      IF( lk_mpp                  )   CALL mpp_sum( sbc_blk_core_alloc ) 
     116      IF( sbc_blk_core_alloc /= 0 )   CALL ctl_warn('sbc_blk_core_alloc: failed to allocate arrays') 
     117   END FUNCTION sbc_blk_core_alloc 
     118 
     119    
    104120   SUBROUTINE sbc_blk_core( kt ) 
    105121      !!--------------------------------------------------------------------- 
     
    151167         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    152168         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    153          &                  sn_tdif, rn_zqt,  rn_zu 
     169         &                  sn_tdif, rn_zqt,  rn_zu    , ln_Cd_L12 
    154170      !!--------------------------------------------------------------------- 
    155171      ! 
     
    157173      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    158174         !                                      ! ====================== ! 
     175         ! 
     176         !                                      ! allocate sbc_blk_core array (clem) 
     177         IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core : unable to allocate standard arrays' ) 
    159178         ! 
    160179         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
     
    319338         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    320339     
     340      ! Make ocean-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     341#if defined key_lim3 
     342      IF( ln_Cd_L12 ) THEN 
     343 
     344         Cd_oce(:,:) = Cd(:,:)    ! record value of pure ocean-atm. drag 
     345 
     346         CALL Cdn10_Lupkes2012( Cd )  ! calculate new drag from Lupkes(2012) equations 
     347 
     348      ENDIF 
     349#endif 
     350       
    321351      ! ... tau module, i and j component 
    322352      DO jj = 1, jpj 
     
    437467      !!--------------------------------------------------------------------- 
    438468      INTEGER  ::   ji, jj    ! dummy loop indices 
    439       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
    440469      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    441470      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     471      REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    442472      !!--------------------------------------------------------------------- 
    443473      ! 
    444474      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
    445475      ! 
    446       ! local scalars ( place there for vector optimisation purposes) 
    447       zcoef_wnorm  = rhoa * Cice 
    448       zcoef_wnorm2 = rhoa * Cice * 0.5 
     476      CALL wrk_alloc( jpi,jpj, Cd ) 
     477 
     478      Cd(:,:) = Cd_ice 
     479       
     480      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     481#if defined key_lim3 
     482      IF( ln_Cd_L12 ) THEN 
     483 
     484         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     485 
     486      ENDIF 
     487#endif 
    449488 
    450489!!gm brutal.... 
     
    467506               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    468507                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    469                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     508               zwnorm_f = rhoa * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    470509               ! ... ice stress at I-point 
    471510               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     
    476515               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    477516                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    478                wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     517               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    479518            END DO 
    480519         END DO 
     
    493532         DO jj = 2, jpjm1 
    494533            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    495                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     534               utau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    496535                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    497                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     536               vtau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    498537                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    499538            END DO 
     
    510549      ENDIF 
    511550 
     551      CALL wrk_dealloc( jpi,jpj, Cd ) 
     552 
    512553      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
    513554       
     
    548589      ! local scalars ( place there for vector optimisation purposes) 
    549590      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    550       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    551       zcoef_dqsb   = rhoa * cpa * Cice 
     591      zcoef_dqla   = -Ls * Cd_ice * 11637800. * (-5897.8) 
     592      zcoef_dqsb   = rhoa * cpa * Cd_ice 
    552593 
    553594      zztmp = 1. / ( 1. - albo ) 
     
    575616               ! ... turbulent heat fluxes 
    576617               ! Sensible Heat 
    577                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     618               z_qsb(ji,jj,jl) = rhoa * cpa * Cd_ice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    578619               ! Latent Heat 
    579                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     620               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cd_ice * wndm_ice(ji,jj)   &                            
    580621                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    581622              ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    820861         ! 
    821862      END DO 
    822  
     863       
    823864      CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
    824865      CALL wrk_dealloc( jpi,jpj, zeta_u, stab ) 
     
    903944   END FUNCTION psi_h 
    904945 
     946 
     947   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     948      !!---------------------------------------------------------------------- 
     949      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
     950      !! 
     951      !! ** Purpose :    Recompute the ice-atm and ocean-atm drags at 10m height to make 
     952      !!                 them dependent on edges at leads, melt ponds and flows. 
     953      !!                 After some approximations, this can be resumed to a dependency 
     954      !!                 on ice concentration. 
     955      !!                 
     956      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
     957      !!                 with the highest level of approximation: level4, eq.(59) 
     958      !!                 The drag can be re-written as follows: 
     959      !! 
     960      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu 
     961      !! 
     962      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59) 
     963      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) 
     964      !!                    A is the concentration of ice minus melt ponds (if any) 
     965      !! 
     966      !!                 This new drag has a parabolic shape (as a function of A) starting at 
     967      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     968      !!                 and going down to Cdi(say 1.4e-3) for A=1 
     969      !! 
     970      !!                 It is theoretically applicable to all ice conditions (not only MIZ) 
     971      !!                 => see Lupkes et al (2013) 
     972      !! 
     973      !! ** References : Lupkes et al. JGR 2012 (theory) 
     974      !!                 Lupkes et al. GRL 2013 (application to GCM) 
     975      !! 
     976      !!---------------------------------------------------------------------- 
     977      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     978      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
     979      REAL(wp), PARAMETER ::   znu   = 1._wp 
     980      REAL(wp), PARAMETER ::   zmu   = 1._wp 
     981      REAL(wp), PARAMETER ::   zbeta = 1._wp 
     982      REAL(wp)            ::   zcoef 
     983      !!---------------------------------------------------------------------- 
     984      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
     985 
     986      Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
     987         &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     988         &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     989       
     990   END SUBROUTINE Cdn10_Lupkes2012 
     991       
    905992   !!====================================================================== 
    906993END MODULE sbcblk_core 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6746 r6763  
    617617      u_ice_b(:,:)     = u_ice(:,:) 
    618618      v_ice_b(:,:)     = v_ice(:,:) 
     619      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 ) 
    619620       
    620621   END SUBROUTINE sbc_lim_bef 
Note: See TracChangeset for help on using the changeset viewer.