Changeset 13536


Ignore:
Timestamp:
2020-09-29T11:22:04+02:00 (4 months ago)
Author:
vancop
Message:

VP rheology, first code version finalised

Location:
NEMO/branches/2020/SI3-03_VP_rheology
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/cfgs/SHARED/field_def_nemo-ice.xml

    r13064 r13536  
    7575          <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" /> 
    7676          <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" /> 
    77           <field id="isig1"        long_name="1st principal stress component for EVP rhg"                                                                        unit=""     /> 
    78           <field id="isig2"        long_name="2nd principal stress component for EVP rhg"                                                                        unit=""     /> 
    79           <field id="isig3"        long_name="convergence measure for EVP rheology (must be around 1)"                                                           unit=""     /> 
     77          <field id="sig1_pnorm"   long_name="P-normalized 1st principal stress component"                                                                       unit=""     /> 
     78          <field id="sig2_pnorm"   long_name="P-normalized 2nd principal stress component"                                                                       unit=""     /> 
    8079          <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  /> 
    8180          <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  /> 
     
    398397     <field field_ref="normstr"          name="normstr" /> 
    399398     <field field_ref="sheastr"          name="sheastr" /> 
    400      <field field_ref="isig1"            name="isig1"   /> 
    401      <field field_ref="isig2"            name="isig2"   /> 
    402      <field field_ref="isig3"            name="isig3"   /> 
     399     <field field_ref="sig1_pnorm"       name="sig1_pnorm"   /> 
     400     <field field_ref="sig2_pnorm"       name="sig2_pnorm"   /> 
    403401      
    404402     <!-- heat fluxes --> 
  • NEMO/branches/2020/SI3-03_VP_rheology/cfgs/SHARED/namelist_ice_ref

    r12920 r13536  
    9191&namdyn_rhg     !   Ice rheology 
    9292!------------------------------------------------------------------------------ 
    93    ln_rhg_EVP       = .true.          !  EVP rheology 
     93   ln_rhg_EVP       = .false.         !  EVP rheology 
    9494      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    9595      rn_creepl     =   2.0e-9        !     creep limit [1/s] 
     
    9999                                      !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
    100100   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg) 
     101   ln_rhg_VP        = .true.          !  VP rheology 
     102   nn_nout_vp       = 10              !     number of outer iterations 
     103   nn_ninn_vp       = 1500            !     number of inner iterations 
     104   ln_zebra_vp      = .false.         !     activate zebra solver 
     105   rn_relaxu_vp     = 0.95            !     relaxation factor for u-velocity 
     106   rn_relaxv_vp     = 0.95            !     relaxation factor for v-velocity 
     107   rn_uerr_max_vp   = 0.80            !     maximum error on velocity 
     108   rn_uerr_min_vp   = 1.0e-04         !     velocity to decide convergence 
     109   nn_cvgchk_vp     = 5               !     iteration step for convergence check 
    101110/ 
    102111!------------------------------------------------------------------------------ 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90

    r13522 r13536  
    160160   INTEGER , PUBLIC ::   nn_nout_vp       !: Number of outer iterations 
    161161   INTEGER , PUBLIC ::   nn_ninn_vp       !: Number of inner iterations (linear system solver) 
    162    LOGICAL , PUBLIC ::   ln_smooth_vp     !: zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009) 
    163162   LOGICAL , PUBLIC ::   ln_zebra_vp      !: activate zebra (solve the linear system problem every odd j-band, then one every even one) 
    164163   REAL(wp), PUBLIC ::   rn_relaxu_vp     !: U-relaxation factor (MV: can probably be merged with V-factor once ok) 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90

    r13522 r13536  
    5656      !! 
    5757      !! ** Action  : comupte - ice velocity (u_ice, v_ice) 
    58       !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) 
     58      !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) - EVP only 
    5959      !!                      - shear, divergence and delta (shear_i, divu_i, delta_i) 
    6060      !!-------------------------------------------------------------------- 
     
    8888      CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    ! 
    8989         !                             !---------------------------------------------! 
    90          CALL ice_dyn_rhg_vp ( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 
     90         CALL ice_dyn_rhg_vp ( kt, shear_i, divu_i, delta_i ) 
    9191         !          
    9292      END SELECT 
     
    9494      IF( lrst_ice ) THEN                       !* write EVP fields in the restart file 
    9595         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) 
    96          IF( ln_rhg_VP  )   CALL rhg_lsr_rst( 'WRITE', kt ) 
     96         ! MV note: no restart needed for VP as there is no time equation for stress tensor 
    9797      ENDIF 
    9898      ! 
     
    121121      !! 
    122122      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, &           !- evp 
    123     &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_smooth_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp 
     123    &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp 
    124124      !!------------------------------------------------------------------- 
    125125       
     
    149149         WRITE(numout,*) '         number of outer iterations                        nn_nout_vp    = ', nn_nout_vp 
    150150         WRITE(numout,*) '         number of inner iterations                        nn_ninn_vp    = ', nn_ninn_vp 
    151          WRITE(numout,*) '         tanh viscosity smooothing                         ln_smooth_vp  = ', ln_smooth_vp 
    152151         WRITE(numout,*) '         activate zebra solver                             ln_zebra_vp   = ', ln_zebra_vp 
    153152         WRITE(numout,*) '         relaxation factor for u                           rn_relaxu_vp  = ', rn_relaxu_vp 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90

    r13279 r13536  
    137137      ! 
    138138      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
     139      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
    139140      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    140141      ! 
     
    168169      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    169170      !! --- diags 
    170       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
     171      REAL(wp) ::   zsig1, zsig2, zsig12 
     172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i          
    171173      !! --- SIMIP diags 
    172174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     
    406408               ! delta at T points 
    407409               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     410                
     411               zdeltastar_t(ji,jj) = zdelta + rn_creepl                     ! store delta* at previous iterate (for ellipse computation) 
    408412 
    409413               ! P/delta at T points 
     
    695699                        &           )   * zmsk00y(ji,jj) 
    696700                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     701                                     !---> MV : comment above is misleading as EVP is a purely explicit method 
     702                                     !---> MV : I found code below quite unreadable :-) 
     703                                     !---> MV : two maskings + landfast / no landfast options, I think this is a bit too much 
     704                                     !---> MV : I would suggest to split at least the landfast / masking steps if possible 
     705                                     !---> MV : Should also comment why 1% of ocean velocity is assumed - not sure it makes sense btw 
     706                                     !---> MV : I would say 100% of ocean current makes much more sense for free-drift 
     707                                     !---> MV : of very low concentration sea ice.. 
     708                                     !---> MV : but maybe it is a numerical trick... in brief it's n 
    697709                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    698710                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     
    715727 
    716728         ! convergence test 
    717          IF( ln_rhg_chkcvg )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 
     729         IF( ln_rhg_chkcvg )   CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 
    718730         ! 
    719731         !                                                ! ==================== ! 
     
    745757            zdt2 = zdt * zdt 
    746758             
     759            zten_i(ji,jj) = zdt 
     760             
    747761            ! shear**2 at T points (doc eq. A16) 
    748762            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     
    797811      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    798812 
    799       ! --- stress tensor --- ! 
    800       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    801          ! 
    802          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     813      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     814      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     815         ! 
     816         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    803817         !          
    804          DO jj = 2, jpjm1 
    805             DO ji = 2, jpim1 
    806                zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    807                   &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    808                   &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    809  
    810                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    811  
    812                zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    813  
    814 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    815 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    816 !!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    817 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    818                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    819                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    820                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    821             END DO 
    822          END DO 
    823          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    824          ! 
    825          CALL iom_put( 'isig1' , zsig1 ) 
    826          CALL iom_put( 'isig2' , zsig2 ) 
    827          CALL iom_put( 'isig3' , zsig3 ) 
    828          ! 
    829          ! Stress tensor invariants (normal and shear stress N/m) 
    830          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    831          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    832  
    833          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    834       ENDIF 
    835  
     818         DO jj = 2, jpj - 1 
     819            DO ji = 2, jpi - 1 
     820             
     821               ! Ice stresses 
     822               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     823               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
     824               ! I know, this can be confusing... 
     825               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     826               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     827               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     828               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     829                
     830               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
     831               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure 
     832               zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress 
     833                
     834            END DO 
     835         END DO 
     836 
     837         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) 
     838          
     839         ! 
     840         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
     841         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress 
     842         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     843          
     844         DEALLOCATE ( zsig_I, zsig_II ) 
     845          
     846      ENDIF 
     847 
     848      ! --- Normalized stress tensor principal components --- ! 
     849      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 
     850      ! Recommendation 1 : we use ice strength, not replacement pressure 
     851      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
     852      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     853         ! 
     854         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )          
     855         !          
     856         DO jj = 2, jpj - 1 
     857            DO ji = 2, jpi - 1 
     858             
     859               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     860               !                        and **deformations** at current iterates 
     861               !                        following Lemieux & Dupont (2020) 
     862               zfac             =   zp_delt(ji,jj) 
     863               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 
     864               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     865               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     866                
     867               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     868               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure 
     869               zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress 
     870 
     871               ! Normalized  principal stresses (used to display the ellipse) 
     872               z1_strength      =   1._wp / strength(ji,jj) 
     873               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  - zsig_II(ji,jj) ) * z1_strength 
     874               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     875            END DO 
     876         END DO 
     877          
     878         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 
     879                
     880         ! 
     881         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     882         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     883 
     884         DEALLOCATE( zsig1_p , zsig2_p ) 
     885          
     886      ENDIF 
    836887      ! --- SIMIP --- ! 
    837888      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     
    907958 
    908959 
    909    SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 
     960   SUBROUTINE rhg_cvg_evp( kt, kiter, kitermax, pu, pv, pub, pvb ) 
    910961      !!---------------------------------------------------------------------- 
    911       !!                    ***  ROUTINE rhg_cvg  *** 
     962      !!                    ***  ROUTINE rhg_cvg_evp  *** 
    912963      !!                      
    913       !! ** Purpose :   check convergence of oce rheology 
     964      !! ** Purpose :   check convergence of evp ice rheology 
    914965      !! 
    915966      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity 
     
    935986         IF( lwp ) THEN 
    936987            WRITE(numout,*) 
    937             WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 
    938             WRITE(numout,*) '~~~~~~~' 
     988            WRITE(numout,*) 'rhg_cvg_evp : ice rheology convergence control' 
     989            WRITE(numout,*) '~~~~~~~~~~~' 
    939990         ENDIF 
    940991         ! 
     
    9741025      ENDIF 
    9751026       
    976    END SUBROUTINE rhg_cvg 
     1027   END SUBROUTINE rhg_cvg_evp 
    9771028 
    9781029 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r13522 r13536  
    1 ! TODOLIST 
    2 ! 
    3 ! Define all symbols 
    4 ! - Do viscosity smoothing with sum (differentiable rheology) 
    5 ! - Re-calculate internal force diagnostic (end of the routine) 
    6 ! - Do proper masking of output fileds with ice mass (zmsk00 see EVP routine) 
    7 ! 
    8 ! quality control - compare code to mitGCM 
    9 ! quality control - comparing EVP and VP codes 
    10 !  
     1! TODOLOIST 
    112! 
    123! Commit 
    134! Compile 
    145! 
     6! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine) 
     7! - QC:compare code to mitGCM 
     8! - QC: comparing EVP and VP codes  
     9! 
    1510! Clarify 
    1611! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 
    17 ! --- Swap of independent term in the U-V solvers ? 
    18 ! --- Is stress tensor for restart needed ? 
    19 ! --- Is stress tensor calculated at the end of the time step 
     12! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe 
     13! --- Is stress tensor for restart needed ? No!!! 
    2014! 
    2115! Test 
    22 ! --- Can we add the 15% mask in the convergence criteria ? 
    2316! --- Try ADI for u-v solver 
    2417! 
    2518! Add 
    2619! - Charge ellipse as an output (good one from Lemieux and Dupont) 
    27 ! - Bottom drag 
     20! - Bottom drag (landfast ice) 
    2821! - Tensile strength 
    2922! - agrif 
     
    5245   !!---------------------------------------------------------------------- 
    5346   !!   ice_dyn_rhg_vp : computes ice velocities from VP rheolog with LSR solvery 
    54    !!   rhg_vp_rst     : read/write EVP fields in ice restart 
    5547   !!---------------------------------------------------------------------- 
    5648   USE phycst         ! Physical constants 
     
    7971 
    8072   PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90 
    81    PUBLIC   rhg_vp_rst       ! called by icedyn_rhg.F90 
    8273 
    8374   !! for convergence tests 
     
    9990CONTAINS 
    10091 
    101    SUBROUTINE ice_dyn_rhg_vp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) 
     92   SUBROUTINE ice_dyn_rhg_vp( kt, pshear_i, pdivu_i, pdelta_i ) 
    10293      !!------------------------------------------------------------------- 
    10394      !! 
     
    164155      !! 
    165156      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step 
    166       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   ! 
    167157      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    168158      !! 
     
    180170      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    181171      REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume 
    182       REAL(wp) ::   zdeltat, zdeltat_star, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 
     172      REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2              ! temporary scalars 
    183173      REAL(wp) ::   zp_deltastar_f 
    184174      REAL(wp) ::   zfac, zfac1, zfac2, zfac3 
     
    192182      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step 
    193183      ! 
    194       REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta at T points 
     184      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
     185      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points 
    195186      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points 
    196187      REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point 
     
    205196      ! 
    206197      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    207       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    208198      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    209199      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
     
    225215!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
    226216     ! 
    227       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
    228       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
     217      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! mask for lots of ice (1), little ice (0) 
     218      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence (1), no ice (0) 
    229219      ! 
    230220      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
     
    232222      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    233223      !! --- diags 
    234       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3          
    235       !! --- SIMIP diags 
    236       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    237       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
    238       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 
    239       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 
    240       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
    241       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)       
     224      REAL(wp) ::   zsig1, zsig2, zsig12 
     225      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components 
     226      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
     227      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP) 
     228      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP) 
     229      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP) 
     230      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP) 
     231      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s, SIMIP) 
     232      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s, SIMIP)       
    242233                         
    243234      !!---------------------------------------------------------------------------------------------------------------------- 
     
    269260      ecc2    = rn_ecc * rn_ecc 
    270261      z1_ecc2 = 1._wp / ecc2 
    271           
    272       ! Initialise stress tensor from previous time step 
    273       zs1 (:,:) = pstress1_i (:,:)  
    274       zs2 (:,:) = pstress2_i (:,:) 
    275       zs12(:,:) = pstress12_i(:,:) 
    276        
     262                
    277263      ! Initialise convergence checks 
    278264      IF( ln_rhg_chkcvg ) THEN 
     
    389375            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 
    390376 
    391             ! masks 
     377            ! Mask for ice presence (1) or absence (0) 
    392378            zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
    393379            zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    394380 
    395             ! switches  
     381            ! Mask for lots of ice (1) or little ice (0) 
    396382            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    397383            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    398384            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
    399             ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
     385            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    400386 
    401387         END DO 
     
    473459                
    474460               ! delta* at T points (following Lemieux and Dupont, GMD 2020) 
    475                zdeltat_star = MAX( zdeltat, rn_delta_min ) 
    476                 
    477                IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min 
    478                 
     461               zdeltastar_t(ji,jj) = zdeltat + rn_creepl 
     462               
    479463               ! P/delta at T-points 
    480                zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltat_star 
     464               zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj) 
    481465                
    482466               ! Temporary zzt and zet factors at T-points 
     
    486470            END DO 
    487471         END DO 
    488  
     472          
    489473         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 
    490474 
     
    502486          
    503487         CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 
    504           
     488 
    505489         !--------------------------------------------------- 
    506490         ! -- Ocean-ice drag and Coriolis RHS contributions 
     
    799783                        !--------------- 
    800784                        DO ji = 2, jpi - 1     
     785                         
    801786                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) 
     787                            
     788                           ! velocity masking for little-ice and no-ice cases 
     789                           ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge 
     790                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &  
     791                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             & 
     792                                          + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
     793                            
    802794                        END DO 
    803795 
     
    871863                        DO jj = 2, jpj - 1 
    872864                           v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) 
     865 
     866                           ! mask velocity for no-ice and little-ice cases                            
     867                           v_ice(ji,jj) =   zmsk00y(ji,jj)                                        &  
     868                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)             & 
     869                                          + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
     870 
    873871                        END DO 
    874872 
     
    897895                  DO jj = 2, jpj - 1 
    898896                     DO ji = 2, jpi - 1 
    899                         zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) ! * zmask15 ---> MV TEST mask at 15% concentration 
     897                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 
    900898                     END DO 
    901899                  END DO 
     
    975973      !------------------------------------------------------------------------------! 
    976974      ! 
    977       ! OPT: subroutinize ? 
     975      ! MV OPT: subroutinize ? 
    978976             
    979977      DO jj = 1, jpj - 1 
     
    997995            zdt2 = zdt * zdt 
    998996             
     997            zten_i(ji,jj) = zdt 
     998             
    999999            ! shear**2 at T points (doc eq. A16) 
    10001000            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     
    10131013            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    10141014            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    1015             pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     1015             
     1016            !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     1017            pdelta_i(ji,jj) = zdelta + rn_creepl 
    10161018 
    10171019         END DO 
     
    10191021       
    10201022      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
    1021       CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    1022        
    1023       ! --- Store the stress tensor for the next time step --- ! 
    1024       ! MV OPT: Are they computed at the end of the tucking fime step ??? 
    1025       ! MV Is stress at previous time step needed for VP, normally no, because they equation is not tucking fime dependent!!! 
    1026       ! 
    1027       pstress1_i (:,:) = zs1 (:,:) 
    1028       pstress2_i (:,:) = zs2 (:,:) 
    1029       pstress12_i(:,:) = zs12(:,:) 
    1030       ! 
    1031  
     1023       
    10321024      !------------------------------------------------------------------------------! 
    10331025      ! 
     
    10351027      ! 
    10361028      !------------------------------------------------------------------------------! 
    1037       ! MV OPT: subroutinize 
    1038       ! 
    1039       ! --- Ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
     1029      ! 
     1030      ! MV OPT: subroutinize ? 
     1031      ! 
     1032      !---------------------------------- 
     1033      ! --- Recompute stresses if needed  
     1034      !---------------------------------- 
     1035      ! 
     1036      ! ---- Sea ice stresses at T-points 
     1037      IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     1038       
     1039         DO jj = 2, jpj - 1 
     1040            DO ji = 2, jpi - 1 
     1041               zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj)  
     1042               zfac                    =   zp_deltastar_t(ji,jj)  
     1043               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     1044               zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj) 
     1045               zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     1046            END DO 
     1047         END DO 
     1048 
     1049         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. ) 
     1050       
     1051      ENDIF 
     1052       
     1053      ! ---- s12 at F-points       
     1054      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     1055 
     1056         DO jj = 1, jpj - 1 
     1057            DO ji = 1, jpi - 1 
     1058             
     1059               ! P/delta* at F points 
     1060               zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) ) 
     1061                
     1062               ! s12 at F-points  
     1063               zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj) 
     1064                
     1065            END DO 
     1066         END DO 
     1067 
     1068         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. ) 
     1069          
     1070      ENDIF 
     1071      ! 
     1072      !----------------------- 
     1073      ! --- Store diagnostics 
     1074      !----------------------- 
     1075      ! 
     1076      ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- ! 
    10401077      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
    10411078         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
    1042  
    10431079         ! 
    10441080         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 
    1045             &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     1081            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
    10461082         ! 
    10471083         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
     
    10591095      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    10601096 
    1061       ! --- Stress tensor --- ! 
    1062       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    1063          ! 
    1064          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     1097      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     1098      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     1099         ! 
     1100         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags. 
     1101         ! Definitions following Coon (1974) and Feltham (2008) 
     1102         ! 
     1103         ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     1104         ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components 
     1105         ! 
     1106         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    10651107         !          
    10661108         DO jj = 2, jpj - 1 
    10671109            DO ji = 2, jpi - 1 
    1068              
    1069                zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    1070                   &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    1071                   &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    1072  
    1073                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    1074  
    1075                zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    1076  
    1077                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    1078                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    1079                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    1080                 
     1110               ! Stress invariants 
     1111               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                        ! 1st invariant, aka average normal stress aka negative pressure 
     1112               zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress                
    10811113            END DO 
    10821114         END DO 
    10831115 
    1084          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    1085          ! 
    1086          CALL iom_put( 'isig1' , zsig1 ) 
    1087          CALL iom_put( 'isig2' , zsig2 ) 
    1088          CALL iom_put( 'isig3' , zsig3 ) 
    1089          ! 
    1090          ! Stress tensor invariants (normal and shear stress N/m) 
    1091          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    1092          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    1093  
    1094          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     1116         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) 
     1117          
     1118         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress 
     1119         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     1120          
     1121         DEALLOCATE ( zsig_I, zsig_II ) 
     1122          
     1123      ENDIF 
     1124 
     1125      ! --- Normalized stress tensor principal components --- ! 
     1126      ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020) 
     1127      ! To plot the yield curve and evaluate physical convergence, they have two recommendations 
     1128      ! Recommendation 1 : Use ice strength, not replacement pressure 
     1129      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765) 
     1130      ! R2 means we need to recompute stresses 
     1131      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     1132         ! 
     1133         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )          
     1134         !          
     1135         DO jj = 2, jpj - 1 
     1136            DO ji = 2, jpi - 1 
     1137               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     1138               !                        and **deformations** at current iterates 
     1139               !                        following Lemieux & Dupont (2020) 
     1140               zfac             =   zp_deltastar_t(ji,jj) 
     1141               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 
     1142               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     1143               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     1144                
     1145               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     1146               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant 
     1147               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd invariant 
     1148 
     1149               ! Normalized  principal stresses (used to display the ellipse) 
     1150               z1_strength      =   1._wp / strength(ji,jj) 
     1151               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  + zsig_II(ji,jj) ) * z1_strength 
     1152               zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     1153            END DO 
     1154         END DO 
     1155         ! 
     1156         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 
     1157         !       
     1158         ! 
     1159         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     1160         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     1161 
     1162         DEALLOCATE( zsig1_p , zsig2_p ) 
    10951163          
    10961164      ENDIF 
     
    10981166      ! --- SIMIP, terms of tendency for momentum equation  --- ! 
    10991167      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
    1100          & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
    1101          ! 
    1102 !!!!!!!!!         ATTENTION LA FORCE INTERNE DOIT ETTRE RECALCIULEEE ICCI !!!!!!!!!!!!!!!! 
     1168         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN 
     1169         ! 
    11031170         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 
    1104             &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
     1171            &                                 zCorU, 'U', -1., zCorV, 'V', -1. ) 
    11051172 
    11061173         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
     
    11081175         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
    11091176         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
     1177      ENDIF 
     1178       
     1179      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     1180 
     1181 
     1182         ! Recalculate internal forces (divergence of stress tensor) 
     1183         DO jj = 2, jpj - 1 
     1184            DO ji = 2, jpi - 1 
     1185               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     1186                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     1187                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     1188                  &                  + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     1189                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     1190                  &                  ) * r1_e1e2u(ji,jj) 
     1191               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     1192                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     1193                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     1194                  &                  + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     1195                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     1196                  &                  ) * r1_e1e2v(ji,jj) 
     1197            END DO 
     1198         END DO 
     1199             
     1200         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. ) 
     1201          
    11101202         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
    11111203         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
     1204          
    11121205      ENDIF 
    11131206 
     1207      ! --- Ice & snow mass and ice area transports 
    11141208      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 
    11151209         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 
     
    11371231 
    11381232         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 
    1139             &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 
    1140             &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. ) 
     1233            &                                 zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 
     1234            &                                 zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. ) 
    11411235 
    11421236         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s) 
     
    11721266                  &       prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 
    11731267    
    1174 !                  CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
    1175 !                          &        zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
    11761268      !!---------------------------------------------------------------------- 
    11771269      !!                    ***  ROUTINE rhg_cvg_vp  *** 
     
    13021394    
    13031395 
    1304  
    1305    SUBROUTINE rhg_vp_rst( cdrw, kt ) 
    1306       !!--------------------------------------------------------------------- 
    1307       !!                   ***  ROUTINE rhg_vp_rst  *** 
    1308       !!                      
    1309       !! ** Purpose :   Read or write RHG file in restart file 
    1310       !! 
    1311       !! ** Method  :   use of IOM library 
    1312       !!---------------------------------------------------------------------- 
    1313       CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    1314       INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
    1315       ! 
    1316       INTEGER  ::   iter            ! local integer 
    1317       INTEGER  ::   id1, id2, id3   ! local integers 
    1318       !!---------------------------------------------------------------------- 
    1319       ! 
    1320       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize 
    1321          !                                   ! --------------- 
    1322          IF( ln_rstart ) THEN                   !* Read the restart file 
    1323             ! 
    1324             id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. ) 
    1325             id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) 
    1326             id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) 
    1327             ! 
    1328             IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist 
    1329                CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  ) 
    1330                CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
    1331                CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 
    1332             ELSE                                     ! start rheology from rest 
    1333                IF(lwp) WRITE(numout,*) 
    1334                IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0' 
    1335                stress1_i (:,:) = 0._wp 
    1336                stress2_i (:,:) = 0._wp 
    1337                stress12_i(:,:) = 0._wp 
    1338             ENDIF 
    1339          ELSE                                   !* Start from rest 
    1340             IF(lwp) WRITE(numout,*) 
    1341             IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0' 
    1342             stress1_i (:,:) = 0._wp 
    1343             stress2_i (:,:) = 0._wp 
    1344             stress12_i(:,:) = 0._wp 
    1345          ENDIF 
    1346          ! 
    1347       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    1348          !                                   ! ------------------- 
    1349          IF(lwp) WRITE(numout,*) '---- rhg-rst ----' 
    1350          iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
    1351          ! 
    1352          CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  ) 
    1353          CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  ) 
    1354          CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 
    1355          ! 
    1356       ENDIF 
    1357       ! 
    1358    END SUBROUTINE rhg_vp_rst 
    1359  
    13601396    
    13611397#else 
Note: See TracChangeset for help on using the changeset viewer.