Changeset 3187

Timestamp:

2011-11-28T17:44:46+01:00 (12 years ago)

Author:

spickles2

Message:

Stephen Pickles, 28 Nov 2011.
First commit of dCSE NEMO project work, part 1 - index re-ordering,
OPA_SRC top level only. Includes fix for sub-optimal auto-partitioning
in nemogcm.F90.

Location:

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC

Files:

• lib_fortran.F90 (modified) (8 diffs)
• nemogcm.F90 (modified) (3 diffs)
• oce.F90 (modified) (1 diff)
• oce_ftrans.h90 (added)
• step.F90 (modified) (5 diffs)
• step_oce.F90 (modified) (1 diff)
• stpctl.F90 (modified) (5 diffs)
• trc_oce.F90 (modified) (2 diffs)
• trc_oce_ftrans.h90 (added)
• vectopt_loop_substitute.h90 (modified) (1 diff)
• wrk_nemo.F90 (modified) (1 diff)
• wrk_nemo_ftrans.h90 (added)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90

@@ -28 +28 @@
 
    INTERFACE glob_sum
-      MODULE PROCEDURE glob_sum_2d, glob_sum_3d,glob_sum_2d_a, glob_sum_3d_a 
+      MODULE PROCEDURE glob_sum_2d, glob_sum_3d, glob_sum_2d_a, glob_sum_3d_a 
    END INTERFACE
 
@@ -39 +39 @@
 #endif
 
+   !! * Control permutation of array indices
+#  include "dom_oce_ftrans.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -69 +72 @@
       !! ** Purpose : perform a masked sum on the inner global domain of a 3D array
       !!-----------------------------------------------------------------------
+!FTRANS ptab :I :I :z 
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab       ! input 3D array
       REAL(wp)                               ::   glob_sum   ! global masked sum
       !!
       INTEGER :: jk
+#if defined key_z_first
+      INTEGER :: ji, jj
+      REAL(wp) :: ztmask
+#endif
       !!-----------------------------------------------------------------------
       !
       glob_sum = 0.e0
+#if defined key_z_first
+      DO jj = 1, jpj
+        DO ji = 1, jpi
+          ztmask = tmask_i(ji,jj)
+          DO jk = 1, jpk
+            glob_sum = glob_sum + ptab(ji,jj,jk)*ztmask
+          END DO
+        END DO
+      END DO
+#else
       DO jk = 1, jpk
          glob_sum = glob_sum + SUM( ptab(:,:,jk)*tmask_i(:,:) )
       END DO
+#endif
       IF( lk_mpp )   CALL mpp_sum( glob_sum )
       !
@@ -107 +126 @@
       !! ** Purpose : perform a masked sum on the inner global domain of two 3D array
       !!-----------------------------------------------------------------------
+!FTRANS ptab1 :I :I :z 
+!FTRANS ptab2 :I :I :z 
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2   ! input 3D array
       REAL(wp)            , DIMENSION(2)     ::   glob_sum       ! global masked sum
       !!
       INTEGER :: jk
+#if defined key_z_first
+      INTEGER :: ji, jj
+      REAL(wp) :: ztmask
+#endif
       !!-----------------------------------------------------------------------
       !
       glob_sum(:) = 0.e0
+#if defined key_z_first
+      DO jj = 1, jpj
+        DO ji = 1, jpi
+          ztmask = tmask_i(ji,jj)
+          DO jk = 1, jpk
+            glob_sum(1) = glob_sum(1) + ptab1(ji,jj,jk)*ztmask
+            glob_sum(2) = glob_sum(2) + ptab2(ji,jj,jk)*ztmask
+          END DO
+        END DO
+      END DO
+#else
       DO jk = 1, jpk
          glob_sum(1) = glob_sum(1) + SUM( ptab1(:,:,jk)*tmask_i(:,:) )
          glob_sum(2) = glob_sum(2) + SUM( ptab2(:,:,jk)*tmask_i(:,:) )
       END DO
+#endif
       IF( lk_mpp )   CALL mpp_sum( glob_sum, 2 )
       !
@@ -161 +198 @@
       !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab
-      REAL(wp)                                     ::   glob_sum   ! global masked sum
+!FTRANS ptab :I :I :z 
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_sum   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
@@ -171 +209 @@
       ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+#if defined key_z_first
+      DO jj = 1, jpj
+         DO ji =1, jpi
+            DO jk = 1, jpk
+#else
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji =1, jpi
+#endif
             ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
@@ -221 +265 @@
       !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab1, ptab2
-      REAL(wp)                                     ::   glob_sum   ! global masked sum
+      REAL(wp), INTENT(in), DIMENSION(:,:,:)   ::   ptab1, ptab2
+!FTRANS ptab1 :I :I :z 
+!FTRANS ptab2 :I :I :z 
+      REAL(wp)                                 ::   glob_sum   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
@@ -231 +277 @@
       ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+#if defined key_z_first
+      DO jj = 1, jpj
+         DO ji =1, jpi
+            DO jk = 1, jpk
+#else
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji =1, jpi
+#endif
             ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -37 +37 @@
    !!   nemo_alloc     : dynamical allocation
    !!   nemo_partition : calculate MPP domain decomposition
-   !!   factorise      : calculate the factors of the no. of MPI processes
+   !!   sqfact         : calculate factors of the no. of MPI processes
    !!----------------------------------------------------------------------
    USE step_oce        ! module used in the ocean time stepping module
@@ -508 +508 @@
       !!                 ***  ROUTINE nemo_partition  ***
       !!
-      !! ** Purpose :   
-      !!
+      !! ** Purpose : Work out a sensible factorisation of the number of
+      !!              processors for the x and y dimensions.
       !! ** Method  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
       !
-      INTEGER, PARAMETER :: nfactmax = 20
-      INTEGER :: nfact ! The no. of factors returned
-      INTEGER :: ierr  ! Error flag
-      INTEGER :: ji
-      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
-      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
-      !!----------------------------------------------------------------------
-
-      ierr = 0
-
-      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
-
-      IF( nfact <= 1 ) THEN
+      INTEGER             :: ifact1, ifact2 ! factors of num_pes, ifact1 <= ifact2
+      !!----------------------------------------------------------------------
+
+      ! Factorise the number of processors into ifact1*ifact2, such that
+      ! ifact1 and ifact2 are as nearly equal as possible.
+
+      CALL sqfact( num_pes, ifact1, ifact2 )
+
+      ! Make sure that the smaller dimension of the processor grid
+      ! is given the smaller dimension of the global domain
+      IF( jpiglo <= jpjglo) THEN
+         jpni = ifact1
+         jpnj = ifact2
+      ELSE
+         jpni = ifact2
+         jpnj = ifact1
+      ENDIF
+
+      ! This should never happen
+      IF( (jpni*jpnj) /= num_pes) THEN
+         WRITE (numout, *) 'WARNING: internal error - factorisation of number of PEs failed'
+      ENDIF
+
+      ! This should only happen if num_pes is prime
+      IF( ifact1 <= 1 ) THEN
          WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
-         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
-         jpnj = 1
-         jpni = num_pes
-      ELSE
-         ! Search through factors for the pair that are closest in value
-         mindiff = 1000000
-         imin    = 1
-         DO ji = 1, nfact-1, 2
-            idiff = ABS( ifact(ji) - ifact(ji+1) )
-            IF( idiff < mindiff ) THEN
-               mindiff = idiff
-               imin = ji
-            ENDIF
-         END DO
-         jpnj = ifact(imin)
-         jpni = ifact(imin + 1)
+         WRITE (numout, *) '       : using grid of ',jpni,' x ',jpnj
       ENDIF
       !
@@ -550 +547 @@
    END SUBROUTINE nemo_partition
 
-
-   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE factorise  ***
-      !!
-      !! ** Purpose :   return the prime factors of n.
-      !!                knfax factors are returned in array kfax which is of 
-      !!                maximum dimension kmaxfax.
-      !! ** Method  :
-      !!----------------------------------------------------------------------
-      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
-      INTEGER                    , INTENT(  out) ::   kerr, knfax
-      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
-      !
-      INTEGER :: ifac, jl, inu
-      INTEGER, PARAMETER :: ntest = 14
-      INTEGER :: ilfax(ntest)
-
-      ! lfax contains the set of allowed factors.
-      data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  &
-         &                            128,   64,   32,   16,    8,   4,   2  /
-      !!----------------------------------------------------------------------
-
-      ! Clear the error flag and initialise output vars
-      kerr = 0
-      kfax = 1
-      knfax = 0
-
-      ! Find the factors of n.
-      IF( kn == 1 )   GOTO 20
-
-      ! nu holds the unfactorised part of the number.
-      ! knfax holds the number of factors found.
-      ! l points to the allowed factor list.
-      ! ifac holds the current factor.
-
-      inu   = kn
-      knfax = 0
-
-      DO jl = ntest, 1, -1
-         !
-         ifac = ilfax(jl)
-         IF( ifac > inu )   CYCLE
-
-         ! Test whether the factor will divide.
-
-         IF( MOD(inu,ifac) == 0 ) THEN
-            !
-            knfax = knfax + 1            ! Add the factor to the list
-            IF( knfax > kmaxfax ) THEN
-               kerr = 6
-               write (*,*) 'FACTOR: insufficient space in factor array ', knfax
-               return
-            ENDIF
-            kfax(knfax) = ifac
-            ! Store the other factor that goes with this one
-            knfax = knfax + 1
-            kfax(knfax) = inu / ifac
-            !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
-         ENDIF
-         !
-      END DO
-
-   20 CONTINUE      ! Label 20 is the exit point from the factor search loop.
-      !
-   END SUBROUTINE factorise
+   SUBROUTINE sqfact ( kn, kna, knb )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE sqfact  ***
+      !!
+      !! ** Purpose :   return factors (kna, knb) of kn, such that
+      !!                (1) kna*knb=kn
+      !!                (2) kna and knb are as near equal as possible
+      !!                (3) kna < knb
+      !! ** Method  :   Search backwards from the square root of kn,
+      !!                until we find an integer that cleanly divides kn
+      !! ** Preconditions : kn must be positive
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kn
+      INTEGER, INTENT(  out) ::   kna, knb
+       
+      ! Search backwards from the square root of n.
+
+      fact_loop: DO kna=SQRT(REAL(kn)),1,-1
+        IF ( kn/kna*kna == kn ) THEN
+          EXIT fact_loop
+        ENDIF
+      END DO fact_loop
+
+      IF( kna < 1 ) kna = 1 
+
+      ! kna divides kn cleanly. Work out the other factor.
+      knb = kn/kna
+
+   END SUBROUTINE sqfact
 
    !!======================================================================

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/oce.F90

@@ -47 +47 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point
 
+   !! * Control permutation of array indices
+#  include "oce_ftrans.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -23 +23 @@
    !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
+   !!                 !  2011-05  (S. Pickles) dCSE NEMO optimisations - z index first
    !!----------------------------------------------------------------------
 
@@ -41 +42 @@
 
    PUBLIC   stp   ! called by opa.F90
+
+   !! * Control permutation of array indices
+   !! DCSE_NEMO: warning! dom_oce and zdf_oce public variables are made available
+   !! through the use of step_oce
+#  include "dom_oce_ftrans.h90"
+#  include "zdf_oce_ftrans.h90"
 
    !! * Substitutions
@@ -77 +84 @@
       INTEGER ::   jk       ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
+#if defined key_z_first
+      INTEGER ::   ji, jj   ! dummy loop indices
+#endif
       !! ---------------------------------------------------------------------
 
@@ -126 +136 @@
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
+#if defined key_z_first
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+               DO jk = 2, nkrnf 
+                  avt(ji,jj,jk) = avt(ji,jj,jk) + 2.e0 * rn_avt_rnf * rnfmsk(ji,jj) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+#else
+         DO jk = 2, nkrnf
+            avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)
+         END DO
+#endif
       ENDIF
       IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
@@ -235 +257 @@
 
       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
-      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics
+                               !                            ! (call after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

@@ -28 +28 @@
    USE trabbc           ! bottom boundary condition        (tra_bbc routine)
    USE trabbl           ! bottom boundary layer            (tra_bbl routine)
-   USE tradmp           ! internal damping                 (tra_dmp routine)
+   USE tradmp, ONLY : lk_tradmp, tra_dmp_init, tra_dmp
+                        ! internal damping                 (tra_dmp routine)
    USE traadv           ! advection scheme control     (tra_adv_ctl routine)
    USE traldf           ! lateral mixing                   (tra_ldf routine)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90

@@ -27 +27 @@
 
    PUBLIC stp_ctl           ! routine called by step.F90
+
+   !! * Control permutation of array indices
+#  include "oce_ftrans.h90"
+#  include "dom_oce_ftrans.h90"
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -67 +72 @@
       ENDIF
 
-      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
-      IF(lwp) REWIND( numstp )                       !  --------------------------
+!! DCSE_NEMO: commenting out these two lines. Do they mess up the profile?
+!     IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
+!     IF(lwp) REWIND( numstp )                       !  --------------------------
 
       !                                              !* Test maximum of velocity (zonal only)
@@ -74 +80 @@
       !! zumax = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
       zumax = 0.e0
+#if defined key_z_first
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 1, jpk
+#else
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
+#endif
                zumax = MAX(zumax,ABS(un(ji,jj,jk)))
           END DO 
@@ -112 +124 @@
       DO jj = 2, jpjm1
          DO ji = 1, jpi
+#if defined key_z_first
+            IF( tmask_1(ji,jj) == 1)  zsmin = MIN(zsmin,sn(ji,jj,1))
+#else
             IF( tmask(ji,jj,1) == 1) zsmin = MIN(zsmin,sn(ji,jj,1))
+#endif
          END DO
       END DO
@@ -121 +137 @@
       IF( zsmin < 0.) THEN 
          IF (lk_mpp) THEN
+#if defined key_z_first
+            CALL mpp_minloc ( sn(:,:,1),tmask_1(:,:), zsmin, ii,ij )
+#else
             CALL mpp_minloc ( sn(:,:,1),tmask(:,:,1), zsmin, ii,ij )
+#endif
          ELSE
+#if defined key_z_first
+            ilocs = MINLOC( sn(:,:,1), mask = tmask_1(:,:) == 1.e0 )
+#else
             ilocs = MINLOC( sn(:,:,1), mask = tmask(:,:,1) == 1.e0 )
+#endif
             ii = ilocs(1) + nimpp - 1
             ij = ilocs(2) + njmpp - 1

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90

@@ -49 +49 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_offline = .FALSE.   !: offline flag
 #endif
+
+   !! * Array index permutations
+#  include "dom_oce_ftrans.h90"
+#  include "trc_oce_ftrans.h90"
 
    !! * Substitutions
@@ -249 +253 @@
       !!----------------------------------------------------------------------
       !
-      ! It is not necessary to compute anything bellow the following depth
+      ! It is not necessary to compute anything below the following depth
       zhext = prldex * ( LOG(10._wp) * zprec + LOG(pqsr_frc) )
       !
       ! Level of light extinction
+!FTRANS fsdepw :I :I :z
+!FTRANS tmask  :I :I :z
       pjl = jpkm1
       DO jk = jpkm1, 1, -1

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/vectopt_loop_substitute.h90

@@ -2 +2 @@
    !!                   ***  vectopt_loop_substitute  ***
    !!----------------------------------------------------------------------
-   !! ** purpose :   substitute the inner loop starting and inding indices 
+   !! ** purpose :   substitute the inner loop starting and ending indices 
    !!      to allow unrolling of do-loop using CPP macro.
    !!----------------------------------------------------------------------

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/wrk_nemo.F90

@@ -98 +98 @@
    CHARACTER(LEN=*), PARAMETER ::   cform_err2 = "(/,' ===>>> : E R R O R',     /,'         ===========',/)"       !:
    CHARACTER(LEN=*), PARAMETER ::   cform_war2 = "(/,' ===>>> : W A R N I N G', /,'         ===============',/)"   !:
+
+   !! * Array index permutations
+#  include "wrk_nemo_ftrans.h90"
 
    !!----------------------------------------------------------------------