Changeset 9019 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA

Timestamp:

2017-12-13T15:58:53+01:00 (6 years ago)

Author:

timgraham

Message:

Merge of dev_CNRS_2017 into branch

Location:

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA

Files:

• eosbn2.F90 (modified) (20 diffs)
• traadv.F90 (modified) (16 diffs)
• traadv_cen.F90 (modified) (8 diffs)
• traadv_fct.F90 (modified) (11 diffs)
• traadv_mle.F90 (modified) (3 diffs)
• traadv_mus.F90 (modified) (3 diffs)
• traadv_qck.F90 (modified) (13 diffs)
• traadv_ubs.F90 (modified) (5 diffs)
• trabbc.F90 (modified) (7 diffs)
• trabbl.F90 (modified) (16 diffs)
• tradmp.F90 (modified) (4 diffs)
• traldf.F90 (modified) (6 diffs)
• traldf_iso.F90 (modified) (3 diffs)
• traldf_lap_blp.F90 (modified) (6 diffs)
• traldf_triad.F90 (modified) (3 diffs)
• tranpc.F90 (modified) (6 diffs)
• tranxt.F90 (modified) (12 diffs)
• traqsr.F90 (modified) (11 diffs)
• trasbc.F90 (modified) (5 diffs)
• trazdf.F90 (modified) (8 diffs)
• trazdf_exp.F90 (deleted)
• trazdf_imp.F90 (deleted)
• zpshde.F90 (modified) (5 diffs)

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

@@ -46 +46 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE lbclnk         ! ocean lateral boundary conditions
    USE timing         ! Timing
@@ -231 +230 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('eos-insitu')
+      IF( ln_timing )   CALL timing_start('eos-insitu')
       !
       SELECT CASE( neos )
@@ -298 +297 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', ovlap=1, kdim=jpk )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('eos-insitu')
+      IF( ln_timing )   CALL timing_stop('eos-insitu')
       !
    END SUBROUTINE eos_insitu
@@ -329 +328 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('eos-pot')
+      IF( ln_timing )   CALL timing_start('eos-pot')
       !
       SELECT CASE ( neos )
@@ -465 +464 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('eos-pot')
+      IF( ln_timing )   CALL timing_stop('eos-pot')
       !
    END SUBROUTINE eos_insitu_pot
@@ -491 +490 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('eos2d')
+      IF( ln_timing )   CALL timing_start('eos2d')
       !
       prd(:,:) = 0._wp
@@ -560 +559 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('eos2d')
+      IF( ln_timing )   CALL timing_stop('eos2d')
       !
    END SUBROUTINE eos_insitu_2d
@@ -583 +582 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('rab_3d')
+      IF( ln_timing )   CALL timing_start('rab_3d')
       !
       SELECT CASE ( neos )
@@ -674 +673 @@
          &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('rab_3d')
+      IF( ln_timing )   CALL timing_stop('rab_3d')
       !
    END SUBROUTINE rab_3d
@@ -696 +695 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
+      IF( ln_timing )   CALL timing_start('rab_2d')
       !
       pab(:,:,:) = 0._wp
@@ -791 +790 @@
          &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
+      IF( ln_timing )   CALL timing_stop('rab_2d')
       !
    END SUBROUTINE rab_2d
@@ -812 +811 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('rab_2d')
+      IF( ln_timing )   CALL timing_start('rab_2d')
       !
       pab(:) = 0._wp
@@ -888 +887 @@
       END SELECT
       !
-      IF( nn_timing == 1 )   CALL timing_stop('rab_2d')
+      IF( ln_timing )   CALL timing_stop('rab_2d')
       !
    END SUBROUTINE rab_0d
@@ -915 +914 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bn2')
+      IF( ln_timing )   CALL timing_start('bn2')
       !
       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
@@ -928 +927 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
             END DO
          END DO
@@ -935 +934 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('bn2')
+      IF( ln_timing )   CALL timing_stop('bn2')
       !
    END SUBROUTINE bn2
@@ -963 +962 @@
       !!----------------------------------------------------------------------
       !
-      IF ( nn_timing == 1 )   CALL timing_start('eos_pt_from_ct')
+      IF( ln_timing )   CALL timing_start('eos_pt_from_ct')
       !
       zdeltaS = 5._wp
@@ -994 +993 @@
       END DO
       !
-      IF( nn_timing == 1 )   CALL timing_stop('eos_pt_from_ct')
+      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct')
       !
    END FUNCTION eos_pt_from_ct
@@ -1128 +1127 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('eos_pen')
+      IF( ln_timing )   CALL timing_start('eos_pen')
       !
       SELECT CASE ( neos )
@@ -1222 +1221 @@
       END SELECT
       !
-      IF( nn_timing == 1 )   CALL timing_stop('eos_pen')
+      IF( ln_timing )   CALL timing_stop('eos_pen')
       !
    END SUBROUTINE eos_pen

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

@@ -14 +14 @@
    !!----------------------------------------------------------------------
    !!   tra_adv       : compute ocean tracer advection trend
-   !!   tra_adv_ctl   : control the different options of advection scheme
+   !!   tra_adv_init  : control the different options of advection scheme
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
    USE domvvl         ! variable vertical scale factors
+   USE sbcwave        ! wave module
+   USE sbc_oce        ! surface boundary condition: ocean
    USE traadv_cen     ! centered scheme           (tra_adv_cen  routine)
    USE traadv_fct     ! FCT      scheme           (tra_adv_fct  routine)
@@ -27 +29 @@
    USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
    USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces
-   USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trd_oce        ! trends: ocean variables
+   USE trdtra         ! trends manager: tracers 
+   USE diaptr         ! Poleward heat transport 
    !
    USE in_out_manager ! I/O manager
@@ -34 +37 @@
    USE prtctl         ! Print control
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE sbcwave        ! wave module
-   USE sbc_oce        ! surface boundary condition: ocean
-   USE diaptr         ! Poleward heat transport 
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_adv        ! routine called by step module
-   PUBLIC   tra_adv_init   ! routine called by opa module
+   PUBLIC   tra_adv        ! called by step.F90
+   PUBLIC   tra_adv_init   ! called by nemogcm.F90
 
    !                            !!* Namelist namtra_adv *
+   LOGICAL ::   ln_traadv_NONE   ! no advection on T and S
    LOGICAL ::   ln_traadv_cen    ! centered scheme flag
    INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme
    LOGICAL ::   ln_traadv_fct    ! FCT scheme flag
    INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme
-   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping
    LOGICAL ::   ln_traadv_mus    ! MUSCL scheme flag
    LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths
@@ -58 +57 @@
    LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag
 
-   INTEGER ::              nadv             ! choice of the type of advection scheme
-   !
-   !                                        ! associated indices:
+   INTEGER ::   nadv             ! choice of the type of advection scheme
+   !                             ! associated indices:
    INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection
    INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme
    INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
-   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping
-   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme
-   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme
-   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme
+   INTEGER, PARAMETER ::   np_MUS     = 3   ! MUSCL scheme
+   INTEGER, PARAMETER ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme
+   INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
    
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -86 +83 @@
       !! ** Method  : - Update (ua,va) with the advection term following nadv
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
       INTEGER ::   jk   ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zun, zvn, zwn )
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        :: zun, zvn, zwn   ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_adv')
       !
       !                                          ! set time step
-      zun(:,:,:) = 0.0
-      zvn(:,:,:) = 0.0
-      zwn(:,:,:) = 0.0
-      !    
-      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
-         r2dt = rdt                                 ! = rdt (restarting with Euler time stepping)
-      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
-         r2dt = 2._wp * rdt                         ! = 2 rdt (leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
       !
       !                                         !==  effective transport  ==!
+      zun(:,:,jpk) = 0._wp
+      zvn(:,:,jpk) = 0._wp
+      zwn(:,:,jpk) = 0._wp
       IF( ln_wave .AND. ln_sdw )  THEN
          DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
@@ -146 +138 @@
       !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
@@ -153 +145 @@
       SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==!
       !
-      CASE ( np_CEN )                                    ! Centered scheme : 2nd / 4th order
+      CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order
          CALL tra_adv_cen    ( kt, nit000, 'TRA',         zun, zvn, zwn     , tsn, tsa, jpts, nn_cen_h, nn_cen_v )
-      CASE ( np_FCT )                                    ! FCT scheme      : 2nd / 4th order
+      CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
          CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v )
-      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting
-         CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_fct_zts )
-      CASE ( np_MUS )                                    ! MUSCL
+      CASE ( np_MUS )                                 ! MUSCL
          CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts        , ln_mus_ups ) 
-      CASE ( np_UBS )                                    ! UBS
+      CASE ( np_UBS )                                 ! UBS
          CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_ubs_v   )
-      CASE ( np_QCK )                                    ! QUICKEST
+      CASE ( np_QCK )                                 ! QUICKEST
          CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts                     )
       !
@@ -175 +165 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+         DEALLOCATE( ztrdt, ztrds )
       ENDIF
       !                                              ! print mean trends (used for debugging)
@@ -181 +171 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv' )
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zun, zvn, zwn )
-      !                                          
+      IF( ln_timing )   CALL timing_stop( 'tra_adv' )
+      !
    END SUBROUTINE tra_adv
 
@@ -197 +185 @@
       INTEGER ::   ioptio, ios   ! Local integers
       !
-      NAMELIST/namtra_adv/ ln_traadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN
-         &                 ln_traadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT
-         &                 ln_traadv_mus,                     ln_mus_ups,   &   ! MUSCL
-         &                 ln_traadv_ubs,           nn_ubs_v,               &   ! UBS
-         &                 ln_traadv_qck                                        ! QCK
+      NAMELIST/namtra_adv/ ln_traadv_NONE,                       &   ! No advection
+         &                 ln_traadv_cen , nn_cen_h, nn_cen_v,   &   ! CEN
+         &                 ln_traadv_fct , nn_fct_h, nn_fct_v,   &   ! FCT
+         &                 ln_traadv_mus , ln_mus_ups,           &   ! MUSCL
+         &                 ln_traadv_ubs ,           nn_ubs_v,   &   ! UBS
+         &                 ln_traadv_qck                             ! QCK
       !!----------------------------------------------------------------------
       !
@@ -217 +206 @@
          WRITE(numout,*)
          WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_adv : chose a advection scheme for tracers'
+         WRITE(numout,*) '      No advection on T & S                     ln_traadv_NONE= ', ln_traadv_NONE
          WRITE(numout,*) '      centered scheme                           ln_traadv_cen = ', ln_traadv_cen
          WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h
@@ -225 +215 @@
          WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h
          WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v
-         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts
          WRITE(numout,*) '      MUSCL scheme                              ln_traadv_mus = ', ln_traadv_mus
          WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups
@@ -233 +222 @@
       ENDIF
       !
-      ioptio = 0                       !==  Parameter control  ==!
-      IF( ln_traadv_cen )   ioptio = ioptio + 1
-      IF( ln_traadv_fct )   ioptio = ioptio + 1
-      IF( ln_traadv_mus )   ioptio = ioptio + 1
-      IF( ln_traadv_ubs )   ioptio = ioptio + 1
-      IF( ln_traadv_qck )   ioptio = ioptio + 1
-      !
-      IF( ioptio == 0 ) THEN
-         nadv = np_NO_adv
-         CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' )
-      ENDIF
-      IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' )
+      !                                !==  Parameter control & set nadv ==!
+      ioptio = 0                       
+      IF( ln_traadv_NONE ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_NO_adv   ;   ENDIF
+      IF( ln_traadv_cen  ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_CEN      ;   ENDIF
+      IF( ln_traadv_fct  ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_FCT      ;   ENDIF
+      IF( ln_traadv_mus  ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_MUS      ;   ENDIF
+      IF( ln_traadv_ubs  ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_UBS      ;   ENDIF
+      IF( ln_traadv_qck  ) THEN   ;   ioptio = ioptio + 1   ;   nadv = np_QCK      ;   ENDIF
+      !
+      IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection option in namelist namtra_adv' )
       !
       IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &          ! Centered
@@ -254 +241 @@
         CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' )
       ENDIF
-      IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) THEN
-         IF( nn_fct_h == 4 ) THEN
-            nn_fct_h = 2
-            CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' )
-         ENDIF
-         IF( .NOT.ln_linssh ) THEN
-            CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' )
-         ENDIF
-         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'tra_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' )
-      ENDIF
       IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS
         CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' )
@@ -275 +252 @@
       ENDIF
       !
-      !                                !==  used advection scheme  ==!  
-      !                                      ! set nadv
-      IF( ln_traadv_cen                      )   nadv = np_CEN
-      IF( ln_traadv_fct                      )   nadv = np_FCT
-      IF( ln_traadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts
-      IF( ln_traadv_mus                      )   nadv = np_MUS
-      IF( ln_traadv_ubs                      )   nadv = np_UBS
-      IF( ln_traadv_qck                      )   nadv = np_QCK
-      !
-      IF(lwp) THEN                           ! Print the choice
+      !                                !==  Print the choice  ==!  
+      IF(lwp) THEN
          WRITE(numout,*)
          SELECT CASE ( nadv )
@@ -292 +261 @@
          CASE( np_FCT     )   ;   WRITE(numout,*) '      ===>>   FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
             &                                                                      ' Vertical   order: ', nn_fct_v
-         CASE( np_FCT_zts )   ;   WRITE(numout,*) '      ===>>   use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
          CASE( np_MUS     )   ;   WRITE(numout,*) '      ===>>   MUSCL    scheme is used'
          CASE( np_UBS     )   ;   WRITE(numout,*) '      ===>>   UBS      scheme is used'

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90

@@ -11 +11 @@
    !!                   NB: on the vertical it is actually a 4th order COMPACT scheme which is used
    !!----------------------------------------------------------------------
-   USE oce      , ONLY: tsn ! now ocean temperature and salinity
    USE dom_oce        ! ocean space and time domain
    USE eosbn2         ! equation of state
@@ -24 +23 @@
    USE trc_oce        ! share passive tracers/Ocean variables
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -30 +28 @@
    PRIVATE
 
-   PUBLIC   tra_adv_cen       ! routine called by step.F90
+   PUBLIC   tra_adv_cen   ! called by traadv.F90
    
    REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
 
-   LOGICAL :: l_trd   ! flag to compute trends
-   LOGICAL :: l_ptr   ! flag to compute poleward transport
-   LOGICAL :: l_hst   ! flag to compute heat/salt transport
+   LOGICAL ::   l_trd   ! flag to compute trends
+   LOGICAL ::   l_ptr   ! flag to compute poleward transport
+   LOGICAL ::   l_hst   ! flag to compute heat/salt transport
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id$
+   !! NEMO/OPA 4.0, NEMO Consortium (2017)
+   !! $Id:$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -48 +46 @@
 
    SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn,     &
-      &                                             ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 
+      &                                        ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_cen  ***
@@ -80 +78 @@
       REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
       REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zwy, zwz, ztu, ztv, ztw
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, ztu, ztv, ztw
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
+      IF( ln_timing )   CALL timing_start('tra_adv_cen')
       !
       IF( kt == kit000 )  THEN
@@ -92 +88 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
       ENDIF
-      !
+      !                          ! set local switches
       l_trd = .FALSE.
       l_hst = .FALSE.
@@ -130 +126 @@
                END DO
             END DO
-            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
+            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond.
             !
             DO jk = 1, jpkm1                       ! Horizontal advective fluxes
@@ -203 +199 @@
          END IF
          !                                 ! "Poleward" heat and salt transports 
-         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !                                 !  heat and salt transport
-         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
+         IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen')
+      IF( ln_timing )   CALL timing_stop('tra_adv_cen')
       !
    END SUBROUTINE tra_adv_cen

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

@@ -9 +9 @@
    !!----------------------------------------------------------------------
    !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
-   !!  tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme 
    !!                   with sub-time-stepping in the vertical direction
    !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm 
@@ -21 +20 @@
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
-   USE phycst, ONLY: rau0_rcp
+   USE phycst  , ONLY : rau0_rcp
    !
    USE in_out_manager ! I/O manager
-   USE iom
+   USE iom            ! 
    USE lib_mpp        ! MPP library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -34 +32 @@
    PRIVATE
 
-   PUBLIC   tra_adv_fct        ! routine called by traadv.F90
-   PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90
-   PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90
+   PUBLIC   tra_adv_fct        ! called by traadv.F90
+   PUBLIC   interp_4th_cpt     ! called by traadv_cen.F90
 
    LOGICAL  ::   l_trd   ! flag to compute trends
@@ -50 +47 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -70 +67 @@
       !!
       !! ** Action : - update pta  with the now advective tracer trends
-      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T)
       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
       !!----------------------------------------------------------------------
@@ -88 +85 @@
       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz, zptry
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_adv_fct')
       !
       IF( kt == kit000 )  THEN
@@ -103 +97 @@
       ENDIF
       !
-      l_trd = .FALSE.
+      l_trd = .FALSE.            ! set local switches
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE.
-      IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE. 
-      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
-         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
+      IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
+      IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
+      IF(   cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  &
+         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
       !
       IF( l_trd .OR. l_hst )  THEN
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
          ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
       IF( l_ptr ) THEN  
-         CALL wrk_alloc( jpi, jpj, jpk, zptry )
+         ALLOCATE( zptry(jpi,jpj,jpk) )
          zptry(:,:,:) = 0._wp
       ENDIF
@@ -184 +178 @@
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
+         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
@@ -308 +302 @@
          END DO
          !
-         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
-            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
-            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
-            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes 
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  !
+            !
+            IF( l_trd ) THEN              ! trend diagnostics
+               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
+               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
+               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
+            ENDIF
+            !                             ! heat/salt transport
+            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
+            !
+            DEALLOCATE( ztrdx, ztrdy, ztrdz )
          ENDIF
-            !
-         IF( l_trd ) THEN 
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
-            !
-         END IF
-         !                                !  heat/salt transport
-         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
-
-         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr ) THEN  
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+         IF( l_ptr ) THEN              ! "Poleward" transports
+            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes
             CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
+            DEALLOCATE( zptry )
          ENDIF
          !
       END DO                     ! end of tracer loop
       !
-                              CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
-      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
+      IF( ln_timing )   CALL timing_stop('tra_adv_fct')
       !
    END SUBROUTINE tra_adv_fct
-
-
-   SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                                  ptb, ptn, pta, kjpt, kn_fct_zts )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_adv_fct_zts  ***
-      !! 
-      !! **  Purpose :   Compute the now trend due to total advection of 
-      !!       tracers and add it to the general trend of tracer equations
-      !!
-      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
-      !!       corrected flux (monotonic correction). This version use sub-
-      !!       timestepping for the vertical advection which increases stability
-      !!       when vertical metrics are small.
-      !!       note: - this advection scheme needs a leap-frog time scheme
-      !!
-      !! ** Action : - update (pta) with the now advective tracer trends
-      !!             - save the trends 
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
-      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps
-      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
-      !
-      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
-      REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep
-      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices  
-      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
-      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
-      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
-      REAL(wp) ::   ztra            ! local scalar
-      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
-      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zwx_sav , zwy_sav
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts')
-      !
-      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
-      CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
-      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-      !
-      l_trd = .FALSE.
-      l_hst = .FALSE.
-      l_ptr = .FALSE.
-      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
-      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
-      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
-         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
-      !
-      IF( l_trd .OR. l_hst )  THEN
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
-         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
-      ENDIF
-      !
-      IF( l_ptr ) THEN  
-         CALL wrk_alloc( jpi, jpj,jpk, zptry )
-         zptry(:,:,:) = 0._wp
-      ENDIF
-      zwi(:,:,:) = 0._wp
-      z_rzts = 1._wp / REAL( kn_fct_zts, wp )
-      zr_p2dt = 1._wp / p2dt
-      !
-      ! surface & Bottom value : flux set to zero for all tracers
-      zwz(:,:, 1 ) = 0._wp
-      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
-      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
-      !
-      !                                                          ! ===========
-      DO jn = 1, kjpt                                            ! tracer loop
-         !                                                       ! ===========
-         !
-         ! Upstream advection with initial mass fluxes & intermediate update
-         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  ! upstream scheme
-                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
-                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
-                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
-                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
-                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
-                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
-               END DO
-            END DO
-         END DO
-         !                       ! upstream tracer flux in the k direction
-         DO jk = 2, jpkm1              ! Interior value
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
-                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
-                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask)
-            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
-                  END DO
-               END DO   
-            ELSE                             ! no cavities, surface value
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
-            ENDIF
-         ENDIF
-         !
-         DO jk = 1, jpkm1         ! total advective trend
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  !                             ! total intermediate advective trends
-                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj)
-                  !                             ! update and guess with monotonic sheme
-                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !                           
-         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
-         !                
-         IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes)
-            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
-         END IF
-         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)
-
-         ! 3. anti-diffusive flux : high order minus low order
-         ! ---------------------------------------------------
-
-         DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes
-            !
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
-                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
-                  !
-                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
-                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
-               END DO
-            END DO
-            !
-            DO jj = 2, jpjm1                    ! partial horizontal divergence
-               DO ji = fs_2, fs_jpim1
-                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
-                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
-               END DO
-            END DO
-            !
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
-                  zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
-               END DO
-            END DO
-         END DO
-         !
-         !                                !* vertical anti-diffusive flux
-         zwz_sav(:,:,:)   = zwz(:,:,:)
-         ztrs   (:,:,:,1) = ptb(:,:,:,jn)
-         ztrs   (:,:,1,2) = ptb(:,:,1,jn)
-         ztrs   (:,:,1,3) = ptb(:,:,1,jn)
-         zwzts  (:,:,:)   = 0._wp
-         !
-         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop
-            !
-            IF( jl == 1 ) THEN                        ! Euler forward to kick things off
-               jtb = 1   ;   jtn = 1   ;   jta = 2
-               zts(:) = p2dt * z_rzts
-               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux
-               !                                            ! starting at jl =1 if kn_fct_zts is odd; 
-               !                                            ! starting at jl =2 otherwise
-            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step
-               jtb = 1   ;   jtn = 2   ;   jta = 3
-               zts(:) = 2._wp * p2dt * z_rzts
-            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps
-               jtb = MOD(jtb,3) + 1
-               jtn = MOD(jtn,3) + 1
-               jta = MOD(jta,3) + 1
-            ENDIF
-            DO jk = 2, jpkm1                          ! interior value
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk)
-                     IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux
-                  END DO
-               END DO
-            END DO
-            IF( ln_linssh ) THEN                    ! top value (only in linear free surface case)
-               IF( ln_isfcav ) THEN                      ! ice-shelf cavities
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
-                     END DO
-                  END DO   
-               ELSE                                      ! no ocean cavities
-                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
-               ENDIF
-            ENDIF
-            !
-            jtaken = MOD( jtaken + 1 , 2 )
-            !
-            DO jk = 2, jpkm1                             ! total advective trends
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 &
-                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
-                        &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-            !
-         END DO
-
-         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
-         CALL lbc_lnk( zwz, 'W',  1. )
-
-         ! 4. monotonicity algorithm
-         ! -------------------------
-         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
-
-
-         ! 5. final trend with corrected fluxes
-         ! ------------------------------------
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       &
-                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-        !
-         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
-            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
-            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
-            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
-         ENDIF
-            !
-         IF( l_trd ) THEN 
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
-            !
-         END IF
-         !                                             ! heat/salt transport
-         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
-
-         !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr ) THEN  
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
-            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
-         ENDIF
-         !
-      END DO
-      !
-                              CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
-                              CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
-                              CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
-      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
-      !
-   END SUBROUTINE tra_adv_fct_zts
 
 
@@ -653 +351 @@
       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('nonosc')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('nonosc')
       !
       zbig  = 1.e+40_wp
@@ -734 +430 @@
       CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('nonosc')
+      IF( ln_timing )   CALL timing_stop('nonosc')
       !
    END SUBROUTINE nonosc

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90

@@ -15 +15 @@
    USE phycst         ! physical constant
    USE zdfmxl         ! mixed layer depth
+   !
    USE lbclnk         ! lateral boundary condition / mpp link
    USE in_out_manager ! I/O manager
    USE iom            ! IOM library
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
 
@@ -86 +86 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ikmax        ! temporary integer
-      REAL(wp) ::   zcuw, zmuw   ! local scalar
-      REAL(wp) ::   zcvw, zmvw   !   -      -
-      REAL(wp) ::   zc                                     !   -      -
-      !
-      INTEGER  ::   ii, ij, ik              ! local integers
-      INTEGER, DIMENSION(3) ::   ilocu      !
-      INTEGER, DIMENSION(2) ::   ilocs      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpsi_uw, zpsi_vw
-      INTEGER, POINTER, DIMENSION(:,:) :: inml_mle
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mle')
-      CALL wrk_alloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
-      CALL wrk_alloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
-      CALL wrk_alloc( jpi, jpj, inml_mle)
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
+      REAL(wp) ::   zcuw, zmuw, zc      ! local scalar
+      REAL(wp) ::   zcvw, zmvw          !   -      -
+      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
+      REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_adv_mle')
       !
       !                                      !==  MLD used for MLE  ==!
@@ -256 +248 @@
          CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
       ENDIF
-      CALL wrk_dealloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
-      CALL wrk_dealloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
-      CALL wrk_dealloc( jpi, jpj, inml_mle)
-
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mle')
+      !
+      IF( ln_timing )   CALL timing_stop('tra_adv_mle')
       !
    END SUBROUTINE tra_adv_mle

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

@@ -26 +26 @@
 
    !
-   USE iom
-   USE wrk_nemo       ! Memory Allocation
+   USE iom            ! XIOS library
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -85 +84 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
       !
-      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices
-      INTEGER  ::   ierr                 ! local integer
-      REAL(wp) ::   zu, z0u, zzwx, zw    ! local scalars
-      REAL(wp) ::   zv, z0v, zzwy, z0w   !   -      -
-      REAL(wp) ::   zalpha               !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      - 
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::   ierr             ! local integer
+      REAL(wp) ::   zu, z0u, zzwx, zw , zalpha   ! local scalars
+      REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      - 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mus')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
+      IF( ln_timing )   CALL timing_start('tra_adv_mus')
       !
       IF( kt == kit000 )  THEN
@@ -279 +275 @@
       END DO                     ! end of tracer loop
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mus')
+      IF( ln_timing )   CALL timing_stop('tra_adv_mus')
       !
    END SUBROUTINE tra_adv_mus

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

@@ -25 +25 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE in_out_manager  ! I/O manager
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -43 +42 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -96 +95 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
+      IF( ln_timing )   CALL timing_start('tra_adv_qck')
       !
       IF( kt == kit000 )  THEN
@@ -118 +117 @@
       CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
+      IF( ln_timing )   CALL timing_stop('tra_adv_qck')
       !
    END SUBROUTINE tra_adv_qck
@@ -138 +137 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
       !----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
@@ -230 +228 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
          !
       END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
       !
    END SUBROUTINE tra_adv_qck_i
@@ -252 +248 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
       !!
-      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
       REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
       !----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
       !
       !                                                          ! ===========
@@ -320 +314 @@
             END DO
          END DO
-         !--- Lateral boundary conditions 
-         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
+         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions 
          !
          ! Tracer flux on the x-direction
@@ -353 +346 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )                     CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
          !
       END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
       !
    END SUBROUTINE tra_adv_qck_j
@@ -377 +368 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
+      !!----------------------------------------------------------------------
       !
       zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
@@ -421 +410 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwz )
-      !
    END SUBROUTINE tra_adv_cen2_k
 
@@ -443 +430 @@
       !----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('quickest')
+      IF( ln_timing )   CALL timing_start('quickest')
       !
       DO jk = 1, jpkm1
@@ -475 +462 @@
       END DO
       !
-      IF( nn_timing == 1 )  CALL timing_stop('quickest')
+      IF( ln_timing )   CALL timing_stop('quickest')
       !
    END SUBROUTINE quickest

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

@@ -20 +20 @@
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
-
    !
-   USE iom
-   USE lib_mpp        ! I/O library
+   USE iom            ! I/O library
+   USE lib_mpp        ! massively parallel library
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
    USE in_out_manager ! I/O manager
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -101 +99 @@
       REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
       REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztu, ztv, zltu, zltv, zti, ztw
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zltu, zltv, zti, ztw   ! 3D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_adv_ubs')
       !
       IF( kt == kit000 )  THEN
@@ -285 +281 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
+      IF( ln_timing )   CALL timing_stop('tra_adv_ubs')
       !
    END SUBROUTINE tra_adv_ubs
@@ -313 +307 @@
       INTEGER  ::   ikm1         ! local integer
       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zbetup, zbetdo )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zbetup, zbetdo     ! 3D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('nonosc_z')
       !
       zbig  = 1.e+40_wp
@@ -387 +379 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zbetup, zbetdo )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
+      IF( ln_timing )   CALL timing_stop('nonosc_z')
       !
    END SUBROUTINE nonosc_z

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

@@ -27 +27 @@
    USE lib_mpp        ! distributed memory computing library
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -77 +76 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt   ! 3D workspace
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_bbc')
+      IF( ln_timing )   CALL timing_start('tra_bbc')
       !
       IF( l_trdtra )   THEN         ! Save the input temperature trend
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt )
+         ALLOCATE( ztrdt(jpi,jpj,jpk) )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
       ENDIF
@@ -98 +97 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt )
+         DEALLOCATE( ztrdt )
       ENDIF
       !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_bbc')
+      IF( ln_timing )   CALL timing_stop('tra_bbc')
       !
    END SUBROUTINE tra_bbc
@@ -130 +129 @@
       TYPE(FLD_N)        ::   sn_qgh    ! informations about the geotherm. field to be read
       CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
-      !
+      !!
       NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
       !!----------------------------------------------------------------------
@@ -161 +160 @@
          !
          CASE ( 1 )                          !* constant flux
-            IF(lwp) WRITE(numout,*) '      *** constant heat flux  =   ', rn_geoflx_cst
+            IF(lwp) WRITE(numout,*) '      ===>>   constant heat flux  =   ', rn_geoflx_cst
             qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
             !
          CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
-            IF(lwp) WRITE(numout,*) '      *** variable geothermal heat flux'
+            IF(lwp) WRITE(numout,*) '      ===>>   variable geothermal heat flux'
             !
             ALLOCATE( sf_qgh(1), STAT=ierror )
@@ -173 +172 @@
             ENDIF
             ALLOCATE( sf_qgh(1)%fnow(jpi,jpj,1)   )
-            IF( sn_qgh%ln_tint )ALLOCATE( sf_qgh(1)%fdta(jpi,jpj,1,2) )
+            IF( sn_qgh%ln_tint )   ALLOCATE( sf_qgh(1)%fdta(jpi,jpj,1,2) )
             ! fill sf_chl with sn_chl and control print
             CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init',   &
@@ -187 +186 @@
          !
       ELSE
-         IF(lwp) WRITE(numout,*) '      *** no geothermal heat flux'
+         IF(lwp) WRITE(numout,*) '      ===>>   no geothermal heat flux'
       ENDIF
       !

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

@@ -13 +13 @@
    !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
    !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
+   !!            4.0  ! 2017-04  (G. Madec)  ln_trabbl namelist variable instead of a CPP key
    !!----------------------------------------------------------------------
-#if   defined key_trabbl
-   !!----------------------------------------------------------------------
-   !!   'key_trabbl'   or                             bottom boundary layer
+
    !!----------------------------------------------------------------------
    !!   tra_bbl_alloc : allocate trabbl arrays
@@ -36 +35 @@
    USE lbclnk         ! ocean lateral boundary conditions
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -49 +47 @@
    PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
-
    !                                !!* Namelist nambbl *
+   LOGICAL , PUBLIC ::   ln_trabbl   !: bottom boundary layer flag
    INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
    INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
@@ -82 +79 @@
       !!                ***  FUNCTION tra_bbl_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
-         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
-         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
-         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
+      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
+         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
+         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
+         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
          !
       IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
@@ -106 +103 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
-      !
-      IF( l_trdtra )   THEN                         !* Save the input trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'tra_bbl')
+      !
+      IF( l_trdtra )   THEN                         !* Save the T-S input trends
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
@@ -150 +147 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
+         DEALLOCATE( ztrdt, ztrds )
+      ENDIF
+      !
+      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
       !
    END SUBROUTINE tra_bbl
@@ -186 +183 @@
       INTEGER  ::   ik           ! local integers
       REAL(wp) ::   zbtr         ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zptb
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
-      !
-      CALL wrk_alloc( jpi, jpj, zptb )
+      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_bbl_dif')
       !
       DO jn = 1, kjpt                                     ! tracer loop
@@ -216 +211 @@
       END DO                                                ! end tracer
       !                                                     ! ===========
-      CALL wrk_dealloc( jpi, jpj, zptb )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
+      !
+      IF( ln_timing )   CALL timing_stop('tra_bbl_dif')
       !
    END SUBROUTINE tra_bbl_dif
@@ -249 +243 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
+      IF( ln_timing )   CALL timing_start( 'tra_bbl_adv')
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
@@ -301 +295 @@
             !
          END DO
-         !                                                       ! ===========
-      END DO                                                     ! end tracer
-      !                                                          ! ===========
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
+         !                                                  ! ===========
+      END DO                                                ! end tracer
+      !                                                     ! ===========
+      !
+      IF( ln_timing )   CALL timing_stop( 'tra_bbl_adv')
       !
    END SUBROUTINE tra_bbl_adv
@@ -349 +344 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
+      IF( ln_timing )   CALL timing_start( 'bbl')
       !
       IF( kt == kit000 )  THEN
@@ -480 +475 @@
       ENDIF
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
+      IF( ln_timing )   CALL timing_stop( 'bbl')
       !
    END SUBROUTINE bbl
@@ -494 +489 @@
       !!              called by nemo_init at the first timestep (kit000)
       !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj               ! dummy loop indices
-      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
-      INTEGER ::   ios                  !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
-      !
-      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
+      INTEGER ::   ji, jj                      ! dummy loop indices
+      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! workspace
+      !!
+      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'tra_bbl_init')
       !
       REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
@@ -519 +513 @@
          WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist nambbl : set bbl parameters'
-         WRITE(numout,*) '      diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
-         WRITE(numout,*) '      advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
-         WRITE(numout,*) '      diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
-         WRITE(numout,*) '      advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
-      ENDIF
-
+         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
+         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
+      ENDIF
+      IF( .NOT.ln_trabbl )   RETURN
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
+         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
+         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
+         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
+      ENDIF
+      !
       !                              ! allocate trabbl arrays
       IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
-
+      !
       IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
       IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
-
+      !
       !                             !* vertical index of  "deep" bottom u- and v-points
       DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
@@ -540 +539 @@
       END DO
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
-      CALL wrk_alloc( jpi, jpj, zmbk )
       zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
-      CALL wrk_dealloc( jpi, jpj, zmbk )
-
+      !
                                         !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
       mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
@@ -570 +567 @@
       ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
       ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
-
-      !
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
+      !
+      IF( ln_timing )   CALL timing_stop( 'tra_bbl_init')
       !
    END SUBROUTINE tra_bbl_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                      No bottom boundary layer scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
-CONTAINS
-   SUBROUTINE tra_bbl_init               ! Dummy routine
-   END SUBROUTINE tra_bbl_init
-   SUBROUTINE tra_bbl( kt )              ! Dummy routine
-      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
-   END SUBROUTINE tra_bbl
-#endif
 
    !!======================================================================

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

@@ -33 +33 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! XIOS
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
-   USE iom
 
    IMPLICIT NONE
@@ -94 +93 @@
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta, ztrdts
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('tra_dmp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,jpts,   zts_dta )
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)     ::  zts_dta
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ztrdts
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_dmp')
+      !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts ) 
+         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
          ztrdts(:,:,:,:) = tsa(:,:,:,:) 
       ENDIF
@@ -154 +153 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
-         CALL wrk_dealloc( jpi,jpj,jpk,jpts,   ztrdts ) 
+         DEALLOCATE( ztrdts ) 
       ENDIF
       !                           ! Control print
@@ -160 +159 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,jpts,   zts_dta )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('tra_dmp')
+      IF( ln_timing )   CALL timing_stop('tra_dmp')
       !
    END SUBROUTINE tra_dmp

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

@@ -30 +30 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
 
@@ -58 +57 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('tra_ldf')
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_ldf')
       !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt, ztrds ) 
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
@@ -85 +84 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt, ztrds ) 
+         DEALLOCATE( ztrdt, ztrds ) 
       ENDIF
       !                                        !* print mean trends (used for debugging)
@@ -91 +90 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('tra_ldf')
+      IF( ln_timing )   CALL timing_stop('tra_ldf')
       !
    END SUBROUTINE tra_ldf
@@ -107 +106 @@
       !!----------------------------------------------------------------------
       !
-      IF(lwp) THEN                     ! Namelist print
+      IF(lwp) THEN                     !==  Namelist print  ==!
          WRITE(numout,*)
          WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
@@ -114 +113 @@
          WRITE(numout,*) '      see ldf_tra_init report for lateral mixing parameters'
       ENDIF
-      !                                   ! use of lateral operator or not
+      !                                !==  use of lateral operator or not  ==!
       nldf   = np_ERROR
       ioptio = 0
-      IF( ln_traldf_lap )   ioptio = ioptio + 1
-      IF( ln_traldf_blp )   ioptio = ioptio + 1
-      IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
-      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral diffusion
-      !
-      IF( nldf /= np_no_ldf ) THEN        ! direction ==>> type of operator  
+      IF( ln_traldf_NONE ) THEN   ;   nldf = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_traldf_lap  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_traldf_blp  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF
+      IF( ioptio /=  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
+      !
+      IF( .NOT.ln_traldf_NONE ) THEN   !==  direction ==>> type of operator  ==!
          ioptio = 0
          IF( ln_traldf_lev )   ioptio = ioptio + 1
          IF( ln_traldf_hor )   ioptio = ioptio + 1
          IF( ln_traldf_iso )   ioptio = ioptio + 1
-         IF( ioptio >  1 )   CALL ctl_stop( 'tra_ldf_init: use only ONE direction (level/hor/iso)' )
+         IF( ioptio /=  1  )   CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' )
          !
          !                                ! defined the type of lateral diffusion from ln_traldf_... logicals

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

@@ -30 +30 @@
    USE phycst         ! physical constants
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -111 +110 @@
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw 
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
-      !
-      CALL wrk_alloc( jpi,jpj,       zdkt, zdk1t, z2d ) 
-      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt , zftu, zftv, ztfw  ) 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_ldf_iso')
       !
       IF( kt == kit000 )  THEN
@@ -386 +382 @@
          !                                                        ! ===============
       END DO                                                      ! end tracer loop
-      !                                                           ! ===============
-      !
-      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw  ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
+      !
+      IF( ln_timing )   CALL timing_stop('tra_ldf_iso')
       !
    END SUBROUTINE tra_ldf_iso

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90

@@ -22 +22 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! distribued memory computing library
    USE timing         ! Timing
-   USE wrk_nemo       ! Memory allocation
-   USE iom
 
    IMPLICIT NONE
@@ -87 +86 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       REAL(wp) ::   zsign            ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztu, ztv, zaheeu, zaheev
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('tra_ldf_lap')
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zaheeu, zaheev
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_ldf_lap')
       !
       IF( kt == nit000 .AND. lwp )  THEN
@@ -97 +96 @@
          WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
       !
       l_hst = .FALSE.
@@ -169 +166 @@
       !                             ! ==================
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_lap')
+      IF( ln_timing )   CALL timing_stop('tra_ldf_lap')
       !
    END SUBROUTINE tra_ldf_lap
@@ -203 +198 @@
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zlap         ! laplacian at t-point
-      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
-      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) :: zlap         ! laplacian at t-point
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_blp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,kjpt,   zlap ) 
-      CALL wrk_alloc( jpi,jpj,    kjpt,   zglu, zglv, zgui, zgvi ) 
+      IF( ln_timing )   CALL timing_start('tra_ldf_blp')
       !
       IF( kt == kit000 .AND. lwp )  THEN
@@ -253 +245 @@
       END SELECT
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,kjpt,   zlap ) 
-      CALL wrk_dealloc( jpi,jpj    ,kjpt,   zglu, zglv, zgui, zgvi ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_blp')
+      IF( ln_timing )   CALL timing_stop('tra_ldf_blp')
       !
    END SUBROUTINE tra_ldf_blp

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90

@@ -27 +27 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -94 +93 @@
       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
       REAL(wp) ::   zah, zah_slp, zaei_slp
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   z2d                                              ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_triad')
-      !
-      CALL wrk_alloc( jpi,jpj,       z2d ) 
-      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  ) 
+      IF( ln_timing )   CALL timing_start('tra_ldf_triad')
       !
       IF( .NOT.ALLOCATED(zdkt3d) )  THEN
@@ -434 +430 @@
       END DO                                                      ! end tracer loop
       !                                                           ! ===============
-      !
-      CALL wrk_dealloc( jpi,jpj,       z2d ) 
-      CALL wrk_dealloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_triad')
+      IF( ln_timing )   CALL timing_stop('tra_ldf_triad')
       !
    END SUBROUTINE tra_ldf_triad

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

@@ -26 +26 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -67 +66 @@
       REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
       REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
-      REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0)
-      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point...
-      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point...
-      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2 
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace
+      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0)
+      REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point...
+      REAL(wp), DIMENSION(        jpk,jpts) ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ) ::   zn2          ! N^2 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zab          ! alpha and beta
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
       !
       LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
@@ -80 +78 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_npc')
+      IF( ln_timing )   CALL timing_start('tra_npc')
       !
       IF( MOD( kt, nn_npc ) == 0 ) THEN
          !
-         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2
-         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta
-         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj
-         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj
-
          IF( l_trdtra )   THEN                    !* Save initial after fields
-            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
             ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
             ztrds(:,:,:) = tsa(:,:,:,jp_sal)
          ENDIF
-
+         !
          IF( l_LB_debug ) THEN
             ! Location of 1 known convection site to follow what's happening in the water column
@@ -101 +94 @@
             klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...          
          ENDIF
-         
+         !
          CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points)
          CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points)
-        
+         !
          inpcc = 0
-
+         !
          DO jj = 2, jpjm1                 ! interior column only
             DO ji = fs_2, fs_jpim1
@@ -313 +306 @@
             CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
-            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+            DEALLOCATE( ztrdt, ztrds )
          ENDIF
          !
@@ -323 +316 @@
          ENDIF
          !
-         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
-         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
-         CALL wrk_dealloc(jpk, zvn2 )
-         CALL wrk_dealloc(jpk, 2, zvts, zvab )
-         !
       ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
+      IF( ln_timing )   CALL timing_stop('tra_npc')
       !
    END SUBROUTINE tra_npc

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

@@ -35 +35 @@
    USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE phycst          ! physical constant
-   USE ldftra          ! lateral physics on tracers
-   USE ldfslp
-   USE bdy_oce   , ONLY: ln_bdy
+   USE ldftra          ! lateral physics : tracers
+   USE ldfslp          ! lateral physics : slopes
+   USE bdy_oce  , ONLY : ln_bdy
    USE bdytra          ! open boundary condition (bdy_tra routine)
    !
@@ -43 +43 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
 #if defined key_agrif
@@ -91 +90 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       REAL(wp) ::   zfact            ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'tra_nxt')
       !
       IF( kt == nit000 ) THEN
@@ -114 +113 @@
  
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =     rdt      ! at nit000             (Euler)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
@@ -120 +119 @@
       ! trends computation initialisation
       IF( l_trdtra )   THEN                    
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,jpk) = 0._wp
          ztrds(:,:,jpk) = 0._wp
@@ -136 +135 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
-         IF( ln_linssh ) THEN 
+         IF( ln_linssh ) THEN       ! linear sea surface height only
             ! Store now fields before applying the Asselin filter 
             ! in order to calculate Asselin filter trend later.
@@ -150 +149 @@
             END DO
          END DO
-         IF (l_trdtra .AND. .NOT. ln_linssh) THEN  ! Zero Asselin filter contribution must be explicitly written out since for vvl
-                                                   ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
+         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
+            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
             ztrdt(:,:,:) = 0._wp
             ztrds(:,:,:) = 0._wp
@@ -181 +180 @@
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
       END IF
-      IF( l_trdtra ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
       !
       !                        ! control print
@@ -187 +186 @@
          &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
+      IF( ln_timing )   CALL timing_stop('tra_nxt')
       !
    END SUBROUTINE tra_nxt
@@ -271 +270 @@
       REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
       REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
       !!----------------------------------------------------------------------
       !
@@ -290 +289 @@
       ENDIF
       !
-      IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) )   THEN
-         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
          ztrd_atf(:,:,:,:) = 0.0_wp
       ENDIF
       zfact = 1._wp / r2dt
+      zfact1 = atfp * p2dt
+      zfact2 = zfact1 * r1_rau0
       DO jn = 1, kjpt      
          DO jk = 1, jpkm1
-            zfact1 = atfp * p2dt
-            zfact2 = zfact1 * r1_rau0
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -357 +356 @@
       END DO
       !
-      IF( l_trdtra .and. cdtype == 'TRA' ) THEN 
-         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
-         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
-         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
-      ENDIF
-      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
-         DO jn = 1, kjpt
-            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
-         END DO
-         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
+         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
+            CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
+            CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
+         ENDIF
+         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
+            DO jn = 1, kjpt
+               CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
+            END DO
+         ENDIF
+         DEALLOCATE( ztrd_atf )
       ENDIF
       !

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -29 +29 @@
    USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
-   USE iom            ! I/O manager
+   USE iom            ! I/O library
    USE fldread        ! read input fields
    USE restart        ! ocean restart
    USE lib_mpp        ! MPP library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -48 +47 @@
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
-   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
@@ -113 +111 @@
       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
       REAL(wp) ::   zlogc, zlogc2, zlogc3 
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot, zchl3d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_qsr')
       !
       IF( kt == nit000 ) THEN
@@ -127 +125 @@
       !
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
+         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
       ENDIF
@@ -161 +159 @@
       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
-         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
+         ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , &
+            &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , &
+            &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   ) 
          !
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
@@ -240 +239 @@
          END DO
          !
-         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
-         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
+         DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) 
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
@@ -269 +267 @@
       END DO
       !
-      IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
-               ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
-               ENDIF
-            END DO
-         END DO
-         ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere
-         CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp )
-      ENDIF
+      ! sea-ice: store the 1st ocean level attenuation coefficient
+      DO jj = 2, jpjm1 
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
+            ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
+            ENDIF
+         END DO
+      END DO
+      CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp )
       !
       IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
-         CALL wrk_alloc( jpi,jpj,jpk,   zetot )
-         !
+         ALLOCATE( zetot(jpi,jpj,jpk) )
          zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
          DO jk = nksr, 1, -1
-            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp
+            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp
          END DO         
          CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
-         !
-         CALL wrk_dealloc( jpi,jpj,jpk,   zetot ) 
+         DEALLOCATE( zetot ) 
       ENDIF
       !
@@ -301 +295 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt ) 
+         DEALLOCATE( ztrdt ) 
       ENDIF
       !                       ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
+      IF( ln_timing )   CALL timing_stop('tra_qsr')
       !
    END SUBROUTINE tra_qsr
@@ -336 +330 @@
       TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read
       !!
-      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  &
+      NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio,  &
          &                  nn_chldta, rn_abs, rn_si0, rn_si1
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init')
+      IF( ln_timing )   CALL timing_start('tra_qsr_init')
       !
       REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
@@ -359 +353 @@
          WRITE(numout,*) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
          WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
-         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice
          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
          WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
@@ -435 +428 @@
       ENDIF
       !
-      IF( nn_timing == 1 )   CALL timing_stop('tra_qsr_init')
+      IF( ln_timing )   CALL timing_stop('tra_qsr_init')
       !
    END SUBROUTINE tra_qsr_init

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

@@ -32 +32 @@
    USE iom            ! xIOS server
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -75 +74 @@
       INTEGER  ::   ikt, ikb              ! local integers
       REAL(wp) ::   zfact, z1_e3t, zdep   ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
+      IF( ln_timing )   CALL timing_start('tra_sbc')
       !
       IF( kt == nit000 ) THEN
@@ -87 +86 @@
       !
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
@@ -232 +231 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
+         DEALLOCATE( ztrdt , ztrds ) 
       ENDIF
       !
@@ -238 +237 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
+      IF( ln_timing )   CALL timing_stop('tra_sbc')
       !
    END SUBROUTINE tra_sbc

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

@@ -4 +4 @@
    !! Ocean active tracers:  vertical component of the tracer mixing trend
    !!==============================================================================
-   !! History :  1.0  ! 2005-11  (G. Madec)  Original code
-   !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
+   !! History :  1.0  !  2005-11  (G. Madec)  Original code
+   !!            3.0  !  2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            4.0  !  2017-06  (G. Madec)  remove explict time-stepping option
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   tra_zdf       : Update the tracer trend with the vertical diffusion
-   !!   tra_zdf_init  : initialisation of the computation
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -20 +20 @@
    USE ldftra         ! lateral diffusion: eddy diffusivity
    USE ldfslp         ! lateral diffusion: iso-neutral slope 
-   USE trazdf_exp     ! vertical diffusion: explicit (tra_zdf_exp routine)
-   USE trazdf_imp     ! vertical diffusion: implicit (tra_zdf_imp routine)
    USE trd_oce        ! trends: ocean variables
    USE trdtra         ! trends: tracer trend manager
@@ -29 +27 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
 
@@ -35 +32 @@
    PRIVATE
 
-   PUBLIC   tra_zdf        ! routine called by step.F90
-   PUBLIC   tra_zdf_init   ! routine called by nemogcm.F90
-
-   INTEGER ::   nzdf = 0   ! type vertical diffusion algorithm used (defined from ln_zdf...  namlist logicals)
+   PUBLIC   tra_zdf       ! called by step.F90
+   PUBLIC   tra_zdf_imp   ! called by trczdf.F90
 
    !! * Substitutions
-#  include "zdfddm_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -56 +50 @@
       !! ** Purpose :   compute the vertical ocean tracer physics.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !
-      INTEGER  ::   jk                   ! Dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_zdf')
-      !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
-         r2dt =  rdt                          ! = rdt (restarting with Euler time stepping)
-      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
-         r2dt = 2. * rdt                      ! = 2 rdt (leapfrog)
-      ENDIF
-      !
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   jk   ! Dummy loop indices
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_zdf')
+      !
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping)
+      ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog)
+      ENDIF
+      !
+      IF( l_trdtra )   THEN                  !* Save ta and sa trends
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
       !
-      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
-      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
-      END SELECT
+      !                                      !* compute lateral mixing trend and add it to the general trend
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 
+
 !!gm WHY here !   and I don't like that !
       ! DRAKKAR SSS control {
@@ -90 +81 @@
          DO jk = 1, jpkm1
             ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &
-                 & / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)
+               &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)
             ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &
-                 & / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)
+              &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
@@ -100 +91 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+         DEALLOCATE( ztrdt , ztrds )
       ENDIF
       !                                          ! print mean trends (used for debugging)
@@ -106 +97 @@
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf')
+      IF( ln_timing )   CALL timing_stop('tra_zdf')
       !
    END SUBROUTINE tra_zdf
 
-
-   SUBROUTINE tra_zdf_init
+ 
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 
       !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE tra_zdf_init  ***
-      !!
-      !! ** Purpose :   Choose the vertical mixing scheme
-      !!
-      !! ** Method  :   Set nzdf from ln_zdfexp
-      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T)
-      !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F)
-      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes,
-      !!          an implicit scheme is required.
-      !!----------------------------------------------------------------------
-      USE zdftke
-      USE zdfgls
-      !!----------------------------------------------------------------------
-      !
-      ! Choice from ln_zdfexp already read in namelist in zdfini module
-      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
-      ELSE                   ;   nzdf = 1           ! use implicit scheme
-      ENDIF
-      !
-      ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics
-      IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
-      IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
-      IF( ln_zdfexp .AND. nzdf == 1 )   CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator',   &
-            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
-            !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
-      ENDIF
-      !
-   END SUBROUTINE tra_zdf_init
+      !!                  ***  ROUTINE tra_zdf_imp  ***
+      !!
+      !! ** Purpose :   Compute the after tracer through a implicit computation
+      !!     of the vertical tracer diffusion (including the vertical component 
+      !!     of lateral mixing (only for 2nd order operator, for fourth order 
+      !!     it is already computed and add to the general trend in traldf) 
+      !!
+      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
+      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
+      !!      It is computed using a backward time scheme (t=after field)
+      !!      which provide directly the after tracer field.
+      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers
+      !!      Surface and bottom boundary conditions: no diffusive flux on
+      !!      both tracers (bottom, applied through the masked field avt).
+      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
+      !!
+      !! ** Action  : - pta  becomes the after tracer
+      !!---------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
+      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field
+      !
+      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::  zrhs             ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('tra_zdf_imp')
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
+         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
+      !                                               ! ============= !
+      DO jn = 1, kjpt                                 !  tracer loop  !
+         !                                            ! ============= !
+         !  Matrix construction
+         ! --------------------
+         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
+         !
+         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   &
+            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
+            !
+            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
+            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
+            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk)
+            ENDIF
+            zwt(:,:,1) = 0._wp
+            !
+            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
+               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
+                        END DO
+                     END DO
+                  END DO
+               ELSE                          ! standard or triad iso-neutral operator
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
+                        END DO
+                     END DO
+                  END DO
+               ENDIF
+            ENDIF
+            !
+            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that
+                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
+                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                 END DO
+               END DO
+            END DO
+            !
+            !! Matrix inversion from the first level
+            !!----------------------------------------------------------------------
+            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+            !
+            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+            !        (        ...               )( ...  ) ( ...  )
+            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+            !
+            !   m is decomposed in the product of an upper and lower triangular matrix.
+            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
+            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
+            !   and "superior" (above diagonal) components of the tridiagonal system.
+            !   The solution will be in the 4d array pta.
+            !   The 3d array zwt is used as a work space array.
+            !   En route to the solution pta is used a to evaluate the rhs and then 
+            !   used as a work space array: its value is modified.
+            !
+            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
+               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
+                  zwt(ji,jj,1) = zwd(ji,jj,1)
+               END DO
+            END DO
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF 
+         !         
+         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+            DO ji = fs_2, fs_jpim1
+               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side
+                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
+               END DO
+            END DO
+         END DO
+         !
+         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
+            DO ji = fs_2, fs_jpim1
+               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+            END DO
+         END DO
+         DO jk = jpk-2, 1, -1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
+                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                                            ! ================= !
+      END DO                                          !  end tracer loop  !
+      !                                               ! ================= !
+      !
+      IF( ln_timing )   CALL timing_stop('tra_zdf_imp')
+      !
+   END SUBROUTINE tra_zdf_imp
 
    !!==============================================================================

branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

@@ -22 +22 @@
    USE lbclnk          ! lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
 
@@ -99 +98 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde')
-      !
-      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp
-      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp
+      IF( ln_timing )   CALL timing_start( 'zps_hde')
+      !
+      pgtu(:,:,:) = 0._wp   ;   zti (:,:,:) = 0._wp   ;   zhi (:,:) = 0._wp
+      pgtv(:,:,:) = 0._wp   ;   ztj (:,:,:) = 0._wp   ;   zhj (:,:) = 0._wp
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -188 +187 @@
       END IF
       !
-      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde')
+      IF( ln_timing )   CALL timing_stop( 'zps_hde')
       !
    END SUBROUTINE zps_hde
-   !
+
+
    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
       &                          prd, pgru, pgrv, pgrui, pgrvi )
@@ -256 +256 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
+      IF( ln_timing )   CALL timing_start( 'zps_hde_isf')
       !
       pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
@@ -453 +453 @@
       END IF  
       !
-      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf')
+      IF( ln_timing )   CALL timing_stop( 'zps_hde_isf')
       !
    END SUBROUTINE zps_hde_isf
+
    !!======================================================================
 END MODULE zpshde