Changeset 4205 for branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

Timestamp:

2013-11-14T18:02:06+01:00 (10 years ago)

Author:

clem

Message:

 

File:

• branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 (modified) (5 diffs)

branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

@@ -123 +123 @@
 
       REAL(wp) ::   dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, ecc2       ! square of yield ellipse eccenticity
+      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
       REAL(wp) ::   zu_ice2, zv_ice1   !
@@ -199 +199 @@
             zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
 #if defined key_lim3
-            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp
+            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
 #endif
 #if defined key_lim2
@@ -325 +325 @@
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
       ecc2 = ecc * ecc
+      ecci = 1. / ecc2
 
       !-Initialise stress tensor 
@@ -425 +426 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) &
-                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
-                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
-                  * zpresh(ji,jj) ) )                          &       
-                  &        / ( 1.0 + alphaevp * dtotel )
-
-               zs2(ji,jj) = ( zs2(ji,jj)   &
-                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
-                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
-                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+                  &          / ( 1._wp + alphaevp * dtotel )
+
+               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+               !   &         / ( 1._wp + dtotel )
 
             END DO
@@ -468 +472 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj)      &
-                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
-                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
-                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1.0 + dtotel ) 
 
             END DO ! ji