Changeset 15742 for NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90

Timestamp:

2022-03-08T16:20:40+01:00 (2 years ago)

Author:

emmafiedler

Message:

Merging in R75 strength changes, diagnostics updates and latest changes

File:

• NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90 (modified) (8 diffs)

NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90

@@ -137 +137 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i, zshear                  ! tension, shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -773 +773 @@
                &   ) * 0.25_wp * r1_e1e2t(ji,jj)
             
-            ! shear at T points
+            ! maximum shear rate at T points (includes tension, output only)
             pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! shear at T-points
+            zshear(ji,jj)   = SQRT( zds2 )
 
             ! divergence at T points
@@ -789 +792 @@
       END DO
       CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., &
-         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. )
+         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1., zshear, 'T', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
@@ -835 +838 @@
                zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(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_II(ji,jj)   =   SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )        ! 2nd  ''       ''    , aka maximum shear stress
                
             END DO
@@ -866 +869 @@
                !                        and **deformations** at current iterates
                !                        following Lemieux & Dupont (2020)
-               zfac             =   zp_delt(ji,jj)
-               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+               zfac             =   strength(ji,jj) / ( pdelta_i(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_II(ji,jj)   =   SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )        ! 2nd  ''       ''    , aka maximum shear stress
       
                ! Normalized  principal stresses (used to display the ellipse)
@@ -882 +885 @@
          END DO               
          !
-         IF( iom_use('sig1_pnorm') )  CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) 
+         IF( iom_use('sig1_pnorm') )  CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) )  
          IF( iom_use('sig2_pnorm') )  CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) 
       
@@ -1006 +1009 @@
 
       ! time
-      it = ( kt - 1 ) * kitermax + kiter
+      it = ( kt - nit000 ) * kitermax + kiter
       
       ! convergence
@@ -1026 +1029 @@
          istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
          ! close file
-         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
+         IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax )   istatus = NF90_CLOSE(ncvgid)
       ENDIF