Changeset 8568 for branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN

Timestamp:

2017-09-27T16:29:24+02:00 (7 years ago)

Author:

gm

Message:

#1911 (ENHANCE-09): PART I.2 - _NONE option + remove zts + see associated wiki page

Location:

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN

Files:

• divhor.F90 (modified) (5 diffs)
• dynadv.F90 (modified) (6 diffs)
• dynadv_cen2.F90 (modified) (4 diffs)
• dynadv_ubs.F90 (modified) (4 diffs)
• dynbfr.F90 (modified) (2 diffs)
• dynhpg.F90 (modified) (26 diffs)
• dynkeg.F90 (modified) (9 diffs)
• dynldf.F90 (modified) (7 diffs)
• dynldf_iso.F90 (modified) (6 diffs)
• dynldf_lap_blp.F90 (modified) (7 diffs)
• dynnxt.F90 (modified) (6 diffs)
• dynspg.F90 (modified) (9 diffs)
• dynspg_exp.F90 (modified) (2 diffs)
• dynspg_ts.F90 (modified) (2 diffs)
• dynvor.F90 (modified) (15 diffs)
• dynzad.F90 (modified) (8 diffs)
• dynzdf.F90 (modified) (2 diffs)
• sshwzv.F90 (modified) (10 diffs)
• wet_dry.F90 (modified) (7 diffs)

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90

@@ -29 +29 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
 
@@ -40 +39 @@
 #  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)
@@ -64 +63 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('div_hor')
+      IF( ln_timing )   CALL timing_start('div_hor')
       !
       IF( kt == nit000 ) THEN
@@ -75 +74 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               hdivn(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * un(ji  ,jj,jk)        &
-                  &               - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * un(ji-1,jj,jk)        &
-                  &               + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vn(ji,jj  ,jk)        &
-                  &               - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk)   )    &
-                  &            / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
+               hdivn(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * un(ji  ,jj,jk)      &
+                  &               - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * un(ji-1,jj,jk)      &
+                  &               + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vn(ji,jj  ,jk)      &
+                  &               - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk)  )   &
+                  &            * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
             END DO  
          END DO  
@@ -90 +89 @@
       END DO
       !
-      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field)
+      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )              !==  runoffs    ==!   (update hdivn field)
       !
-      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL sbc_isf_div( hdivn )              !==  ice shelf  ==!   (update hdivn field)
       !
-      IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn )  !==  ice sheet  ==!   (update hdivn field)
+      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn ) !==  ice sheet  ==!   (update hdivn field)
       !
-      CALL lbc_lnk( hdivn, 'T', 1. )                !==  lateral boundary cond.  ==!   (no sign change)
+      CALL lbc_lnk( hdivn, 'T', 1. )   !   (no sign change)
       !
-      IF( nn_timing == 1 )  CALL timing_stop('div_hor')
+      IF( ln_timing )   CALL timing_stop('div_hor')
       !
    END SUBROUTINE div_hor

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

@@ -7 +7 @@
    !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
    !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option 
+   !!            4.0  !  2017-07  (G. Madec)  add a linear dynamics option
    !!----------------------------------------------------------------------
 
@@ -30 +31 @@
  
    !                                    !* namdyn_adv namelist *
-   LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form flag
-   INTEGER, PUBLIC ::   nn_dynkeg       !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth
+   LOGICAL, PUBLIC ::   ln_dynadv_NONE  !: linear dynamics (no momentum advection)
+   LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form
+   INTEGER, PUBLIC ::      nn_dynkeg       !: scheme of grad(KE): =0 C2 ; =1 Hollingsworth
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
-   LOGICAL, PUBLIC ::   ln_dynzad_zts   !: vertical advection with sub-timestepping (requires vector form)
    
-   INTEGER ::   nadv   ! choice of the formulation and scheme for the advection
+   INTEGER, PUBLIC ::   n_dynadv   !: choice of the formulation and scheme for momentum advection
+   !                               !  associated indices:
+   INTEGER, PUBLIC, PARAMETER ::   np_LIN_dyn = 0   ! no advection: linear dynamics
+   INTEGER, PUBLIC, PARAMETER ::   np_VEC_c2  = 1   ! vector form : 2nd order centered scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_FLX_c2  = 2   ! flux   form : 2nd order centered scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_FLX_ubs = 3   ! flux   form : 3rd order Upstream Biased Scheme
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -53 +59 @@
       !! ** Purpose :   compute the ocean momentum advection trend.
       !!
-      !! ** Method  : - Update (ua,va) with the advection term following nadv
+      !! ** Method  : - Update (ua,va) with the advection term following n_dynadv
+      !!
       !!      NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T) 
       !!      a metric term is add to the coriolis term while in vector form 
@@ -62 +69 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_adv')
+      IF( ln_timing )   CALL timing_start( 'dyn_adv' )
       !
-      SELECT CASE ( nadv )                  ! compute advection trend and add it to general trend
-      CASE ( 0 )     
-                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad     ( kt )               ! vector form : vertical advection
-      CASE ( 1 )     
-                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad_zts ( kt )               ! vector form : vertical advection with sub-timestepping
-      CASE ( 2 ) 
-                      CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme
-      CASE ( 3 )   
-                      CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme
+      SELECT CASE( n_dynadv )    !==  compute advection trend and add it to general trend  ==!
+      CASE( np_VEC_c2  )     
+         CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
+         CALL dyn_zad     ( kt )               ! vector form : vertical advection
+      CASE( np_FLX_c2  ) 
+         CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme
+      CASE( np_FLX_ubs )   
+         CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme (UP3)
       END SELECT
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv')
+      IF( ln_timing )   CALL timing_stop( 'dyn_adv' )
       !
    END SUBROUTINE dyn_adv
@@ -87 +91 @@
       !!                
       !! ** Purpose :   Control the consistency between namelist options for 
-      !!              momentum advection formulation & scheme and set nadv
+      !!              momentum advection formulation & scheme and set n_dynadv
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio, ios   ! Local integer
       !
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
+      NAMELIST/namdyn_adv/ ln_dynadv_NONE, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs
       !!----------------------------------------------------------------------
       !
@@ -108 +112 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
-         WRITE(numout,*) '      Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
-         WRITE(numout,*) '      = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
-         WRITE(numout,*) '      2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
-         WRITE(numout,*) '      3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
-         WRITE(numout,*) '      Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
+         WRITE(numout,*) '      linear dynamics : no momentum advection          ln_dynadv_NONE = ', ln_dynadv_NONE
+         WRITE(numout,*) '      Vector form: 2nd order centered scheme           ln_dynadv_vec  = ', ln_dynadv_vec
+         WRITE(numout,*) '         with Hollingsworth scheme (=1) or not (=0)       nn_dynkeg   = ', nn_dynkeg
+         WRITE(numout,*) '      flux form: 2nd order centred scheme              ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '                 3rd order UBS scheme                  ln_dynadv_ubs  = ', ln_dynadv_ubs
       ENDIF
 
-      ioptio = 0                      ! Parameter control
-      IF( ln_dynadv_vec  )   ioptio = ioptio + 1
-      IF( ln_dynadv_cen2 )   ioptio = ioptio + 1
-      IF( ln_dynadv_ubs  )   ioptio = ioptio + 1
+      ioptio = 0                      ! parameter control and set n_dynadv
+      IF( ln_dynadv_NONE ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_LIN_dyn   ;   ENDIF
+      IF( ln_dynadv_vec  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_VEC_c2    ;   ENDIF
+      IF( ln_dynadv_cen2 ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_FLX_c2    ;   ENDIF
+      IF( ln_dynadv_ubs  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_FLX_ubs   ;   ENDIF
 
-      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
-      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
-         CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
-      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   &  
-         CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
+      IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' )
+      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
 
-      !                               ! Set nadv
-      IF( ln_dynadv_vec  )   nadv =  0 
-      IF( ln_dynzad_zts  )   nadv =  1
-      IF( ln_dynadv_cen2 )   nadv =  2
-      IF( ln_dynadv_ubs  )   nadv =  3
 
       IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '      ===>>   vector form : keg + zad + vor is used' 
-         IF( nadv ==  1 )   WRITE(numout,*) '      ===>>   vector form : keg + zad_zts + vor is used'
-         IF( nadv ==  0 .OR. nadv ==  1 ) THEN
+         SELECT CASE( n_dynadv )
+         CASE( np_LIN_dyn )   ;   WRITE(numout,*) '      ===>>   linear dynamics : no momentum advection used'
+         CASE( np_VEC_c2  )   ;   WRITE(numout,*) '      ===>>   vector form : keg + zad + vor is used' 
             IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) '              with Centered standard keg scheme'
             IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) '              with Hollingsworth keg scheme'
-         ENDIF
-         IF( nadv ==  2 )   WRITE(numout,*) '      ===>>   flux form   : 2nd order scheme is used'
-         IF( nadv ==  3 )   WRITE(numout,*) '      ===>>   flux form   : UBS       scheme is used'
+         CASE( np_FLX_c2  )   ;   WRITE(numout,*) '      ===>>   flux form   : 2nd order scheme is used'
+         CASE( np_FLX_ubs )   ;   WRITE(numout,*) '      ===>>   flux form   : UBS       scheme is used'
+         END SELECT
       ENDIF
       !

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90

@@ -20 +20 @@
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -31 +30 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -51 +50 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zfu, zfv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfu_t, zfu_f, zfu_uw, zfu
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zfv_t, zfv_f, zfv_vw, zfv, zfw
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_cen2')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
+      IF( ln_timing )   CALL timing_start('dyn_adv_cen2')
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -148 +145 @@
          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_cen2')
+      IF( ln_timing )   CALL timing_stop('dyn_adv_cen2')
       !
    END SUBROUTINE dyn_adv_cen2

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

@@ -23 +23 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -37 +36 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -74 +73 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu, zfv
-      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zlu_uu, zlv_vv, zlu_uv, zlv_vu
+      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
+      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_ubs')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,        zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
-      CALL wrk_alloc( jpi,jpj,jpk,jpts,   zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
+      IF( ln_timing )   CALL timing_start('dyn_adv_ubs')
       !
       IF( kt == nit000 ) THEN
@@ -241 +238 @@
          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,        zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
-      CALL wrk_dealloc( jpi,jpj,jpk,jpts,   zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
+      IF( ln_timing )   CALL timing_stop('dyn_adv_ubs')
       !
    END SUBROUTINE dyn_adv_ubs

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

@@ -57 +57 @@
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_bfr')
+      IF( ln_timing )   CALL timing_start('dyn_bfr')
       !
 !!gm bug : time step is only rdt (not 2 rdt if euler start !)
@@ -109 +109 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_bfr')
+      IF( ln_timing )   CALL timing_stop('dyn_bfr')
       !
    END SUBROUTINE dyn_bfr

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

@@ -44 +44 @@
    USE lib_mpp         ! MPP library
    USE eosbn2          ! compute density
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
    USE iom
@@ -84 +83 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_hpg')
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_hpg')
       !
       IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
-         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
@@ -108 +107 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         DEALLOCATE( ztrdu , ztrdv )
       ENDIF
       !
@@ -114 +113 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg')
+      IF( ln_timing )   CALL timing_stop('dyn_hpg')
       !
    END SUBROUTINE dyn_hpg
@@ -134 +133 @@
       INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
       REAL(wp) ::   znad
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zts_top, zrhd   ! hypothesys on isf density
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  zrhdtop_isf    ! density at bottom of ISF
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  ziceload       ! density at bottom of ISF
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
@@ -165 +164 @@
       !
       IF( ln_hpg_djc )   &
-         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method &
-                           & currently disabled (bugs under investigation). Please select &
-                           & either  ln_hpg_sco or  ln_hpg_prj instead')
-      !
-      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )        &
-         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', &
-         &                 '   the standard jacobian formulation hpg_sco    or '    , &
-         &                 '   the pressure jacobian formulation hpg_prj'            )
-
-      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   &
-         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
-      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   &
-         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
+         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method',   &
+         &                 '   currently disabled (bugs under investigation).'        ,   &
+         &                 '   Please select either  ln_hpg_sco or  ln_hpg_prj instead' )
+         !
+      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )          &
+         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   &
+         &                 '   the standard jacobian formulation hpg_sco    or '    ,   &
+         &                 '   the pressure jacobian formulation hpg_prj'             )
+         !
+      IF( ln_hpg_isf ) THEN
+         IF( .NOT. ln_isfcav )   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
+       ELSE
+         IF(       ln_isfcav )   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
+      ENDIF
       !
       !                               ! Set nhpg from ln_hpg_... flags
@@ -197 +197 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
       ! 
-      ! initialisation of ice shelf load
-      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0
-      IF (       ln_isfcav ) THEN
-         CALL wrk_alloc( jpi,jpj, 2,  ztstop) 
-         CALL wrk_alloc( jpi,jpj,jpk, zrhd  )
-         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload) 
+      !                          
+      IF ( .NOT. ln_isfcav ) THEN     !--- no ice shelf load
+         riceload(:,:) = 0._wp
+         !
+      ELSE                            !--- set an ice shelf load
          !
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'  
-
-         ! To use density and not density anomaly
-         znad=1._wp
-         
-         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
-
-         ! compute density of the water displaced by the ice shelf 
-         DO jk = 1, jpk
-            CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk))
-         END DO
-      
-         ! compute rhd at the ice/oce interface (ice shelf side)
-         CALL eos(ztstop,risfdep,zrhdtop_isf)
-
-         ! Surface value + ice shelf gradient
-         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
-         ! divided by 2 later
-         ziceload = 0._wp
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikt=mikt(ji,jj)
+         IF(lwp) WRITE(numout,*) '   ice shelf case: set the ice-shelf load'
+         ALLOCATE( zts_top(jpi,jpj,jpts) , zrhd(jpi,jpj,jpk) , zrhdtop_isf(jpi,jpj) , ziceload(jpi,jpj) ) 
+         !
+         znad = 1._wp                     !- To use density and not density anomaly
+         !
+         !                                !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+         zts_top(:,:,jp_tem) = -1.9_wp   ;   zts_top(:,:,jp_sal) = 34.4_wp
+         !
+         DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
+            CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )
+         END DO
+         !
+         !                                !- compute rhd at the ice/oce interface (ice shelf side)
+         CALL eos( zts_top , risfdep, zrhdtop_isf )
+         !
+         !                                !- Surface value + ice shelf gradient
+         ziceload = 0._wp                       ! compute pressure due to ice shelf load 
+         DO jj = 1, jpj                         ! (used to compute hpgi/j for all the level from 1 to miku/v)
+            DO ji = 1, jpi                      ! divided by 2 later
+               ikt = mikt(ji,jj)
                ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
-               DO jk=2,ikt-1
+               DO jk = 2, ikt-1
                   ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
                      &                              * (1._wp - tmask(ji,jj,jk))
                END DO
                IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
-            END DO
-         END DO
-         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
-
-         CALL wrk_dealloc( jpi,jpj, 2,  ztstop) 
-         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  )
-         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload) 
-      END IF
+                  &                                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+            END DO
+         END DO
+         riceload(:,:) = ziceload(:,:)  ! need to be saved for diaar5
+         !
+         DEALLOCATE( zts_top , zrhd , zrhdtop_isf , ziceload ) 
+      ENDIF
       !
    END SUBROUTINE dyn_hpg_init
@@ -268 +261 @@
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zhpj )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
+      !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN
@@ -315 +306 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zhpi, zhpj )
-      !
    END SUBROUTINE hpg_zco
 
@@ -333 +322 @@
       INTEGER  ::   iku, ikv                         ! temporary integers
       REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zhpj )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
+      !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN
@@ -405 +392 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zhpi, zhpj )
-      !
    END SUBROUTINE hpg_zps
 
@@ -433 +418 @@
       REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
       LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
-      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
+      !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN
@@ -452 +434 @@
       !
       IF( ln_wd ) THEN
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1 
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 
+               ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
-             IF(ll_tmp1) THEN
-               zcpx(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
-             ELSE
-               zcpx(ji,jj) = 0._wp
-             END IF
-      
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+               IF(ll_tmp1) THEN
+                  zcpx(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+                  zcpx(ji,jj) = ABS(   ( sshn(ji+1,jj)+ht_wd(ji+1,jj) - sshn(ji,jj)-ht_wd(ji,jj) )   &
+                     &                / ( sshn(ji+1,jj)                - sshn(ji,jj)              )   )
+               ELSE
+                  zcpx(ji,jj) = 0._wp
+               ENDIF
+               !
+               ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-
-             IF(ll_tmp1) THEN
-               zcpy(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
-             ELSE
-               zcpy(ji,jj) = 0._wp
-             END IF
-           END DO
-        END DO
-        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      END IF
+                  !
+               IF(ll_tmp1) THEN
+                  zcpy(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+                  zcpy(ji,jj) = ABS(   ( sshn(ji,jj+1)+ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj) )   &
+                     &               / ( sshn(ji,jj+1)                - sshn(ji,jj)                )   )
+               ELSE
+                  zcpy(ji,jj) = 0._wp
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
 
       ! Surface value
@@ -507 +490 @@
             zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
                &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
-
-
+            !
             IF( ln_wd ) THEN
-
-              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
-              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
-              zuap = zuap * zcpx(ji,jj)
-              zvap = zvap * zcpy(ji,jj)
+               zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
+               zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+               zuap = zuap * zcpx(ji,jj)
+               zvap = zvap * zcpy(ji,jj)
             ENDIF
-
+            !
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
@@ -539 +520 @@
                zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
-
+               !
                IF( ln_wd ) THEN
-                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
-                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
-                 zuap = zuap * zcpx(ji,jj)
-                 zvap = zvap * zcpy(ji,jj)
-               ENDIF
-
+                  zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
+                  zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+                  zuap = zuap * zcpx(ji,jj)
+                  zvap = zvap * zcpy(ji,jj)
+               ENDIF
+               !
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
@@ -554 +535 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
-      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd )   DEALLOCATE( zcpx , zcpy )
       !
    END SUBROUTINE hpg_sco
@@ -583 +563 @@
       INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
       REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
-      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,  2, ztstop) 
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj)
-      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce )
-      !
-      ! Local constant initialization
-      zcoef0 = - grav * 0.5_wp
-  
-      ! To use density and not density anomaly
-      znad=1._wp
-
-      ! iniitialised to 0. zhpi zhpi 
-      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
+      REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
+      REAL(wp), DIMENSION(jpi,jpj)      ::  zrhdtop_oce
+      !!----------------------------------------------------------------------
+      !
+      zcoef0 = - grav * 0.5_wp   ! Local constant initialization
+      !
+      znad=1._wp                 ! To use density and not density anomaly
+      !
+      !                          ! iniitialised to 0. zhpi zhpi 
+      zhpi(:,:,:) = 0._wp   ;   zhpj(:,:,:) = 0._wp
 
       ! compute rhd at the ice/oce interface (ocean side)
       ! usefull to reduce residual current in the test case ISOMIP with no melting
-      DO ji=1,jpi
-        DO jj=1,jpj
-          ikt=mikt(ji,jj)
-          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
-          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
+      DO ji = 1, jpi
+        DO jj = 1, jpj
+          ikt = mikt(ji,jj)
+          zts_top(ji,jj,1) = tsn(ji,jj,ikt,1)
+          zts_top(ji,jj,2) = tsn(ji,jj,ikt,2)
         END DO
       END DO
-      CALL eos( ztstop, risfdep, zrhdtop_oce )
+      CALL eos( zts_top, risfdep, zrhdtop_oce )
 
 !==================================================================================     
@@ -667 +641 @@
          END DO
       END DO
-     !
-      CALL wrk_dealloc( jpi,jpj,2  , ztstop)
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj)
-      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce )
       !
    END SUBROUTINE hpg_isf
@@ -690 +660 @@
       REAL(wp) ::   z1_12, cffv, cffy   !    "         "
       LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )
-      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
-      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
-      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
-      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
+      !!----------------------------------------------------------------------
       !
       IF( ln_wd ) THEN
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1 
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 
+               ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
-             IF(ll_tmp1) THEN
-               zcpx(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
-             ELSE
-               zcpx(ji,jj) = 0._wp
-             END IF
+               IF(ll_tmp1) THEN
+                  zcpx(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
+               ELSE
+                  zcpx(ji,jj) = 0._wp
+               ENDIF
       
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+               ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
-             IF(ll_tmp1) THEN
-               zcpy(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
-             ELSE
-               zcpy(ji,jj) = 0._wp
-             END IF
-           END DO
-        END DO
-        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      END IF
+               IF(ll_tmp1) THEN
+                  zcpy(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
+               ELSE
+                  zcpy(ji,jj) = 0._wp
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -903 +868 @@
          END DO
       END DO
-      CALL lbc_lnk(rho_k,'W',1.)
-      CALL lbc_lnk(rho_i,'U',1.)
-      CALL lbc_lnk(rho_j,'V',1.)
+      CALL lbc_lnk( rho_k, 'W', 1. )
+      CALL lbc_lnk( rho_i, 'U', 1. )
+      CALL lbc_lnk( rho_j, 'V', 1. )
 
 
@@ -949 +914 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )
-      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
-      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
-      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd )   DEALLOCATE( zcpx, zcpy )
       !
    END SUBROUTINE hpg_djc
@@ -980 +942 @@
       REAL(wp) :: zrhdt1
       REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
-      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh )
-      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n )
-      IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
+      REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
+      !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN
@@ -1003 +960 @@
 
       IF( ln_wd ) THEN
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1 
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 
+               ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
-             IF(ll_tmp1) THEN
-               zcpx(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
-             ELSE
-               zcpx(ji,jj) = 0._wp
-             END IF
+               IF(ll_tmp1) THEN
+                  zcpx(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                             &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
+               ELSE
+                  zcpx(ji,jj) = 0._wp
+               ENDIF
       
-             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+               ll_tmp1 = MIN(   sshn(ji,jj)             ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
                   &                                                         > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
+               ll_tmp2 = ( ABS( sshn(ji,jj)             -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
                   &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
-             IF(ll_tmp1) THEN
-               zcpy(ji,jj) = 1.0_wp
-             ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+               IF(ll_tmp1) THEN
+                  zcpy(ji,jj) = 1.0_wp
+               ELSE IF(ll_tmp2) THEN
+                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
                            &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
-             ELSE
-               zcpy(ji,jj) = 0._wp
-             END IF
-           END DO
-        END DO
-        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
-      END IF
+               ELSE
+                  zcpy(ji,jj) = 0._wp
+               ENDIF
+            END DO
+         END DO
+         CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
 
       ! Clean 3-D work arrays
@@ -1298 +1256 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
-      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh )
-      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n )
-      IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
+      IF( ln_wd )   DEALLOCATE( zcpx, zcpy )
       !
    END SUBROUTINE hpg_prj
@@ -1353 +1308 @@
            !!Simply geometric average
                DO jk = 2, jpkm1-1
-                  zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))
-                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))
+                  zdf1 = (fsp(ji,jj,jk  ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk  ) - xsp(ji,jj,jk-1))
+                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk  )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk  ))
 
                   IF(zdf1 * zdf2 <= 0._wp) THEN
@@ -1403 +1358 @@
             END DO
          END DO
-
+         !
       ELSE
-           CALL ctl_stop( 'invalid polynomial type in cspline' )
-      ENDIF
-
+         CALL ctl_stop( 'invalid polynomial type in cspline' )
+      ENDIF
+      !
    END SUBROUTINE cspline
 

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

@@ -22 +22 @@
    USE lib_mpp         ! MPP library
    USE prtctl          ! Print control
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
    USE bdy_oce         ! ocean open boundary conditions
@@ -39 +38 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -75 +74 @@
       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zu, zv       ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
-      INTEGER  ::   jb                 ! dummy loop indices
-      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
-      INTEGER  ::   fu, fv
+      INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices
+      INTEGER  ::   ii, ifu, ib_bdy   ! local integers
+      INTEGER  ::   ij, ifv, igrd     !   -       -
+      REAL(wp) ::   zu, zv            ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
+      IF( ln_timing )   CALL timing_start('dyn_keg')
       !
       IF( kt == nit000 ) THEN
@@ -94 +90 @@
       ENDIF
 
-      IF( l_trddyn ) THEN           ! Save ua and va trends
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
+      IF( l_trddyn ) THEN           ! Save the input trends
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
@@ -112 +108 @@
                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-                     fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
-                     un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
+                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
+                     un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
                   END DO
                END DO
@@ -122 +118 @@
                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-                     fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
-                     vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
+                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
+                     vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
                   END DO
                END DO
@@ -172 +168 @@
       ENDIF      
 
-
       !
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
@@ -187 +182 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
+         DEALLOCATE( ztrdu , ztrdv )
       ENDIF
       !
@@ -193 +188 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
+      IF( ln_timing )   CALL timing_stop('dyn_keg')
       !
    END SUBROUTINE dyn_keg

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

@@ -27 +27 @@
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -48 +47 @@
 #  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)
@@ -62 +61 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf')
+      IF( ln_timing )   CALL timing_start('dyn_ldf')
       !
       IF( l_trddyn )   THEN                      ! temporary save of momentum trends
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
@@ -85 +84 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
+         DEALLOCATE ( ztrdu , ztrdv )
       ENDIF
       !                                          ! print sum trends (used for debugging)
@@ -91 +90 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf')
+      IF( ln_timing )   CALL timing_stop('dyn_ldf')
       !
    END SUBROUTINE dyn_ldf
@@ -102 +101 @@
       !! ** Purpose :   initializations of the horizontal ocean dynamics physics
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio, ierr         ! temporary integers 
+      INTEGER ::   ioptio, ierr   ! temporary integers 
       !!----------------------------------------------------------------------
       !
-      !                                   ! Namelist nam_dynldf: already read in ldfdyn module
+      !                                !==  Namelist nam_dynldf  ==!   already read in ldfdyn module
       !
-      IF(lwp) THEN                        ! Namelist print
+      IF(lwp) THEN                     !== Namelist print  ==!
          WRITE(numout,*)
          WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics'
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '      laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
-         WRITE(numout,*) '      bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
-         WRITE(numout,*) '      iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
-         WRITE(numout,*) '      horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
-         WRITE(numout,*) '      iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
+         WRITE(numout,*) '      Type of operator'
+         WRITE(numout,*) '         no explicit diffusion       ln_dynldf_NONE = ', ln_dynldf_NONE
+         WRITE(numout,*) '         laplacian operator          ln_dynldf_lap  = ', ln_dynldf_lap
+         WRITE(numout,*) '         bilaplacian operator        ln_dynldf_blp  = ', ln_dynldf_blp
+         WRITE(numout,*) '      Direction of action'
+         WRITE(numout,*) '         iso-level                   ln_dynldf_lev  = ', ln_dynldf_lev
+         WRITE(numout,*) '         horizontal (geopotential)   ln_dynldf_hor  = ', ln_dynldf_hor
+         WRITE(numout,*) '         iso-neutral                 ln_dynldf_iso  = ', ln_dynldf_iso
       ENDIF
-      !                                   ! use of lateral operator or not
+      !                                !==  use of lateral operator or not  ==!
       nldf = np_ERROR
       ioptio = 0
-      IF( ln_dynldf_lap )   ioptio = ioptio + 1
-      IF( ln_dynldf_blp )   ioptio = ioptio + 1
-      IF( ioptio >  1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on momentum' )
-      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral mixing operator
+      IF( ln_dynldf_NONE ) THEN   ;   nldf = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_dynldf_lap  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_dynldf_blp  ) THEN   ;                          ioptio = ioptio + 1   ;   ENDIF
+      IF( ioptio /= 1    )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
       !
-      IF( nldf /= np_no_ldf ) THEN        ! direction ==>> type of operator  
+      IF(.NOT.ln_dynldf_NONE ) THEN    !==  direction ==>> type of operator  ==!
          ioptio = 0
          IF( ln_dynldf_lev )   ioptio = ioptio + 1
          IF( ln_dynldf_hor )   ioptio = ioptio + 1
          IF( ln_dynldf_iso )   ioptio = ioptio + 1
-         IF( ioptio >  1   )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
-         IF( ioptio == 0   )   CALL ctl_stop( '          use at least ONE direction (level/hor/iso)' )
+         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' )
          !
-         !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
+         !                             ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals
          ierr = 0
-         IF ( ln_dynldf_lap ) THEN      ! laplacian operator
-            IF ( ln_zco ) THEN                ! z-coordinate
+         IF( ln_dynldf_lap ) THEN         ! laplacian operator
+            IF( ln_zco ) THEN                ! z-coordinate
                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation)
                IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level = horizontal (no rotation)
                IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation)
             ENDIF
-            IF ( ln_zps ) THEN             ! z-coordinate with partial step
+            IF( ln_zps ) THEN                ! z-coordinate with partial step
                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level              (no rotation)
                IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level              (no rotation)
                IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation)
             ENDIF
-            IF ( ln_sco ) THEN             ! s-coordinate
+            IF( ln_sco ) THEN                ! s-coordinate
                IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation)
                IF ( ln_dynldf_hor )   nldf = np_lap_i   ! horizontal             (   rotation)
@@ -154 +155 @@
          ENDIF
          !
-         IF( ln_dynldf_blp ) THEN          ! bilaplacian operator
-            IF ( ln_zco ) THEN                ! z-coordinate
-               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level = horizontal (no rotation)
-               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level = horizontal (no rotation)
-               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator
+            IF( ln_zco ) THEN                ! z-coordinate
+               IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_hor )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
             ENDIF
-            IF ( ln_zps ) THEN             ! z-coordinate with partial step
-               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
-               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level              (no rotation)
-               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+            IF( ln_zps ) THEN                ! z-coordinate with partial step
+               IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   nldf = np_blp     ! iso-level              (no rotation)
+               IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
             ENDIF
-            IF ( ln_sco ) THEN             ! s-coordinate
-               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
-               IF ( ln_dynldf_hor )   ierr = 2          ! horizontal             (   rotation)
-               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
+            IF( ln_sco ) THEN                ! s-coordinate
+               IF( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   ierr = 2          ! horizontal             (   rotation)
+               IF( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation)
             ENDIF
          ENDIF

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

@@ -28 +28 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
 
@@ -45 +44 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -108 +107 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
-      REAL(wp) ::   zmskt, zmskf                                     !   -      -
-      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
-      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
-      !
-      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
+      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
+      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
+      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4            !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso')
-      !
-      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
+      IF( ln_timing )   CALL timing_start('dyn_ldf_iso')
       !
       IF( kt == nit000 ) THEN
@@ -343 +339 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk)
+               zcof0 = 0.5_wp * rn_aht_0 * umask(ji,jj,jk)
                !
-               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
-               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
+               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
+               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
                !
-               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
-                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
-               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   &
-                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. )
-
-               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
-               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
+               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
+                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
+               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
+                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
+
+               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
+               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
                ! vertical flux on u field
-               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
-                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
-                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
-                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
+               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
+                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
+                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
+                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
                ! vertical mixing coefficient (akzu)
-               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
                akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
             END DO
@@ -369 +365 @@
          DO jk = 2, jpkm1
             DO ji = 2, jpim1
-               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk)
-
-               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
-               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
-
-               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
-                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
-               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
-                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
-
-               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
-               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
+               zcof0 = 0.5_wp * rn_aht_0 * vmask(ji,jj,jk)
+               !
+               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
+               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
+               !
+               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
+                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
+               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
+                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
+
+               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
+               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
                ! vertical flux on v field
-               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
-                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
-                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
-                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
+               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
+                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
+                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
+                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
                ! vertical mixing coefficient (akzv)
-               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
                akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
             END DO
@@ -404 +400 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso')
+      IF( ln_timing )   CALL timing_stop('dyn_ldf_iso')
       !
    END SUBROUTINE dyn_ldf_iso

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90

@@ -19 +19 @@
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -31 +30 @@
 #  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)
@@ -57 +56 @@
       REAL(wp) ::   zsign        ! local scalars
       REAL(wp) ::   zua, zva     ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zcur, zdiv
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
       !!----------------------------------------------------------------------
       !
@@ -66 +65 @@
       ENDIF
       !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_ldf_lap')
-      !
-      CALL wrk_alloc( jpi, jpj, zcur, zdiv ) 
+      IF( ln_timing )   CALL timing_start('dyn_ldf_lap')
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
@@ -107 +104 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zcur, zdiv ) 
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap')
+      IF( ln_timing )   CALL timing_stop('dyn_ldf_lap')
       !
    END SUBROUTINE dyn_ldf_lap
@@ -131 +127 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zulap, zvlap   ! laplacian at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_blp')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zulap, zvlap ) 
+      IF( ln_timing )   CALL timing_start('dyn_ldf_blp')
       !
       IF( kt == nit000 )  THEN
@@ -154 +148 @@
       CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zulap, zvlap ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_blp')
+      IF( ln_timing )   CALL timing_stop('dyn_ldf_blp')
       !
    END SUBROUTINE dyn_ldf_blp

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -44 +44 @@
    USE lbclnk         ! lateral boundary condition (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE prtctl         ! Print control
    USE timing         ! Timing
@@ -57 +56 @@
 
    !!----------------------------------------------------------------------
-   !! 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)
@@ -97 +96 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
-      !
-      IF( ln_dynspg_ts       )   CALL wrk_alloc( jpi,jpj,       zue, zve)
-      IF( l_trddyn           )   CALL wrk_alloc( jpi,jpj,jpk,   zua, zva)
+      IF( ln_timing    )   CALL timing_start('dyn_nxt')
+      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
+      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
       !
       IF( kt == nit000 ) THEN
@@ -253 +251 @@
             ELSE                          ! Asselin filter applied on thickness weighted velocity
                !
-               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
+               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
                ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
@@ -280 +278 @@
                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
                !
-               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
+               DEALLOCATE( ze3u_f , ze3v_f )
             ENDIF
             !
@@ -346 +344 @@
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
-      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
+      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
+      IF( l_trddyn     )   DEALLOCATE( zua, zva )
+      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
       !
    END SUBROUTINE dyn_nxt

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -28 +28 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -47 +46 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
+   !! NEMO/OPA 4.0 , LODYC-IPSL  (2017)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -71 +70 @@
       !!             period is used to prevent the divergence of odd and even time step.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
-      !
-      INTEGER  ::   ji, jj, jk                             ! dummy loop indices
-      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r             ! temporary scalar
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpice
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_spg')
+      INTEGER, INTENT(in   ) ::   kt   ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r   ! local scalars
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_spg')
       !
       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
-         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
@@ -124 +123 @@
          !
          IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
-            CALL wrk_alloc( jpi,jpj,   zpice )
-            !                                            
+            ALLOCATE( zpice(jpi,jpj) )
             zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
             zgrau0r     = - grav * r1_rau0
@@ -135 +133 @@
                END DO
             END DO
-            !
-            CALL wrk_dealloc( jpi,jpj,   zpice )         
+            DEALLOCATE( zpice )         
          ENDIF
          !
@@ -161 +158 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
+         DEALLOCATE( ztrdu , ztrdv ) 
       ENDIF
       !                                      ! print mean trends (used for debugging)
@@ -167 +164 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg')
+      IF( ln_timing )   CALL timing_stop('dyn_spg')
       !
    END SUBROUTINE dyn_spg
@@ -186 +183 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_init')
+      IF( ln_timing )   CALL timing_start('dyn_spg_init')
       !
       REWIND( numnam_ref )              ! Namelist namdyn_spg in reference namelist : Free surface
@@ -227 +224 @@
       ENDIF
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init')
+      IF( ln_timing )   CALL timing_stop('dyn_spg_init')
       !
    END SUBROUTINE dyn_spg_init

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

@@ -61 +61 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_exp')
+      IF( ln_timing )   CALL timing_start('dyn_spg_exp')
       !
       IF( kt == nit000 ) THEN
@@ -93 +93 @@
       ENDIF
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_exp')
+      IF( ln_timing )   CALL timing_stop('dyn_spg_exp')
       !
    END SUBROUTINE dyn_spg_exp

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -162 +162 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_spg_ts')
+      IF( ln_timing )   CALL timing_start('dyn_spg_ts')
       !
       IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
@@ -1125 +1125 @@
       IF( ln_wd )   DEALLOCATE( zcpx, zcpy )
       !
-      IF ( ln_diatmb ) THEN
+      IF( ln_diatmb ) THEN
          CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity
          CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity
       ENDIF
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
+      IF( ln_timing )   CALL timing_stop('dyn_spg_ts')
       !
    END SUBROUTINE dyn_spg_ts

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

@@ -14 +14 @@
    !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
    !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
-   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
-   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
-   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory
+   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
+   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
    !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
+   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   dyn_vor      : Update the momentum trend with the vorticity trend
-   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
-   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
-   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
-   !!   dyn_vor_init : set and control of the different vorticity option
+   !!   dyn_vor       : Update the momentum trend with the vorticity trend
+   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
+   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
+   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
+   !!   dyn_vor_init  : set and control of the different vorticity option
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
    USE dom_oce        ! ocean space and time domain
    USE dommsk         ! ocean mask
-   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
+   USE dynadv         ! momentum advection
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
@@ -40 +41 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-
 
    IMPLICIT NONE
@@ -80 +79 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -98 +97 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_vor')
-      !
-      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
-      !
-      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==!
-      !
-      CASE ( np_ENE )                                 !* energy conserving scheme
-         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_vor')
+      !
+      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
+         !
+         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
+         !
+         ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
+         ztrdv(:,:,:) = va(:,:,:)
+         SELECT CASE( nvor_scheme )
+         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme
+            IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+         CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme
+            IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme
+            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+         END SELECT
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
+         !
+         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, nrvm, un , vn , ua, va )                    ! relative vorticity or metric trend
+            SELECT CASE( nvor_scheme )
+            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme
+            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme
+            CASE( np_EEN )           ;   CALL vor_een( kt, nrvm, un , vn , ua, va )  ! energy & enstrophy scheme
+            END SELECT
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, un , vn , ua, va )                    ! planetary vorticity trend
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE                                               ! total vorticity trend
+         ENDIF
+         !
+         DEALLOCATE( ztrdu, ztrdv )
+         !
+      ELSE              !==  total vorticity trend added to the general trend  ==!
+         !
+         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
+         CASE( np_ENE )                        !* energy conserving scheme
                              CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend
             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         ENDIF
-         !
-      CASE ( np_ENS )                                 !* enstrophy conserving scheme
-         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two    
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE                                               ! total vorticity trend
+         CASE( np_ENS )                        !* enstrophy conserving scheme
                              CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend
             IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
-         ENDIF
-         !
-      CASE ( np_MIX )                                 !* mixed ene-ens scheme
-         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend (ens)
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend (ene)
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE                                               ! total vorticity trend
+         CASE( np_MIX )                        !* mixed ene-ens scheme
                              CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
                              CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-        ENDIF
-         !
-      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme
-         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
-         ELSE                                               ! total vorticity trend
+         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
                              CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend
             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         ENDIF
-         !
-      END SELECT
+         END SELECT
+         !
+      ENDIF
       !
       !                       ! print sum trends (used for debugging)
@@ -190 +159 @@
          &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor')
+      IF( ln_timing )   CALL timing_stop('dyn_vor')
       !
    END SUBROUTINE dyn_vor
@@ -217 +184 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
-      !                                                                ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
+      INTEGER                          , INTENT(in   )::   kt          ! ocean time-step index
+      INTEGER                          , INTENT(in   )::   kvor        ! total, planetary, relative, or metric
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
       REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('vor_ene')
-      !
-      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('vor_ene')
       !
       IF( kt == nit000 ) THEN
@@ -264 +228 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
+                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   &
                      &                   * r1_e1e2f(ji,jj)
                END DO
@@ -311 +275 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
+      !
+      IF( ln_timing )  CALL timing_stop('vor_ene')
       !
    END SUBROUTINE vor_ene
@@ -338 +301 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
-         !                                                             ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
+      INTEGER                          , INTENT(in   )::   kt          ! ocean time-step index
+      INTEGER                          , INTENT(in   )::   kvor        ! total, planetary, relative, or metric
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zuav, zvau   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
-      !
-      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('vor_ens')
       !
       IF( kt == nit000 ) THEN
@@ -431 +391 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
+      !
+      IF( ln_timing )   CALL timing_stop('vor_ens')
       !
    END SUBROUTINE vor_ens
@@ -455 +414 @@
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
-         !                                                             ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
+      INTEGER                          , INTENT(in   )::   kt          ! ocean time-step index
+      INTEGER                          , INTENT(in   )::   kvor        ! total, planetary, relative, or metric
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -465 +423 @@
       REAL(wp) ::   zua, zva     ! local scalars
       REAL(wp) ::   zmsk, ze3    ! local scalars
-      !
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f
-      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('vor_een')
-      !
-      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
-      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
+      REAL(wp), DIMENSION(jpi,jpj)   :: zwx , zwy , zwz , z1_e3f
+      REAL(wp), DIMENSION(jpi,jpj)   :: ztnw, ztne, ztsw, ztse
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('vor_een')
       !
       IF( kt == nit000 ) THEN
@@ -599 +553 @@
       !                                                ! ===============
       !
-      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
-      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
+      IF( ln_timing )   CALL timing_stop('vor_een')
       !
    END SUBROUTINE vor_een
@@ -618 +569 @@
       INTEGER ::   ios             ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
+      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix,   &
+         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_msk
       !!----------------------------------------------------------------------
 
@@ -672 +624 @@
       !                      
       IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
-      ncor = np_COR
-      IF( ln_dynadv_vec ) THEN     
-         IF(lwp) WRITE(numout,*) '      ===>>   Vector form advection : vorticity = Coriolis + relative vorticity'
+      ncor = np_COR                       ! planetary vorticity
+      SELECT CASE( n_dynadv )
+      CASE( np_LIN_dyn )
+         IF(lwp) WRITE(numout,*) '      ===>>   linear dynamics : total vorticity = Coriolis'
+         nrvm = np_COR        ! planetary vorticity
+         ntot = np_COR        !     -         -
+      CASE( np_VEC_c2  )
+         IF(lwp) WRITE(numout,*) '      ===>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
          nrvm = np_RVO        ! relative vorticity
-         ntot = np_CRV        ! relative + planetary vorticity
-      ELSE                        
-         IF(lwp) WRITE(numout,*) '      ===>>   Flux form advection   : vorticity = Coriolis + metric term'
+         ntot = np_CRV        ! relative + planetary vorticity         
+      CASE( np_FLX_c2 , np_FLX_ubs  )
+         IF(lwp) WRITE(numout,*) '      ===>>   flux form dynamics : total vorticity = Coriolis + metric term'
          nrvm = np_MET        ! metric term
          ntot = np_CME        ! Coriolis + metric term
-      ENDIF
+      END SELECT
       
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '      ===>>   energy conserving scheme'
-         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '      ===>>   enstrophy conserving scheme'
-         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '      ===>>   mixed enstrophy/energy conserving scheme'
-         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '      ===>>   energy and enstrophy conserving scheme'
+         SELECT CASE( nvor_scheme )
+         CASE( np_ENE )   ;   WRITE(numout,*) '      ===>>   energy conserving scheme'
+         CASE( np_ENS )   ;   WRITE(numout,*) '      ===>>   enstrophy conserving scheme'
+         CASE( np_MIX )   ;   WRITE(numout,*) '      ===>>   mixed enstrophy/energy conserving scheme'
+         CASE( np_EEN )   ;   WRITE(numout,*) '      ===>>   energy and enstrophy conserving scheme'
+         END SELECT         
       ENDIF
       !

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

@@ -5 +5 @@
    !!======================================================================
    !! History :  OPA  ! 1991-01  (G. Madec) Original code
-   !!            7.0  ! 1991-11  (G. Madec)
-   !!            7.5  ! 1996-01  (G. Madec) statement function for e3
    !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
    !!----------------------------------------------------------------------
@@ -22 +20 @@
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -29 +26 @@
    
    PUBLIC   dyn_zad       ! routine called by dynadv.F90
-   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
 
    !! * Substitutions
 #  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)
@@ -58 +54 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
       !
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      REAL(wp) ::   zua, zva        ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zua, zva     ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zww
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuw, zwvw
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_zad')
-      !
-      CALL wrk_alloc( jpi,jpj, zww ) 
-      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 
+      IF( ln_timing )   CALL timing_start('dyn_zad')
       !
       IF( kt == nit000 ) THEN
-         IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
       ENDIF
 
       IF( l_trddyn )   THEN         ! Save ua and va trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
@@ -96 +89 @@
       !
       ! Surface and bottom advective fluxes set to zero
-      IF ( ln_isfcav ) THEN
+      IF( ln_isfcav ) THEN
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1           ! vector opt.
@@ -119 +112 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1       ! vector opt.
-               !                         ! vertical momentum advective trends
-               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
-               !                         ! add the trends to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zva
+               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
             END DO  
          END DO  
@@ -133 +122 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+         DEALLOCATE( ztrdu, ztrdv ) 
       ENDIF
       !                             ! Control print
@@ -139 +128 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj, zww ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
+      IF( ln_timing )   CALL timing_stop('dyn_zad')
       !
    END SUBROUTINE dyn_zad
 
-
-   SUBROUTINE dyn_zad_zts ( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dynzad_zts  ***
-      !! 
-      !! ** Purpose :   Compute the now vertical momentum advection trend and 
-      !!      add it to the general trend of momentum equation. This version
-      !!      uses sub-timesteps for improved numerical stability with small
-      !!      vertical grid sizes. This is especially relevant when using 
-      !!      embedded ice with thin surface boxes.
-      !!
-      !! ** Method  :   The now vertical advection of momentum is given by:
-      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
-      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
-      !!      Add this trend to the general trend (ua,va):
-      !!         (ua,va) = (ua,va) + w dz(u,v)
-      !!
-      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
-      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
-      !
-      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
-      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
-      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
-      REAL(wp) ::   zua, zva        ! temporary scalars
-      REAL(wp) ::   zr_rdt          ! temporary scalar
-      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
-      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
-      CALL wrk_alloc( jpi,jpj,jpk,3,   zus , zvs ) 
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
-      ENDIF
-
-      IF( l_trddyn )   THEN         ! Save ua and va trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
-         ztrdu(:,:,:) = ua(:,:,:) 
-         ztrdv(:,:,:) = va(:,:,:) 
-      ENDIF
-      
-      IF( neuler == 0 .AND. kt == nit000 ) THEN
-          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
-      ELSE
-          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
-      ENDIF
-      
-      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
-         DO jj = 2, jpj                   
-            DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
-         DO ji = fs_2, fs_jpim1           ! vector opt.
- !!gm missing ISF boundary condition
-            zwuw(ji,jj, 1 ) = 0._wp
-            zwvw(ji,jj, 1 ) = 0._wp
-            zwuw(ji,jj,jpk) = 0._wp
-            zwvw(ji,jj,jpk) = 0._wp
-         END DO  
-      END DO
-
-! Start with before values and use sub timestepping to reach after values
-
-      zus(:,:,:,1) = ub(:,:,:)
-      zvs(:,:,:,1) = vb(:,:,:)
-
-      DO jl = 1, jnzts                   ! Start of sub timestepping loop
-
-         IF( jl == 1 ) THEN              ! Euler forward to kick things off
-           jtb = 1   ;   jtn = 1   ;   jta = 2
-           zts = z2dtzts
-         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
-           jtb = 1   ;   jtn = 2   ;   jta = 3
-           zts = 2._wp * z2dtzts
-         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           ! Vertical momentum advection at level w and u- and v- vertical
-            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
-               DO ji = fs_2, fs_jpim1        ! vector opt.
-                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)
-                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)
-               END DO  
-            END DO   
-         END DO
-         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1       ! vector opt.
-                  !                         ! vertical momentum advective trends
-                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
-                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
-                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
-               END DO  
-            END DO  
-         END DO
-
-      END DO      ! End of sub timestepping loop
-
-      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
-      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1       ! vector opt.
-               !                         ! vertical momentum advective trends
-               !                         ! add the trends to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
-               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
-            END DO  
-         END DO  
-      END DO
-
-      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
-      ENDIF
-      !                             ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
-      CALL wrk_dealloc( jpi,jpj,jpk,3,   zus , zvs ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
-      !
-   END SUBROUTINE dyn_zad_zts
-
    !!======================================================================
 END MODULE dynzad

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

@@ -76 +76 @@
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf')
+      IF( ln_timing )   CALL timing_start('dyn_zdf')
       !
       IF( kt == nit000 ) THEN       !* initialization
@@ -392 +392 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
          !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf')
+      IF( ln_timing )   CALL timing_stop('dyn_zdf')
       !
    END SUBROUTINE dyn_zdf

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -22 +22 @@
    USE divhor         ! horizontal divergence
    USE phycst         ! physical constants
-   USE bdy_oce   , ONLY: ln_bdy, bdytmask
+   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
    USE bdydyn2d       ! bdy_ssh routine
 #if defined key_agrif
@@ -36 +36 @@
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE wet_dry         ! Wetting/Drying flux limting
+   USE wet_dry        ! Wetting/Drying flux limting
 
    IMPLICIT NONE
@@ -74 +73 @@
       INTEGER  ::   jk            ! dummy loop indice
       REAL(wp) ::   z2dt, zcoef   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
-      !
-      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('ssh_nxt')
       !
       IF( kt == nit000 ) THEN
@@ -134 +131 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
       !
-      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
+      IF( ln_timing )   CALL timing_stop('ssh_nxt')
       !
    END SUBROUTINE ssh_nxt
@@ -160 +155 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   z1_2dt       ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('wzv')
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
@@ -180 +174 @@
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
-         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
+         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
          !
          DO jk = 1, jpkm1
@@ -200 +194 @@
          END DO
          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
-         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
+         DEALLOCATE( zhdiv ) 
       ELSE   ! z_star and linear free surface cases
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
@@ -215 +209 @@
       ENDIF
       !
-      IF( nn_timing == 1 )  CALL timing_stop('wzv')
+      IF( ln_timing )   CALL timing_stop('wzv')
       !
    END SUBROUTINE wzv
@@ -244 +238 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
+      IF( ln_timing )  CALL timing_start('ssh_swp')
       !
       IF( kt == nit000 ) THEN
@@ -271 +265 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-      IF( nn_timing == 1 )   CALL timing_stop('ssh_swp')
+      IF( ln_timing )   CALL timing_stop('ssh_swp')
       !
    END SUBROUTINE ssh_swp

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

@@ -11 +11 @@
 
    !!----------------------------------------------------------------------
-   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
-   !!                when wetting and drying happens 
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
-   USE sbcrnf          ! river runoff 
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   !!   wad_init      : initialisation of wetting and drying
+   !!   wad_lmt       : horizontal flux limiter and limited velocity when wetting and drying happens
+   !!   wad_lmt_bt    : same as wad_lmt for the barotropic stepping (dynspg_ts)
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce  , ONLY: ln_rnf   ! surface boundary condition: ocean
+   USE sbcrnf         ! river runoff 
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -31 +32 @@
    !! ---------------------------------------------------------------------
 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd          !: wetting and drying t-pt depths
-                                                                     !  (can include negative depths)
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask   !: u- and v- limiter 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd    !: wetting and drying t-pt depths
+   !                                                           !  (can include negative depths)
 
    LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
    REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
    REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
-   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying 
-                                           !: will be considered
+   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying will be considered
    INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
 
@@ -48 +48 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -58 +59 @@
       !! ** input   : - namwad namelist
       !!----------------------------------------------------------------------
+      INTEGER  ::   ios, ierr   ! Local integer
+      !!
       NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      INTEGER  ::   ierr                ! Local integer status array allocation 
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namwad in reference namelist 
-                                        ! : Parameters for Wetting/Drying
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying
       READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
 905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
-      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist 
-                                        ! : Parameters for Wetting/Drying
+      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying
       READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
       IF(lwm) WRITE ( numond, namwad )
-
+      !
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
@@ -103 +102 @@
       !! ** Action  : - calculate flux limiter and W/D flag
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
-      REAL(wp), INTENT(in) :: z2dt
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1        !!gm DOCTOR names: should start with p !
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sshemp
+      REAL(wp)                , INTENT(in   ) ::   z2dt
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
@@ -113 +112 @@
       REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
       REAL(wp) ::   ztmp                ! local scalars
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
-      !!----------------------------------------------------------------------
-      !
-
-      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
-
-      IF(ln_wd) THEN
-
-        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
-        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
-        !
-       
-        !IF(lwp) WRITE(numout,*)
-        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
-       
-        jflag  = 0
-        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
-
-       
-        zflxp(:,:)   = 0._wp
-        zflxn(:,:)   = 0._wp
-        zflxu(:,:)   = 0._wp
-        zflxv(:,:)   = 0._wp
-
-        zwdlmtu(:,:)  = 1._wp
-        zwdlmtv(:,:)  = 1._wp
-       
-        ! Horizontal Flux in u and v direction
-        DO jk = 1, jpkm1  
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
-              END DO  
-           END DO  
-        END DO
-       
-        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
-        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
-       
-        wdmask(:,:) = 1
-        DO jj = 2, jpj
-           DO ji = 2, jpi 
-
-             IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE   ! we don't care about land cells
-             IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
-
-              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
-                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
-              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
-                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
-
-              zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
-              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
-                sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
-                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
-                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
-                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
-                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
-                wdmask(ji,jj) = 0._wp
-              END IF
-           ENDDO
-        END DO
-
-      
-        !! start limiter iterations 
-        DO jk1 = 1, nn_wdit + 1
-       
-          
-           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
-           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
-           jflag = 0     ! flag indicating if any further iterations are needed
-          
-           DO jj = 2, jpj
-              DO ji = 2, jpi 
-        
-                 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 
-                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
-        
-                 ztmp = e1e2t(ji,jj)
-
-                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
-                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
-                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
-                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
-          
-                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
-                 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
-          
-                 IF( zdep1 > zdep2 ) THEN
-                   wdmask(ji, jj) = 0
-                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
-                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
-                   ! flag if the limiter has been used but stop flagging if the only
-                   ! changes have zeroed the coefficient since further iterations will
-                   ! not change anything
-                   IF( zcoef > 0._wp ) THEN
-                      jflag = 1 
-                   ELSE
-                      zcoef = 0._wp
-                   ENDIF
-                   IF(jk1 > nn_wdit) zcoef = 0._wp
-                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
-                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
-                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
-                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
-                 END IF
-              END DO ! ji loop
-           END DO  ! jj loop
-
-           CALL lbc_lnk( zwdlmtu, 'U', 1. )
-           CALL lbc_lnk( zwdlmtv, 'V', 1. )
-
-           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
-
-           IF(jflag == 0) EXIT
-          
-        END DO  ! jk1 loop
-       
-        DO jk = 1, jpkm1
-          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
-          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
-        END DO
-
-        CALL lbc_lnk( un, 'U', -1. )
-        CALL lbc_lnk( vn, 'V', -1. )
-      !
-        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
-        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
-        CALL lbc_lnk( un_b, 'U', -1. )
-        CALL lbc_lnk( vn_b, 'V', -1. )
-       
-        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
-       
-        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
-        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
-        !
-        !
-        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
-        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
-        !
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+      REAL(wp),  DIMENSION(jpi,jpj) ::   zwdlmtu, zwdlmtv   ! W/D flux limiters
+      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxp  ,  zflxn    ! local 2D workspace
+      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxu  ,  zflxv    ! local 2D workspace
+      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxu1 , zflxv1    ! local 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('wad_lmt')
+      !
+      !IF(lwp) WRITE(numout,*)
+      !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
+      !
+      jflag  = 0
+      zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
+      ! 
+      zflxp(:,:)   = 0._wp
+      zflxn(:,:)   = 0._wp
+      zflxu(:,:)   = 0._wp
+      zflxv(:,:)   = 0._wp
+      !
+      zwdlmtu(:,:) = 1._wp
+      zwdlmtv(:,:) = 1._wp
+      ! 
+      ! Horizontal Flux in u and v direction
+      DO jk = 1, jpkm1  
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+            END DO  
+         END DO  
+      END DO
+      !
+      zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+      zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+      ! 
+      wdmask(:,:) = 1
+      DO jj = 2, jpj
+         DO ji = 2, jpi 
+            !
+            IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE   ! we don't care about land cells
+            IF( ht_wd(ji,jj)     > zdepwd )   CYCLE   ! and cells which are unlikely to dry
+            !
+            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )   &
+               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp ) 
+            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )   &
+               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp ) 
+            !
+            zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+            IF( zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
+               sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
+               IF( zflxu(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = 0._wp
+               IF( zflxu(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = 0._wp
+               IF( zflxv(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = 0._wp
+               IF( zflxv(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = 0._wp 
+               wdmask(ji,jj) = 0._wp
+            ENDIF
+         END DO
+      END DO
+      !!
+      !! start limiter iterations 
+      DO jk1 = 1, nn_wdit + 1
+         ! 
+         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+         jflag = 0     ! flag indicating if any further iterations are needed
+         ! 
+         DO jj = 2, jpj
+            DO ji = 2, jpi 
+               !
+               IF( tmask(ji,jj,1) < 0.5_wp )   CYCLE 
+               IF( ht_wd(ji,jj)   > zdepwd )   CYCLE
+               !
+               ztmp = e1e2t(ji,jj)
+               !
+               zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )   &
+                  &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp ) 
+               zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )   &
+                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp ) 
+               !
+               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+               zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
+               !
+               IF( zdep1 > zdep2 ) THEN
+                  wdmask(ji, jj) = 0
+                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                  ! flag if the limiter has been used but stop flagging if the only
+                  ! changes have zeroed the coefficient since further iterations will
+                  ! not change anything
+                  IF( zcoef > 0._wp ) THEN   ;   jflag = 1 
+                  ELSE                       ;   zcoef = 0._wp
+                  ENDIF
+                  IF( jk1 > nn_wdit )   zcoef = 0._wp
+                  IF( zflxu1(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = zcoef
+                  IF( zflxu1(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = zcoef
+                  IF( zflxv1(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = zcoef
+                  IF( zflxv1(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = zcoef
+               ENDIF
+            END DO ! ji loop
+         END DO  ! jj loop
+         !
+         CALL lbc_lnk( zwdlmtu, 'U', 1. )
+         CALL lbc_lnk( zwdlmtv, 'V', 1. )
+         !
+         IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain
+         !
+         IF(jflag == 0)   EXIT
+         ! 
+      END DO  ! jk1 loop
+       
+      DO jk = 1, jpkm1
+         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:) 
+         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:) 
+      END DO
+
+!!gm ==> Andrew  : the lbclnk below is useless since above lbclnk is applied on zwdlmtu/v
+!!                                             and un, vn always with lbclnk
+      CALL lbc_lnk( un, 'U', -1. )
+      CALL lbc_lnk( vn, 'V', -1. )
+!!gm end
+      !
+      un_b(:,:) = un_b(:,:) * zwdlmtu(:,:)
+      vn_b(:,:) = vn_b(:,:) * zwdlmtv(:,:)
+!!gm ==> Andrew   : probably same as above
+      CALL lbc_lnk( un_b, 'U', -1. )
+      CALL lbc_lnk( vn_b, 'V', -1. )
+!!gm end
+       
+      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
+       
+      !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+      !
+      !
+      !
+      IF( ln_timing )   CALL timing_stop('wad_lmt')
       !
    END SUBROUTINE wad_lmt
@@ -284 +272 @@
       REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
       REAL(wp) ::   ztmp                ! local scalars
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
-      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
-
-      IF(ln_wd) THEN
-
-        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
-        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
-        !
-       
-        !IF(lwp) WRITE(numout,*)
-        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
-       
-        jflag  = 0
-        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
-
-        z2dt = rdtbt   
-       
-        zflxp(:,:)   = 0._wp
-        zflxn(:,:)   = 0._wp
-
-        zwdlmtu(:,:)  = 1._wp
-        zwdlmtv(:,:)  = 1._wp
-       
-        ! Horizontal Flux in u and v direction
-       
-        DO jj = 2, jpj
-           DO ji = 2, jpi 
-
-             IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
-             IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
-
-              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
-                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
-              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
-                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
-
-              zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
-              IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
-                sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
-                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
-                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
-                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
-                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
-              END IF
-           ENDDO
-        END DO
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
+      REAL(wp), DIMENSION(jpi,jpj) ::   zflxp,  zflxn            ! local 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zflxu1, zflxv1           ! local 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('wad_lmt_bt')
+      !       
+      !IF(lwp) WRITE(numout,*)
+      !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
+       
+      jflag  = 0
+      zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
+
+      z2dt = rdtbt   
+       
+      zflxp(:,:)   = 0._wp
+      zflxn(:,:)   = 0._wp
+
+      zwdlmtu(:,:) = 1._wp
+      zwdlmtv(:,:) = 1._wp
+       
+      ! Horizontal Flux in u and v direction
+       
+      DO jj = 2, jpj
+         DO ji = 2, jpi 
+            !
+            IF( tmask(ji,jj,1) < 0.5_wp )   CYCLE   ! we don't care about land cells
+            IF( ht_wd(ji,jj)   > zdepwd )   CYCLE   ! and cells which are unlikely to dry
+            !
+            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )   &
+               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp ) 
+            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )   &
+               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp ) 
+
+            zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
+            IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
+               sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
+               IF( zflxu(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = 0._wp
+               IF( zflxu(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = 0._wp
+               IF( zflxv(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = 0._wp
+               IF( zflxv(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = 0._wp 
+            ENDIF
+         END DO
+      END DO
 
       
-        !! start limiter iterations 
-        DO jk1 = 1, nn_wdit + 1
-       
+      !! start limiter iterations 
+      DO jk1 = 1, nn_wdit + 1
+         ! 
+         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+         jflag = 0     ! flag indicating if any further iterations are needed
+         !
+         DO jj = 2, jpj
+            DO ji = 2, jpi 
+               !
+               IF( tmask(ji,jj, 1 ) < 0.5_wp  )   CYCLE 
+               IF( ht_wd(ji,jj)      > zdepwd )   CYCLE
+               !
+               ztmp = e1e2t(ji,jj)
+               !
+               zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )   &
+                  &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp ) 
+               zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )   &
+                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp ) 
           
-           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
-           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
-           jflag = 0     ! flag indicating if any further iterations are needed
+               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+               zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
           
-           DO jj = 2, jpj
-              DO ji = 2, jpi 
-        
-                 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 
-                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
-        
-                 ztmp = e1e2t(ji,jj)
-
-                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
-                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
-                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
-                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
-          
-                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
-                 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
-          
-                 IF(zdep1 > zdep2) THEN
-                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
-                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
-                   ! flag if the limiter has been used but stop flagging if the only
-                   ! changes have zeroed the coefficient since further iterations will
-                   ! not change anything
-                   IF( zcoef > 0._wp ) THEN
+               IF(zdep1 > zdep2) THEN
+                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                  ! flag if the limiter has been used but stop flagging if the only
+                  ! changes have zeroed the coefficient since further iterations will
+                  ! not change anything
+                  IF( zcoef > 0._wp ) THEN
                       jflag = 1 
-                   ELSE
+                  ELSE
                       zcoef = 0._wp
-                   ENDIF
-                   IF(jk1 > nn_wdit) zcoef = 0._wp
-                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
-                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
-                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
-                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
-                 END IF
-              END DO ! ji loop
-           END DO  ! jj loop
-
-           CALL lbc_lnk( zwdlmtu, 'U', 1. )
-           CALL lbc_lnk( zwdlmtv, 'V', 1. )
-
-           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
-
-           IF(jflag == 0) EXIT
-          
-        END DO  ! jk1 loop
-       
-        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
-        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
-
-        CALL lbc_lnk( zflxu, 'U', -1. )
-        CALL lbc_lnk( zflxv, 'V', -1. )
-       
-        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
-       
-        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
-        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
-        !
-        !
-        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
-        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
-        !
-      END IF
-
-      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+                  ENDIF
+                  IF( jk1 > nn_wdit )   zcoef = 0._wp
+                  IF( zflxu1(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = zcoef
+                  IF( zflxu1(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = zcoef
+                  IF( zflxv1(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = zcoef
+                  IF( zflxv1(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = zcoef
+               ENDIF
+            END DO ! ji loop
+         END DO  ! jj loop
+         !
+         CALL lbc_lnk( zwdlmtu, 'U', 1. )
+         CALL lbc_lnk( zwdlmtv, 'V', 1. )
+         !
+         IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain
+         !
+         IF( jflag == 0 )   EXIT
+         !    
+      END DO  ! jk1 loop
+      ! 
+      zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
+      zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
+      !
+      CALL lbc_lnk( zflxu, 'U', -1. )
+      CALL lbc_lnk( zflxv, 'V', -1. )
+      !
+      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
+       
+      !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+      !
+      IF( ln_timing )  CALL timing_stop('wad_lmt')
+      !
    END SUBROUTINE wad_lmt_bt