From afe038d1e119968c942526c8a00dd789648fe81e Mon Sep 17 00:00:00 2001
From: cetlod <Christian.Ethe@ipsl.fr>
Date: Thu, 7 Jul 2022 16:23:46 +0200
Subject: [PATCH] First step of TOP halo cleanup

---
 src/TOP/trc.F90        | 30 +++++++++++++++---------------
 src/TOP/trcais.F90     |  8 ++++----
 src/TOP/trcbc.F90      |  6 +++---
 src/TOP/trcnam.F90     |  7 ++++++-
 src/TOP/trcopt.F90     |  2 +-
 src/TOP/trcstp.F90     | 26 +++++---------------------
 src/TOP/trcstp_rk3.F90 | 37 +++++++++----------------------------
 src/TOP/trcwri.F90     | 36 ++++++++++++++++++++----------------
 8 files changed, 63 insertions(+), 89 deletions(-)

diff --git a/src/TOP/trc.F90 b/src/TOP/trc.F90
index d5b85e041..10c40818b 100644
--- a/src/TOP/trc.F90
+++ b/src/TOP/trc.F90
@@ -158,26 +158,26 @@ CONTAINS
       !!-------------------------------------------------------------------
       ierr(:) = 0
       !
-      ALLOCATE( tr(jpi,jpj,jpk,jptra,jpt)                                             ,       &  
-         &      trc_i(jpi,jpj,jptra)  , trc_o(jpi,jpj,jptra)                          ,       &
-         &      gtru (jpi,jpj,jptra)  , gtrv (jpi,jpj,jptra)                          ,       &
-         &      gtrui(jpi,jpj,jptra)  , gtrvi(jpi,jpj,jptra)                          ,       &
-         &      trc_ice_ratio(jptra)  , trc_ice_prescr(jptra) , cn_trc_o(jptra)       ,       &
-         &      neln(jpi,jpj)         , heup(jpi,jpj)         , heup_01(jpi,jpj)      ,       &
-         &      etot(jpi,jpj,jpk)     , etot_ndcy(jpi,jpj,jpk)                        ,       &
-         &      sbc_trc_b(jpi,jpj,jptra), sbc_trc(jpi,jpj,jptra)                      ,       &  
-         &      cvol(jpi,jpj,jpk)     , trai(jptra)                                   ,       &
-         &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       &
-         &      ln_trc_ini(jptra)     ,                                                       &
-         &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       &
-         &      ln_trc_ais(jptra)     ,                                                       &
+      ALLOCATE( tr(jpi,jpj,jpk,jptra,jpt)                                         ,       &  
+         &      gtru (jpi,jpj,jptra) , gtrv (jpi,jpj,jptra)                       ,       &
+         &      gtrui(jpi,jpj,jptra) , gtrvi(jpi,jpj,jptra)                       ,       &
+         &      trc_i(A2D(0),jptra)  , trc_o(A2D(0),jptra)                        ,       &
+         &      trc_ice_ratio(jptra) , trc_ice_prescr(jptra) , cn_trc_o(jptra)    ,       &
+         &      neln(A2D(0))         , heup(A2D(0))         , heup_01(A2D(0))     ,       &
+         &      etot(A2D(0),jpk)     , etot_ndcy(A2D(0),jpk)                      ,       &
+         &      sbc_trc_b(A2D(0),jptra), sbc_trc(A2D(0),jptra)                    ,       &  
+         &      cvol(jpi,jpj,jpk)    , trai(jptra)                                ,       &
+         &      ctrcnm(jptra)        , ctrcln(jptra)         , ctrcun(jptra)      ,       &
+         &      ln_trc_ini(jptra)    ,                                                    &
+         &      ln_trc_sbc(jptra)    , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)  ,       &
+         &      ln_trc_ais(jptra)    ,                                                    &
          &      STAT = ierr(1)  )
       !
       IF( ln_bdy       )   ALLOCATE( trcdta_bdy(jptra, jp_bdy)  , STAT = ierr(2) )
       !
-      IF (jp_dia3d > 0 )   ALLOCATE( trc3d(jpi,jpj,jpk,jp_dia3d), STAT = ierr(3) )
+      IF (jp_dia3d > 0 )   ALLOCATE( trc3d(A2D(0),jpk,jp_dia3d), STAT = ierr(3) )
       !
-      IF (jp_dia2d > 0 )   ALLOCATE( trc2d(jpi,jpj,jp_dia2d)    , STAT = ierr(4) )
+      IF (jp_dia2d > 0 )   ALLOCATE( trc2d(A2D(0),jp_dia2d)    , STAT = ierr(4) )
       ! 
       trc_alloc = MAXVAL( ierr )
       IF( trc_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trc_alloc: failed to allocate arrays' )
diff --git a/src/TOP/trcais.F90 b/src/TOP/trcais.F90
index 5bd051387..010e9fce1 100644
--- a/src/TOP/trcais.F90
+++ b/src/TOP/trcais.F90
@@ -169,7 +169,7 @@ CONTAINS
             DO jn = 1, jptra
                IF( ln_trc_ais(jn) ) THEN
                   jl = n_trc_indais(jn)
-                  DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+                  DO_2D( 0, 0, 0, 0 )
                      zfact = 1. / e3t(ji,jj,1,Kmm)
                      ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + fwficb(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) * zfact
                   END_2D
@@ -181,7 +181,7 @@ CONTAINS
             DO jn = 1, jptra
                IF( ln_trc_ais(jn) ) THEN
                   jl = n_trc_indais(jn)
-                  DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+                  DO_2D( 0, 0, 0, 0 )
                      IF( ln_isfpar_mlt ) THEN
                         zcalv = fwfisf_par(ji,jj) * r1_rho0 / rhisf_tbl_par(ji,jj)
                         ikt = misfkt_par(ji,jj)
@@ -213,7 +213,7 @@ CONTAINS
             DO jn = 1, jptra
                IF( ln_trc_ais(jn) ) THEN
                   jl = n_trc_indais(jn)
-                  DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+                  DO_2D( 0, 0, 0, 0 )
                      DO jk = 1, icblev
                         zcalv  =  fwficb(ji,jj) * r1_rho0 
                         ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + rf_trafac(jl) * zcalv / gdepw(ji,jj,icblev+1,Kmm)
@@ -228,7 +228,7 @@ CONTAINS
             DO jn = 1, jptra
                IF( ln_trc_ais(jn) ) THEN
                   jl = n_trc_indais(jn)
-                  DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+                  DO_2D( 0, 0, 0, 0 )
                      IF( ln_isfpar_mlt ) THEN
                         zcalv = - fwfisf_par(ji,jj) * r1_rho0 / rhisf_tbl_par(ji,jj)
                         ikt = misfkt_par(ji,jj)
diff --git a/src/TOP/trcbc.F90 b/src/TOP/trcbc.F90
index fc829d8db..141ecf03c 100644
--- a/src/TOP/trcbc.F90
+++ b/src/TOP/trcbc.F90
@@ -414,7 +414,7 @@ CONTAINS
          !
          ! Remove river dilution for tracers with absent river load
          IF( ln_rnf_ctl .AND. .NOT.ln_trc_cbc(jn) ) THEN
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                DO jk = 1, nk_rnf(ji,jj)
 #if defined key_RK3
                   zrnf =  rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj)
@@ -432,7 +432,7 @@ CONTAINS
          IF( ln_trc_sbc(jn) ) THEN
             jl = n_trc_indsbc(jn)
             sf_trcsbc(jl)%fnow(:,:,1) = MAX( rtrn, sf_trcsbc(jl)%fnow(:,:,1) ) ! avoid nedgative value due to interpolation
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                zfact = 1. / ( e3t(ji,jj,1,Kmm) * rn_sbc_time )
                ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) * zfact
             END_2D
@@ -443,7 +443,7 @@ CONTAINS
             IF( l_offline )   rn_rfact = 1._wp
             jl = n_trc_indcbc(jn)
             sf_trccbc(jl)%fnow(:,:,1) = MAX( rtrn, sf_trccbc(jl)%fnow(:,:,1) ) ! avoid nedgative value due to interpolation
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                DO jk = 1, nk_rnf(ji,jj)
                   zfact = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) 
                   ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zfact
diff --git a/src/TOP/trcnam.F90 b/src/TOP/trcnam.F90
index cacbe6171..80a181e05 100644
--- a/src/TOP/trcnam.F90
+++ b/src/TOP/trcnam.F90
@@ -254,7 +254,12 @@ CONTAINS
          WRITE(numout,*) '   Namelist : namtrc_dcy                    '
          WRITE(numout,*) '      Diurnal cycle for TOP ln_trcdc2dm    = ', ln_trcdc2dm
       ENDIF
-
+      !      ! Define logical parameter ton control dirunal cycle in TOP
+      l_trcdm2dc = ( ln_trcdc2dm .AND. .NOT. ln_dm2dc  ) 
+      !
+      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
+         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
+      !
    END SUBROUTINE trc_nam_dcy
 
    SUBROUTINE trc_nam_trd
diff --git a/src/TOP/trcopt.F90 b/src/TOP/trcopt.F90
index 7620a4c00..d07b4500d 100644
--- a/src/TOP/trcopt.F90
+++ b/src/TOP/trcopt.F90
@@ -348,7 +348,7 @@ CONTAINS
       !!                     ***  ROUTINE trc_opt_alloc  ***
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
+      ALLOCATE( ekb(jpi,jpj,jpk),ekr(jpi,jpj,jpk),  &
                 ekg(jpi,jpj,jpk),zeps(jpi,jpj,jpk),  STAT= trc_opt_alloc  ) 
       !
       IF( trc_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_opt_alloc : failed to allocate arrays.' )
diff --git a/src/TOP/trcstp.F90 b/src/TOP/trcstp.F90
index 9d7e7574b..065d7f73d 100644
--- a/src/TOP/trcstp.F90
+++ b/src/TOP/trcstp.F90
@@ -37,6 +37,8 @@ MODULE trcstp
    REAL(wp) ::   rsecfst, rseclast       ! ???
    REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr   ! save qsr during TOP time-step
 
+      !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -74,17 +76,13 @@ CONTAINS
       ll_trcstat  = ( sn_cfctl%l_trcstat ) .AND. &
      &              ( ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) )
 
-      IF( kt == nittrc000 )                      CALL trc_stp_ctl   ! control 
       IF( kt == nittrc000 .AND. lk_trdmxl_trc )  CALL trd_mxl_trc_init    ! trends: Mixed-layer
       !
       IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
          DO jk = 1, jpk
             cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
-         IF ( ll_trcstat .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend )   &
-            & .OR. iom_use( "pno3tot" ) .OR. iom_use( "ppo4tot" ) .OR. iom_use( "psiltot" )   &
-            & .OR. iom_use( "palktot" ) .OR. iom_use( "pfertot" ) )                           &
-            &     areatot = glob_sum( 'trcstp', cvol(:,:,:) )
+         IF ( ll_trcstat .OR. kt == nitrst )  areatot = glob_sum( 'trcstp', cvol(:,:,:) )
       ENDIF
       !
       IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
@@ -141,20 +139,6 @@ CONTAINS
    END SUBROUTINE trc_stp
 
 
-   SUBROUTINE trc_stp_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_stp_ctl  ***
-      !!----------------------------------------------------------------------
-      !
-      ! Define logical parameter ton control dirunal cycle in TOP
-      l_trcdm2dc = ( ln_trcdc2dm .AND. .NOT. ln_dm2dc  ) 
-      !
-      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
-         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
-      !
-   END SUBROUTINE trc_stp_ctl
-
-
    SUBROUTINE trc_mean_qsr( kt )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE trc_mean_qsr  ***
@@ -188,7 +172,7 @@ CONTAINS
             WRITE(numout,*) 
          ENDIF
          !
-         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
+         ALLOCATE( qsr_arr(A2D(0),nb_rec_per_day ) )
          !
          !                                            !* Restart: read in restart file
          IF( ln_rsttr .AND. nn_rsttr /= 0 .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  &
@@ -239,7 +223,7 @@ CONTAINS
              qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
           ENDDO
           qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
-          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
+          qsr_mean(:,:) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
       ENDIF
       !
       IF( lrst_trc ) THEN    !* Write the mean of qsr in restart file 
diff --git a/src/TOP/trcstp_rk3.F90 b/src/TOP/trcstp_rk3.F90
index 450d4b5ed..cdf89df57 100644
--- a/src/TOP/trcstp_rk3.F90
+++ b/src/TOP/trcstp_rk3.F90
@@ -41,6 +41,8 @@ MODULE trcstp_rk3
    REAL(wp) ::   rsecfst, rseclast       ! ???
    REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr   ! save qsr during TOP time-step
 
+      !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -71,15 +73,14 @@ CONTAINS
       l_trcstat  = ( sn_cfctl%l_trcstat ) .AND. &
            &       ( ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) )
       !
-      IF( kt == nittrc000 )                      CALL trc_stp_ctl   ! control 
+      IF( kt == nittrc000 )                      CALL trc_stpsctl   ! control 
       IF( kt == nittrc000 .AND. lk_trdmxl_trc )  CALL trd_mxl_trc_init    ! trends: Mixed-layer
       !
       IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
          DO jk = 1, jpk
             cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
-         IF( l_trcstat .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend ) )   &
-            &     areatot = glob_sum( 'trcstp', cvol(:,:,:) )
+         IF( l_trcstat .OR. kt == nitrst ) areatot = glob_sum( 'trcstp', cvol(:,:,:) )
       ENDIF
       !
       IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
@@ -146,22 +147,6 @@ CONTAINS
    END SUBROUTINE trc_stp_end
 
 
-   SUBROUTINE trc_stp_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_stp_ctl  ***
-      !! ** Purpose :        Control  + ocean volume
-      !!----------------------------------------------------------------------
-      !
-      ! Define logical parameter ton control dirunal cycle in TOP
-      l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 .AND. ncpl_qsr_freq /= 0 )
-      l_trcdm2dc = l_trcdm2dc .AND. .NOT. l_offline
-      !
-      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
-         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
-      !
-   END SUBROUTINE trc_stp_ctl
-
-
    SUBROUTINE trc_mean_qsr( kt )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE trc_mean_qsr  ***
@@ -185,13 +170,9 @@ CONTAINS
       IF( ln_timing )   CALL timing_start('trc_mean_qsr')
       !
       IF( kt == nittrc000 ) THEN
-         IF( ln_cpl )  THEN  
-            rdt_sampl = rday / ncpl_qsr_freq
-            nb_rec_per_day = ncpl_qsr_freq
-         ELSE  
-            rdt_sampl = MAX( 3600., rn_Dt )
-            nb_rec_per_day = INT( rday / rdt_sampl )
-         ENDIF
+         !
+         rdt_sampl = REAL( ncpl_qsr_freq )
+         nb_rec_per_day = INT( rday / ncpl_qsr_freq )
          !
          IF(lwp) THEN
             WRITE(numout,*) 
@@ -199,7 +180,7 @@ CONTAINS
             WRITE(numout,*) 
          ENDIF
          !
-         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
+         ALLOCATE( qsr_arr(A2D(0),nb_rec_per_day ) )
          !
          !                                            !* Restart: read in restart file
          IF( ln_rsttr .AND. nn_rsttr /= 0 .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  &
@@ -250,7 +231,7 @@ CONTAINS
              qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
           END DO
           qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
-          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
+          qsr_mean(:,:) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
       ENDIF
       !
       IF( lrst_trc ) THEN    !* Write the mean of qsr in restart file 
diff --git a/src/TOP/trcwri.F90 b/src/TOP/trcwri.F90
index 39e3123fe..8849c2256 100644
--- a/src/TOP/trcwri.F90
+++ b/src/TOP/trcwri.F90
@@ -42,12 +42,12 @@ CONTAINS
       INTEGER, INTENT( in )     :: kt
       INTEGER, INTENT( in )     :: Kmm  ! time level indices
       !
-      INTEGER                   :: jk, jn
+      INTEGER                   :: ji,jj,jk,jn
       CHARACTER (len=20)        :: cltra
       CHARACTER (len=40)        :: clhstnam
       INTEGER ::   inum = 11            ! temporary logical unit
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
-      !!---------------------------------------------------------------------
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3d   ! 3D workspace
+      !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_wri')
       !
@@ -58,6 +58,8 @@ CONTAINS
            WRITE(inum,*) clhstnam
            CLOSE(inum)
          ENDIF
+         
+         ALLOCATE( z3d(A2D(0),jpk) )  ;   z3d(A2D(0),:) = 0._wp
 
          ! Output of initial vertical scale factor
          CALL iom_put( "e3t_0", e3t_0(:,:,:) )
@@ -66,25 +68,27 @@ CONTAINS
          !
          IF( .NOT.ln_linssh )  CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
          !
-         IF ( iom_use("e3t") ) THEN  ! time-varying e3t
-            DO jk = 1, jpk
-               z3d(:,:,jk) =  e3t(:,:,jk,Kmm)
-            END DO
-            CALL iom_put( "e3t", z3d(:,:,:) )
+         ! --- vertical scale factors --- !
+         IF( iom_use("e3t") ) THEN  ! time-varying e3t
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               z3d(ji,jj,jk) =  e3t(ji,jj,jk,Kmm)
+            END_3D
+            CALL iom_put( "e3t", z3d )
          ENDIF
          IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
-            DO jk = 1, jpk
-               z3d(:,:,jk) =  e3u(:,:,jk,Kmm)
-            END DO
-            CALL iom_put( "e3u", z3d(:,:,:) )
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               z3d(ji,jj,jk) =  e3u(ji,jj,jk,Kmm)
+            END_3D 
+            CALL iom_put( "e3u" , z3d )
          ENDIF
          IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
-            DO jk = 1, jpk
-               z3d(:,:,jk) =  e3v(:,:,jk,Kmm)
-            END DO
-            CALL iom_put( "e3v", z3d(:,:,:) )
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               z3d(ji,jj,jk) =  e3v(ji,jj,jk,Kmm)
+            END_3D
+            CALL iom_put( "e3v" , z3d )
          ENDIF
          !
+         DEALLOCATE( z3d )
       ENDIF
       !
       ! write the tracer concentrations in the file
-- 
GitLab