Changeset 11377


Ignore:
Timestamp:
2019-07-31T14:26:56+02:00 (12 months ago)
Author:
clem
Message:

clean rheology, add a couple of outputs and remove landfast_home option

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/field_def_nemo-ice.xml

    r11371 r11377  
    4444          <field id="icefrb"       long_name="Sea-ice freeboard"                                       standard_name="sea_ice_freeboard"                         unit="m"    /> 
    4545          <field id="icealb"       long_name="Sea-ice or snow albedo"                                  standard_name="sea_ice_albedo"                            unit=""    detect_missing_value="true" /> 
    46           <field id="tau_icebfr"   long_name="ice friction on ocean bottom for landfast ice"                                                                     unit="N/2"  /> 
    4746      
    4847     <!-- melt ponds --> 
     
    7170          <field id="utau_oi"      long_name="X-component of ocean stress on sea ice"                  standard_name="sea_ice_base_upward_x_stress"              unit="N/m2" /> 
    7271          <field id="vtau_oi"      long_name="Y-component of ocean stress on sea ice"                  standard_name="sea_ice_base_upward_y_stress"              unit="N/m2" /> 
     72          <field id="utau_bi"      long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress"              unit="N/m2" /> 
     73          <field id="vtau_bi"      long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress"              unit="N/m2" /> 
    7374          <field id="isig1"        long_name="1st principal stress component for EVP rhg"                                                                        unit=""     /> 
    7475          <field id="isig2"        long_name="2nd principal stress component for EVP rhg"                                                                        unit=""     /> 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref

    r10911 r11377  
    5656   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    5757   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
    58    ln_landfast_home = .false.         !  landfast: parameterization from "home made" 
    5958      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
    6059                                      !          recommended range: [0.1 ; 0.25] - L16=0.125 - home=0.15 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90

    r11371 r11377  
    136136   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    137137   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016  
    138    LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made  
    139138   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice 
    140139   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn.F90

    r11371 r11377  
    221221      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
    222222         &             rn_ishlat ,                                                           & 
    223          &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
     223         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    224224      !!------------------------------------------------------------------- 
    225225      ! 
     
    244244         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
    245245         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
    246          WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home 
    247246         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
    248247         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
     
    271270      ENDIF 
    272271      !                                      !--- Landfast ice 
    273       IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp 
    274       ! 
    275       IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 
    276          CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 
    277       ENDIF 
     272      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp 
    278273      ! 
    279274      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90

    r11317 r11377  
    6969         WRITE(numout,*)'ice_dyn_rhg: sea-ice rheology' 
    7070         WRITE(numout,*)'~~~~~~~~~~~' 
    71       ENDIF 
    72       ! 
    73       IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    74          tau_icebfr(:,:) = 0._wp 
    75          DO jl = 1, jpl 
    76             WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    77          END DO 
    7871      ENDIF 
    7972      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg_evp.F90

    r11371 r11377  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 
    115       !                                                                              or only at the main time step (F) 
    116114      INTEGER ::   ji, jj       ! dummy loop indices 
    117115      INTEGER ::   jter         ! local integers 
     
    137135      ! 
    138136      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    139       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points 
     137      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
    140138      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    141139      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    142       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
    143       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param) 
    144       REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
    145140      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    146       REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses 
    147141      ! 
    148142      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     
    152146      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    153147      !                                                                 !    ice bottom surface if ice is embedded    
    154       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    155       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
    156       ! 
    157       REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays 
    158       REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence 
     148      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
     150      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array 
     151      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points 
     152      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points 
     153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
     154      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
     155      ! 
     156      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
     157      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    159158      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice 
    160159 
     
    254253 
    255254      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    256       IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile 
    257       ELSE                                               ;   zkt = 0._wp 
     255      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     256      ELSE                         ;   zkt = 0._wp 
    258257      ENDIF 
    259258      ! 
     
    298297             
    299298            ! Drag ice-atm. 
    300             zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    301             zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     299            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     300            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
    302301 
    303302            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     
    306305 
    307306            ! masks 
    308             zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
    309             zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     307            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     308            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    310309 
    311310            ! switches 
    312             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp 
    313             ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF 
    314             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp 
    315             ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF 
     311            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     312            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     313            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     314            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
    316315 
    317316         END DO 
     
    329328               ! ice-bottom stress at U points 
    330329               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    331                zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     330               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    332331               ! ice-bottom stress at V points 
    333332               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    334                zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     333               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    335334               ! ice_bottom stress at T points 
    336335               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    337                tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     336               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    338337            END DO 
    339338         END DO 
    340339         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 
    341340         ! 
    342       ELSEIF( ln_landfast_home ) THEN          !-- Home made 
     341      ELSE                               !-- no landfast 
    343342         DO jj = 2, jpjm1 
    344343            DO ji = fs_2, fs_jpim1 
    345                zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 
    346                zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 
    347             END DO 
    348          END DO 
    349          ! 
    350       ELSE                                     !-- no landfast 
    351          DO jj = 2, jpjm1 
    352             DO ji = fs_2, fs_jpim1 
    353                zTauU_ib(ji,jj) = 0._wp 
    354                zTauV_ib(ji,jj) = 0._wp 
     344               ztaux_base(ji,jj) = 0._wp 
     345               ztauy_base(ji,jj) = 0._wp 
    355346            END DO 
    356347         END DO 
     
    493484                  !                 !--- tau_bottom/v_ice 
    494485                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    495                   zTauB = - zTauV_ib(ji,jj) / zvel 
     486                  zTauB = ztauy_base(ji,jj) / zvel 
     487                  !                 !--- OceanBottom-to-Ice stress 
     488                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    496489                  ! 
    497490                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    498                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     491                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    499492                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    500493                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    501494                  ! 
    502495                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    503                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     496                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    504497                  ! 
    505498                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    506                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     499                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    507500                  ! 
    508501                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    509                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    510                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    511                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    512                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    513                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    514                      &           ) * zmaskV(ji,jj) 
     502                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     503                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     504                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     505                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     506                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     507                        &           )   * zmsk00y(ji,jj) 
    515508                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    516                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    517                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    518                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    519                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    520                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    521                      &            ) * zmaskV(ji,jj) 
     509                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     510                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     511                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     512                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     513                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     514                        &            )   * zmsk00y(ji,jj) 
    522515                  ENDIF 
    523516               END DO 
     
    529522            CALL agrif_interp_ice( 'V' ) 
    530523#endif 
    531             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     524            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    532525            ! 
    533526            DO jj = 2, jpjm1 
     
    541534                  !                 !--- tau_bottom/u_ice 
    542535                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    543                   zTauB = - zTauU_ib(ji,jj) / zvel 
     536                  zTauB = ztaux_base(ji,jj) / zvel 
     537                  !                 !--- OceanBottom-to-Ice stress 
     538                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    544539                  ! 
    545540                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    546                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     541                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    547542                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    548543                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    549544                  ! 
    550545                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    551                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     546                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    552547                  ! 
    553548                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    554                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     549                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    555550                  ! 
    556551                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    557                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    558                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    559                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    560                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    561                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    562                      &            ) * zmaskU(ji,jj) 
     552                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     553                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     554                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     555                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     556                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     557                        &            )   * zmsk00x(ji,jj) 
    563558                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    564                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    565                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    566                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    567                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    568                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    569                      &            ) * zmaskU(ji,jj) 
     559                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     560                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     561                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     562                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     563                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     564                        &            )   * zmsk00x(ji,jj) 
    570565                  ENDIF 
    571566               END DO 
     
    577572            CALL agrif_interp_ice( 'U' ) 
    578573#endif 
    579             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     574            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    580575            ! 
    581576         ELSE ! odd iterations 
     
    591586                  !                 !--- tau_bottom/u_ice 
    592587                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    593                   zTauB = - zTauU_ib(ji,jj) / zvel 
     588                  zTauB = ztaux_base(ji,jj) / zvel 
     589                  !                 !--- OceanBottom-to-Ice stress 
     590                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    594591                  ! 
    595592                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    596                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     593                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    597594                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    598595                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    599596                  ! 
    600597                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    601                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     598                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    602599                  ! 
    603600                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    604                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     601                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    605602                  ! 
    606603                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    607                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    608                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    609                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    610                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    611                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    612                      &            ) * zmaskU(ji,jj) 
     604                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     605                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     606                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     607                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     608                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     609                        &            )   * zmsk00x(ji,jj) 
    613610                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    614                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    615                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    616                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    617                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    618                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    619                      &            ) * zmaskU(ji,jj) 
     611                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     612                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     613                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     614                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     615                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     616                        &            )   * zmsk00x(ji,jj) 
    620617                  ENDIF 
    621618               END DO 
     
    627624            CALL agrif_interp_ice( 'U' ) 
    628625#endif 
    629             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     626            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    630627            ! 
    631628            DO jj = 2, jpjm1 
     
    639636                  !                 !--- tau_bottom/v_ice 
    640637                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    641                   zTauB = - zTauV_ib(ji,jj) / zvel 
     638                  zTauB = ztauy_base(ji,jj) / zvel 
     639                  !                 !--- OceanBottom-to-Ice stress 
     640                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    642641                  ! 
    643642                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    644                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     643                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    645644                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    646645                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    647646                  ! 
    648647                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    649                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     648                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    650649                  ! 
    651650                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    652                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     651                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    653652                  ! 
    654653                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    655                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    656                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    657                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    658                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    659                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    660                      &           ) * zmaskV(ji,jj) 
     654                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     655                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     656                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     657                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     658                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     659                        &           )   * zmsk00y(ji,jj) 
    661660                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    662                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    663                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    664                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    665                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    666                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    667                      &            ) * zmaskV(ji,jj) 
     661                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     662                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     663                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     664                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     665                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     666                        &            )   * zmsk00y(ji,jj) 
    668667                  ENDIF 
    669668               END DO 
     
    675674            CALL agrif_interp_ice( 'V' ) 
    676675#endif 
    677             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     676            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    678677            ! 
    679678         ENDIF 
     
    690689      END DO                                              !  end loop over jter  ! 
    691690      !                                                   ! ==================== ! 
    692       ! 
    693       IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN 
    694          CALL bdy_ice_dyn( 'U' ) 
    695          CALL bdy_ice_dyn( 'V' ) 
    696       ENDIF 
    697691      ! 
    698692      !------------------------------------------------------------------------------! 
     
    757751      END DO 
    758752 
    759       ! --- ice-ocean and ice-atm. stresses --- ! 
    760       IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') ) THEN 
    761          CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., & 
    762             &                                  zTauU_ia, 'U', -1., zTauV_ia, 'V', -1. ) 
     753      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
     754      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     755         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
     756         ! 
     757         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 
     758            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     759         ! 
    763760         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
    764761         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 
    765          CALL iom_put( 'utau_ai' , zTauU_ia * zmsk00 ) 
    766          CALL iom_put( 'vtau_ai' , zTauV_ia * zmsk00 ) 
     762         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
     763         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
     764         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
     765         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
    767766      ENDIF 
    768   
    769       ! --- ice-oceanbottom stress in case of landfast --- ! 
    770       IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr * zmsk00 ) 
    771        
     767        
    772768      ! --- divergence, shear and strength --- ! 
    773769      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
     
    817813         ! 
    818814         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 
    819             &                                  zCorx, 'U', -1., zCory, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
     815            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
    820816 
    821817         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
    822818         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
    823          CALL iom_put( 'corstrx' , zCorx * zmsk00 )   ! Coriolis force term in force balance (x) 
    824          CALL iom_put( 'corstry' , zCory * zmsk00 )   ! Coriolis force term in force balance (y) 
     819         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
     820         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
    825821         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
    826822         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
Note: See TracChangeset for help on using the changeset viewer.