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

Timestamp:

2021-12-06T18:03:30+01:00 (2 years ago)

Author:

edblockley

Message:

Modifying EVP routine to allow calculation of sig?_pnorm

File:

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

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

@@ -856 +856 @@
       ! Recommendation 1 : we use ice strength, not replacement pressure
       ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
-!!$      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 jj = 1, jpj
-!!$            DO ji = 1, jpi
-!!$            
-!!$               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
-!!$               !                        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 ) )
-!!$               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
-!!$               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
-!!$               
-!!$               ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
-!!$      
-!!$               ! Normalized  principal stresses (used to display the ellipse)
-!!$               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
-!!$               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
-!!$               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
-!!$            END DO
-!!$         END DO               
-!!$         !
-!!$         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
-!!$         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
-!!$      
-!!$         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
-!!$         
-!!$      ENDIF
+      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 jj = 1, jpj
+            DO ji = 1, jpi
+            
+               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+               !                        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 ) )
+               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+               
+               ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
+      
+               ! Normalized  principal stresses (used to display the ellipse)
+               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+            END DO
+         END DO               
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+      
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
+         
+      ENDIF
 
       ! --- SIMIP --- !