Changeset 3791 for branches/2012

Timestamp:

2013-02-10T14:16:41+01:00 (11 years ago)

Author:

gm

Message:

dev_MERGE_2012: #1041 : corrected diagnostic for the ice shear rate

File:

• branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 (modified) (8 diffs)

branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

@@ -123 +123 @@
       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
       REAL(wp) ::   zu_ice2, zv_ice1   !
-      REAL(wp) ::   zddc, zdtc, zdst   ! delta on corners and on centre
+      REAL(wp) ::   zddc, zdtc, zzdst   ! delta on corners and on centre
       REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
       REAL(wp) ::   sigma1, sigma2     ! internal ice stress
@@ -147 +147 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zdd   , zdt      ! Divergence and tension at centre of grid cells
       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdst             ! Shear on centre of grid cells
       REAL(wp), POINTER, DIMENSION(:,:) ::   deltat, deltac   ! Delta at centre and corners of grid cells
       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2 
@@ -153 +154 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
                                                               !   ocean surface (ssh_m) if ice is not embedded
-                                                              !   ice top surface if ice is embedded
-      
+                                                              !   ice top surface if ice is embedded   
       !!-------------------------------------------------------------------
 
       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
       CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
-      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
+      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  )
       CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
 
@@ -405 +405 @@
 
                !- Calculate Delta at centre of grid cells
-               zdst       = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
+               zzdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
                   &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
                   &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
@@ -412 +412 @@
                   &         / area(ji,jj)
 
-               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
-               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
+               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )  
+               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl )
+!!gm faster to replace the line above with simply:
+!!                deltat(ji,jj) = MAX( delta, creepl )
+!!gm end 
 
                !-Calculate stress tensor components zs1 and zs2 
@@ -711 +714 @@
                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
 
-               zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)          &
-                  &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &
-                  &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &
-                  &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)         & 
-                  &          )                                             &
-                  &         / area(ji,jj)
-
-               deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 
-                  &                          ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & 
+               zdst(ji,jj) = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
+                  &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    &
+                  &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    &
+                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
+
+               deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   & 
+                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 
                   &                          ) + creepl
 
@@ -732 +733 @@
       !
       ! * Invariants of the stress tensor are required for limitd_me
-      ! accelerates convergence and improves stability
+      !   (accelerates convergence and improves stability)
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
             divu_i (ji,jj) = zdd   (ji,jj)
             delta_i(ji,jj) = deltat(ji,jj)
-            shear_i(ji,jj) = zds   (ji,jj)
+            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) )
          END DO
       END DO
-
-      ! Lateral boundary condition
-      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
+      CALL lbc_lnk( divu_i (:,:), 'T', 1. )      ! Lateral boundary condition
       CALL lbc_lnk( delta_i(:,:), 'T', 1. )
-      CALL lbc_lnk( shear_i(:,:), 'F', 1. )
+      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
 
       ! * Store the stress tensor for the next time step
@@ -790 +789 @@
       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
       CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
-      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
+      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  )
       CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )