Changeset 11951 for NEMO/branches/2019

Timestamp:

2019-11-22T16:52:49+01:00 (4 years ago)

Author:

stefryn

Message:

fixed bug in stress calculation

Location:

NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE

Files:

• ice.F90 (modified) (2 diffs)
• icedyn_rhg_eap.F90 (modified) (12 diffs)

NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/ice.F90

@@ -226 +226 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field                  [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aniso_11, aniso_12   !: structure tensor elements
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdg_conv
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
@@ -409 +410 @@
          &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      &
          &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      &
-         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , STAT=ierr(ii) )
+         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) )
 
       ii = ii + 1

NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90

@@ -67 +67 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, aniso_11, aniso_12 )
+   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12 )
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE ice_dyn_rhg_eap  ***
@@ -122 +122 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   aniso_11, aniso_12   ! structure tensor components
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
       !!
       INTEGER ::   ji, jj       ! dummy loop indices
@@ -141 +141 @@
       REAL(wp) ::   zfac_x, zfac_y
       REAL(wp) ::   zshear, zdum1, zdum2
-      REAL(wp) ::   stressptmp, stressmtmp                              ! anisotropic stress tensor components
+      REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components
       REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution
       !
@@ -404 +404 @@
                ! --- anisotropic stress calculation --- !
                CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, &
-                                       zdsT, aniso_11(ji,jj), aniso_12(ji,jj), &
+                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), &
                                        stressptmp, stressmtmp, &
                                        stress12tmp(ji,jj), strength(ji,jj), &
                                        alphar, alphas)
+               IF (jter == nn_nevp) THEN
+                  rdg_conv(ji,jj) = -min( alphar, 0._wp)    
+               ENDIF
 
                ! delta at T points
-               !zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
 
                ! P/delta at T points
-               !zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
 
                ! alpha & beta for aEVP
@@ -428 +431 @@
                zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1
                zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1
-             
+               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1
+               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1
+               !WRITE(numout,*) 'stress output', ji, jj, zs1(ji,jj)
             END DO
          END DO
@@ -437 +442 @@
             DO ji = 1, jpim1
                ! stress12tmp at F points 
-               zdsT = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
+               stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
                   &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
                   &   ) * 0.25_wp * r1_e1e2f(ji,jj)
@@ -449 +454 @@
                ! stress at F points (zkt/=0 if landfast)
                !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
-               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmp(ji,jj) * zalph1  ) * z1_alph1
+               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1
+               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1
 
             END DO
@@ -749 +755 @@
       END DO
       CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      !CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1. )
       
       ! --- Store the stress tensor for the next time step --- !
@@ -922 +929 @@
       REAL(wp), PARAMETER ::   kfriction = 0.45_wp
       !!---------------------------------------------------------------------
+      ! Factor to maintain the same stress as in EVP (see Section 3)
+      ! Can be set to 1 otherwise
+      invstressconviso = 1._wp/(1._wp+kfriction*kfriction)
+ 
+      invsin = 1._wp/sin(2._wp*phi) * invstressconviso 
+      !now uses phi as set in higher code
+       
+      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
+      ! rates
+
       ! 1) structure tensor
       a22 = 1._wp-a11
@@ -949 +966 @@
       dtemp12 = shear*0.5_wp
       dtemp22 = 0.5_wp*(divu - tension)
+
+      !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear
 
       ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
@@ -1073 +1092 @@
       sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin
 
+      !WRITE(numout,*) 'strength inside loop', strength
+      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22
+
       ! Back - rotation of the stress from principal axes into general coordinates
 
@@ -1083 +1105 @@
       stress12 = sgprm12
       stressm  = sgprm11 - sgprm22
+
+      !WRITE(numout,*) 'stress output inside loop', ksub, stressp
 
       ! Compute ridging and sliding functions in general coordinates