Merge pull request #19 from TillRasmussen/eap

Eap

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -1939,58 +1939,38 @@ subroutine stepa  (nx_block,   ny_block,       &
          j = indxtj(ij)
 
 ! ne 
-         call calc_ffrac(1, stressp_1(i,j), stressm_1(i,j),  &
-                           stress12_1(i,j),                  &
-                           a11_1(i,j),                       &
-                           mresult11)
-
-         call calc_ffrac(2, stressp_1(i,j), stressm_1(i,j),  &
-                           stress12_1(i,j),                  &
-                           a12_1(i,j),                       &
-                           mresult12)
+        call calc_ffrac(stressp_1(i,j), stressm_1(i,j),  &
+                           stress12_1(i,j),                 &
+                           a11_1(i,j), a12_1(i,j),          &
+                           mresult11, mresult12)
 
          a11_1(i,j) = (a11_1(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
          a12_1(i,j) = (a12_1(i,j)*dtei - mresult12) * dteikth ! implicit
 
 
 ! nw
-         call calc_ffrac(1, stressp_2(i,j), stressm_2(i,j),  &
-                           stress12_2(i,j),                  &
-                           a11_2(i,j),                       &
-                           mresult11)
-
-         call calc_ffrac(2, stressp_2(i,j), stressm_2(i,j),  &
-                           stress12_2(i,j),                  &
-                           a12_2(i,j),                       &
-                           mresult12)
+        call calc_ffrac(stressp_2(i,j), stressm_2(i,j),  &
+                           stress12_2(i,j),                 &
+                           a11_2(i,j), a12_2(i,j),          &
+                           mresult11, mresult12)
 
          a11_2(i,j) = (a11_2(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
          a12_2(i,j) = (a12_2(i,j)*dtei - mresult12) * dteikth ! implicit
 
 ! sw
-         call calc_ffrac(1, stressp_3(i,j), stressm_3(i,j),  &
-                           stress12_3(i,j),                  &
-                           a11_3(i,j),                       &
-                           mresult11)
-
-         call calc_ffrac(2, stressp_3(i,j), stressm_3(i,j),  &
-                           stress12_3(i,j),                  &
-                           a12_3(i,j),                       &
-                           mresult12)
+        call calc_ffrac(stressp_3(i,j), stressm_3(i,j),  &
+                           stress12_3(i,j),                 &
+                           a11_3(i,j), a12_3(i,j),          &
+                           mresult11, mresult12)
 
          a11_3(i,j) = (a11_3(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
          a12_3(i,j) = (a12_3(i,j)*dtei - mresult12) * dteikth ! implicit
 
 ! se
-         call calc_ffrac(1, stressp_4(i,j), stressm_4(i,j),  &
-                           stress12_4(i,j),                  &
-                           a11_4(i,j),                       &
-                           mresult11)
-
-         call calc_ffrac(2, stressp_4(i,j), stressm_4(i,j),  &
-                           stress12_4(i,j),                  &
-                           a12_4(i,j),                       &
-                           mresult12)
+        call calc_ffrac(stressp_4(i,j), stressm_4(i,j),  &
+                           stress12_4(i,j),                 &
+                           a11_4(i,j), a12_4(i,j),          &
+                           mresult11, mresult12)
 
          a11_4(i,j) = (a11_4(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
          a12_4(i,j) = (a12_4(i,j)*dtei - mresult12) * dteikth ! implicit
@@ -2010,19 +1990,17 @@ end subroutine stepa
 ! the ice floe re-orientation due to fracture
 ! Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2
 
-      subroutine calc_ffrac (blockno, stressp, stressm, &
-                             stress12,                  &
-                             a1x,                       &
-                             mresult)
 
-      integer(kind=int_kind), intent(in) :: &
-         blockno
+      subroutine calc_ffrac (stressp, stressm, &
+                             stress12,            &
+                             a1x,  a2x,            &
+                             mresult1, mresult2)
 
       real (kind=dbl_kind), intent(in) :: &
-         stressp, stressm, stress12, a1x
+         stressp, stressm, stress12, a1x, a2x
 
       real (kind=dbl_kind), intent(out) :: &
-         mresult
+         mresult1, mresult2
 
       ! local variables
 
@@ -2042,11 +2020,12 @@ subroutine calc_ffrac (blockno, stressp, stressm, &
        if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
           file=__FILE__, line=__LINE__)
 
-       sigma11 = p5*(stressp+stressm) 
-       sigma12 = stress12 
-       sigma22 = p5*(stressp-stressm) 
+       sigma11 = p5*(stressp+stressm)
+       sigma12 = stress12
+       sigma22 = p5*(stressp-stressm)
 
-       if ((sigma11-sigma22) == c0) then
+!       if ((sigma11-sigma22) == c0) then sigma11-sigma22 == 0 => stressn ==0
+       if (stressm == c0) then
          gamma = p5*(pih)
        else
          gamma = p5*atan2((c2*sigma12),(sigma11-sigma22))
@@ -2068,24 +2047,27 @@ subroutine calc_ffrac (blockno, stressp, stressm, &
 
 ! Pure divergence
        if ((sigma_1 >= c0).and.(sigma_2 >= c0))  then
-         mresult = c0
+         mresult1 = c0
+         mresult2 = c0
 
 ! Unconfined compression: cracking of blocks not along the axial splitting direction
 ! which leads to the loss of their shape, so we again model it through diffusion
        elseif ((sigma_1 >= c0).and.(sigma_2 < c0))  then
-         if (blockno == 1) mresult = kfrac * (a1x - Q12Q12)
-         if (blockno == 2) mresult = kfrac * (a1x + Q11Q12)
+         mresult1 = kfrac * (a1x - Q12Q12)
+         mresult2 = kfrac * (a2x + Q11Q12)
 
 ! Shear faulting
        elseif (sigma_2 == c0) then
-         mresult = c0
+        mresult1 = c0
+        mresult2 = c0
        elseif ((sigma_1 <= c0).and.(sigma_1/sigma_2 <= threshold)) then
-         if (blockno == 1) mresult = kfrac * (a1x - Q12Q12)
-         if (blockno == 2) mresult = kfrac * (a1x + Q11Q12)
+         mresult1 = kfrac * (a1x - Q12Q12)
+         mresult2 = kfrac * (a2x + Q11Q12)
 
 ! Horizontal spalling
-       else 
-         mresult = c0
+       else
+         mresult1 = c0
+         mresult2 = c0
        endif
 
       end subroutine calc_ffrac