Changeset 13549 for NEMO

Timestamp:

2020-10-01T12:01:26+02:00 (4 years ago)

Author:

clem

Message:

4.0-HEAD: remove one communication in rheology, so about 100 overall.

File:

• NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90 (modified) (8 diffs)

NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90

@@ -127 +127 @@
       REAl(wp) ::   zbetau, zbetav
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
-      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
+      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
       REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
@@ -136 +136 @@
       REAL(wp) ::   zshear, zdum1, zdum2
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
@@ -143 +143 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
@@ -372 +372 @@
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+         DO jj = 1, jpjm1
             DO ji = 1, jpim1
 
@@ -382 +382 @@
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
-
-         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
-            DO ji = 2, jpi ! no vector loop
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
 
                ! shear**2 at T points (doc eq. A16)
@@ -405 +404 @@
                
                ! delta at T points
-               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
-
-               ! P/delta at T points
-               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+               zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1._wp )
+         
+         ! P/delta at T points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            END DO
+         END DO
+
+         DO jj = 2, jpj    ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
+            DO ji = 2, jpi ! no vector loop
+
+               ! divergence at T points (duplication to avoid communications)
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               
+               ! tension at T points (duplication to avoid communications)
+               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e1e2t(ji,jj)
 
                ! alpha for aEVP
@@ -426 +446 @@
                
                ! stress at T points (zkt/=0 if landfast)
-               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
-               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+               zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
              
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
 
          ! Save beta at T-points for further computations
@@ -759 +778 @@
             
             ! delta at T points
-            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
-            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
-            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+            zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch
 
          END DO