Merge branch 'cherry-pick-b6f56da1' into 'main'

Merge branch '22-rheology_diags' into 'main'

Closes #22

See merge request nemo/nemo!46

cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml

@@ -48,11 +48,18 @@
        <!-- rheology -->
        <field field_ref="icediv"           name="sidive"  />
        <field field_ref="iceshe"           name="sishea"  />
+       <field field_ref="icedlt"           name="sidelt"  />
        <field field_ref="icestr"           name="sistre"  />
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
        <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
        <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
+
+       <field field_ref="velo_res"         name="velo_res"  />
+       <field field_ref="velo_ero"         name="velo_ero"  />
+       <field field_ref="uice_eri"         name="uice_eri"  />
+       <field field_ref="vice_eri"         name="vice_eri"  />
+       <field field_ref="uice_cvg"         name="uice_cvg"  />
        
        <!-- heat fluxes -->
        <field field_ref="qt_oce_ai"        name="qt_oce_ai"  />
@@ -64,7 +71,7 @@
        <field field_ref="qns_ice"          name="qns_ice"    />
        <field field_ref="qemp_ice"         name="qemp_ice"   />
        <field field_ref="albedo"           name="albedo"     />
-       
+
        <field field_ref="hfxcndtop"        name="hfxcndtop"  />
        <field field_ref="hfxcndbot"        name="hfxcndbot"  />
        <field field_ref="hfxsensib"        name="hfxsensib"  />

cfgs/SHARED/field_def_nemo-ice.xml

@@ -200,6 +200,12 @@
       <!-- rheology convergence tests -->
       <field id="uice_cvg"   long_name="sea ice velocity convergence"      standard_name="sea_ice_velocity_convergence"      unit="m/s" />
 
+      <!-- vp rheology convergence tests -->
+      <field id="velo_res" long_name="sea ice velocity residual"                        standard_name="sea_ice_velocity_residual"         unit="N/m2" />
+      <field id="velo_ero" long_name="sea ice velocity error last outer iteration"      standard_name="sea_ice_velocity_outer_error"      unit="m/s" />
+      <field id="uice_eri" long_name="uice velocity error last inner iteration"         standard_name="sea_ice_u_velocity_inner_error"    unit="m/s" />
+      <field id="vice_eri" long_name="vice velocity error last inner iteration"         standard_name="sea_ice_v_velocity_inner_error"    unit="m/s" />
+
       <!-- ================= -->
       <!-- Add-ons for SIMIP -->
       <!-- ================= -->
@@ -429,7 +435,7 @@
       <field field_ref="sheastr"          name="sheastr" />
       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
-      <field field_ref="icedlt"           name="sidelta" />
+      <field field_ref="icedlt"           name="sidelt" />
 
       <!-- heat fluxes -->
       <field field_ref="qt_oce_ai"        name="qt_oce_ai"  />

src/ICE/icedyn_rhg_evp.F90

@@ -145,7 +145,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
-      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i, zshear                  ! tension, shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
       !                                                                 !    ocean surface (ssh_m) if ice is not embedded
@@ -749,8 +749,11 @@ CONTAINS
             &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
             &   ) * 0.25_wp * r1_e1e2t(ji,jj)
 
-         ! shear at T points
-         pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj)
+         ! maximum shear rate at T points (includes tension, output only)
+         pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj) ! 
+
+         ! shear at T-points
+         zshear(ji,jj)   = SQRT( zds2 ) * zmsk(ji,jj)
 
          ! divergence at T points
          pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
@@ -758,21 +761,25 @@ CONTAINS
             &             ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
 
          ! delta at T points
-         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
-         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
+         zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
+
+         ! delta* at T points (pdelta_i)
+         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
+         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch  
+                           ! it seems that deformation used for advection and mech redistribution is delta*
+                           ! MV in principle adding creep limit is a regularization for viscosity not for delta
+                           ! delta_star should not (in my view) be used in a replacement for delta
 
       END_2D
+
       CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
-         &                            zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp )
+         &                            zshear  , 'T', 1._wp, zdelta , 'T', 1._wp, zs1     , 'T', 1._wp, zs2   , 'T', 1._wp, zs12, 'F', 1._wp )
 
       ! --- Store the stress tensor for the next time step --- !
       pstress1_i (:,:) = zs1 (:,:)
       pstress2_i (:,:) = zs2 (:,:)
       pstress12_i(:,:) = zs12(:,:)
       !
-
-      !------------------------------------------------------------------------------!
       ! 5) diagnostics
       !------------------------------------------------------------------------------!
       ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
@@ -795,7 +802,7 @@ CONTAINS
       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
-      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , zdelta   * zmsk00 )   ! delta
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
       IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
@@ -805,21 +812,19 @@ CONTAINS
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
             ! Ice stresses
-            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
-            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
-            ! I know, this can be confusing...
-            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
-            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+            ! sigma1, sigma2, sigma12 are some recombination of the stresses (HD MWR002, Bouillon et al., OM2013)
+            ! not to be confused with stress tensor components, stress invariants, or stress principal components
+            zfac             =   strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )          ! viscosity
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
             zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
-            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
 
             ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
-            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                      ! 1st stress invariant, aka average normal stress, aka negative pressure
-            zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )  ! 2nd  ''       ''    , aka maximum shear stress
+            zsig_I (ji,jj)   =   0.5_wp * zsig1 
+            zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 )
 
          END_2D
          !
-         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
          IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
 
@@ -830,24 +835,23 @@ CONTAINS
       ! --- Normalized stress tensor principal components --- !
       ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
       ! Recommendation 1 : we use ice strength, not replacement pressure
-      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
+      ! Recommendation 2 : for EVP, no need to use viscosities at last iteration (stress is properly iterated)
       IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
          !
          ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
-            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
-            !                        and **deformations** at current iterates
+            ! For EVP solvers, ice stresses at current iterates can be used
             !                        following Lemieux & Dupont (2020)
-            zfac             =   zp_delt(ji,jj)
-            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+            zfac             =   strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
             zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
-            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
 
             ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
-            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st stress invariant, aka average normal stress, aka negative pressure
-            zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )   ! 2nd  ''       ''    , aka maximum shear stress
+            zsig_I(ji,jj)    =   0.5_wp * zsig1                                         ! normal stress
+            zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 ) ! max shear stress
 
             ! Normalized  principal stresses (used to display the ellipse)
             z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )

tests/ICE_ADV2D/EXPREF/README

+
+!! comment from MV 
+!! not clear why one has both EXP00 and EXPREF directories
+!! end MV
+
 WARNING: For now, the test case ICE_ADV2D only works if the logical "ll_neg" is set to FALSE in the routine icedyn_adv_umx.F90
                   (it is still unclear why)
 -------
@@ -28,8 +33,8 @@ How to run
 ----------
 
 a) Compile and run the model once to get a mesh_mask.nc file with the following command:
-../../../makenemo -r ICE_ADV2D -n ICE_ADV2D -m X64_ADA -j 4
-mpirun ./nemo -np 1
+./makenemo -n ICE_ADV2D -a ICE_ADV2D -m  X64_JEANZAY -j8
+./nemo.exe (on jeanzaypp)
 
 b) Create the initial condition file for sea-ice (initice.nc) by running this python script: 
 python ./make_INITICE.py

tests/ICE_ADV2D/EXPREF/make_INITICE.py

@@ -15,18 +15,18 @@ fcoord='mesh_mask.nc'
 # output file
 fflx='initice.nc'
 
-print '   creating init ice file  ' +fflx
+print ('   creating init ice file  ' +fflx)
 
 # Reading coordinates file
 nccoord=netcdf(fcoord,'r')
-nav_lon=nccoord.variables['nav_lon']
-nav_lat=nccoord.variables['nav_lat']
+nav_lon=nccoord.variables['glamt']
+nav_lat=nccoord.variables['gphit']
 time_counter=1
 LON1= nav_lon.shape[1]
-LAT1= nav_lon.shape[0]
-print 'nav_lon.shape[1]' ,nav_lon.shape[1]
-print 'LON1 ', LON1
-print 'LAT1 ', LAT1
+LAT1= nav_lon.shape[1]
+print ('nav_lon.shape[1]' ,nav_lon.shape[1])
+print ('LON1 ', LON1)
+print ('LAT1 ', LAT1)
 
 # Creating INITICE netcdf file
 nc=netcdf(fflx,'w')